SRA meta data parsing

python으로 전체 SRA data중 원하는 데이터 골라내기

SRA가 뭔지에 대한 설명은 여기 참고

여기에선 SRA Metadata 를 활용해서 원하는 Raw data를 찾는 방법을 소개한다.

🤨⁉️ 지금까지 publish 된 데이터 중에서 single cell관련된 데이터 만을 모으고 싶다.

SRA metadata download

https://ftp-trace.ncbi.nlm.nih.gov/sra/reports/Metadata/ 이쪽으로 가면 지금까지 나온 전체 SRA accession에 해당하는 metadata를 다운 받을 수 있다. 가장 최신&Full버전으로 다운 받아주면 된다.

wget https://ftp-trace.ncbi.nlm.nih.gov/sra/reports/Metadata/NCBI_SRA_Metadata_Full_20240520.tar.gz

ls -l Metadata | grep "SRA" | grep "^d" > SRAaccession_list.txt
sed 's/ /\t/g' SRAaccession_list.txt | cut -f 9 > SRAaccession.txt

Metadata는 xml형식으로 작성되어 있다. 작성일 기준 20210715.tar.gz 파일에는 1,053,754개의 accession number가 있었다.

하나의 SRA number가 폴더형식으로 존재하고 그 아래 experiment.xml이나 study.xml 파일이 들어있는 형태. extperiment, study 둘다 존재 할 수도 있고 둘중 하나만 있을 수도 있다.

일단은 xml 파일 형태로는 parsing이 어려우니 table형태로 재가공 한 후, 원하는 단어가 들어 있는 SRA accession만 추출해내는 방식으로 진행했다.

Converting XML files to Table form

• experiment.xml 으로 정리

import os, glob, datetime

#SRA 284만개
main_dir = "/storage2/Project/subin/Project06_scNDB/01.Data/SRAdb/XML"

SRA_accession_list = open("%s/SRAaccession.txt"%main_dir)

SRA_parsing_out = open('%s/SRA_exp_accession_xml_parsing.txt'%main_dir, 'w')
SRA_no_exp_out = open('%s/SRA_no_exp.txt'%main_dir, 'w')


for  id_list in SRA_accession_list:
	accession_id = id_list.split('\n')[0]
	project_accession = ""
	exp_accession = ""
	sample_accession = ""
	title = ""
	lib_instrument = ""
	lib_strategy = ""
	lib_source = ""
	lib_protocol = ""
	try: 
		exp_xml = open("%s/Metadata/%s/%s.experiment.xml"%(main_dir,accession_id,accession_id))
		for line in exp_xml:
			if '<PRIMARY_ID>SRP' in line :
				project_accession = line.split('<PRIMARY_ID>')[1].split('</PRIMARY_ID>')[0]
				continue
			elif '<PRIMARY_ID>SRX' in line :
				exp_accession = line.split('<PRIMARY_ID>')[1].split('</PRIMARY_ID>')[0]
				continue
			elif '<PRIMARY_ID>SRS' in line :
				sample_accession = line.split('<PRIMARY_ID>')[1].split('</PRIMARY_ID>')[0]
				continue
			elif '<TITLE>' in line :
				title = line.split('<TITLE>')[1].split('</TITLE>')[0]
				continue
			elif '<INSTRUMENT_MODEL>' in line :
				lib_instrument = line.split('<INSTRUMENT_MODEL>')[1].split('</INSTRUMENT_MODEL>')[0]
				continue
			elif '<LIBRARY_STRATEGY>' in line :
				lib_strategy = line.split('<LIBRARY_STRATEGY>')[1].split('</LIBRARY_STRATEGY>')[0]
				continue
			elif '<LIBRARY_SOURCE>' in line :
				lib_source = line.split('<LIBRARY_SOURCE>')[1].split('</LIBRARY_SOURCE>')[0]
				continue
			elif '<LIBRARY_CONSTRUCTION_PROTOCOL>' in line :
				lib_protocol = line.split('<LIBRARY_CONSTRUCTION_PROTOCOL>')[1].split('</LIBRARY_CONSTRUCTION_PROTOCOL>')[0]
				continue

		SRA_parsing_out.write("%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n"%(accession_id,project_accession,exp_accession,sample_accession,title,lib_instrument,lib_strategy,lib_source,lib_protocol))
	except:
		SRA_no_exp_out.write("%s\n"%accession_id)

SRA_parsing_out.close()
SRA_no_exp_out.close()

• study.xml 로 정리

import os, glob, datetime

#SRA 284만개
main_dir = "/storage2/Project/subin/Project06_scNDB/01.Data/SRAdb/XML"

SRA_accession_list = open("%s/SRAaccession.txt"%main_dir)

SRA_study_parsing_out = open('%s/SRA_study_accession_xml_parsing.txt'%main_dir, 'w')
SRA_no_study_out = open('%s/SRA_no_study.txt'%main_dir, 'w')

for  id_list in SRA_accession_list:
	accession_id = id_list.split('\n')[0]
	project_accession = ""
	Bioproject_accession = ""
	GEO_accession = ""
	Study_title = ""
	Study_abstract = ""
	try: 
		study_xml = open("%s/Metadata/%s/%s.study.xml"%(main_dir,accession_id,accession_id))
		for line in study_xml:
			if '<PRIMARY_ID>SRP' in line :
				project_accession = line.split('<PRIMARY_ID>')[1].split('</PRIMARY_ID>')[0]
				continue
			elif '<EXTERNAL_ID namespace="BioProject" label="primary">' in line :
				Bioproject_accession = line.split('<EXTERNAL_ID namespace="BioProject" label="primary">')[1].split('</EXTERNAL_ID>')[0]
				continue
			elif '<EXTERNAL_ID namespace="GEO">' in line :
				GEO_accession = line.split('<EXTERNAL_ID namespace="GEO">')[1].split('</EXTERNAL_ID>')[0]
				continue
			elif '<STUDY_TITLE>' in line :
				Study_title = line.split('<STUDY_TITLE>')[1].split('</STUDY_TITLE>')[0]
				continue
			elif '<STUDY_ABSTRACT>' in line :
				Study_abstract = line.split('<STUDY_ABSTRACT>')[1].split('</STUDY_ABSTRACT>')[0]
				continue
		SRA_study_parsing_out.write("%s\t%s\t%s\t%s\t%s\t%s\t\n"%(accession_id,project_accession,Bioproject_accession,GEO_accession,Study_title,Study_abstract))
	except:
		SRA_no_study_out.write("%s\n"%accession_id)

SRA_study_parsing_out.close()
SRA_no_study_out.close()