XML NCBI Blast 파일에서 첫 번째 적중 요소를 추출하는 방법은 무엇입니까?
https://stackoverflow.com/questions/1935385
-
20-09-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
Russian문제
나는 NCBI XML Blast 파일에서 첫 번째 히트 만 추출하려고합니다. 다음으로 첫 번째 HSP 만 받고 싶습니다. 마지막 단계에서 나는 Best Score를 기준으로 이것을 얻고 싶습니다. 여기에서 XML 파일의 샘플을 명확하게하기 위해 :
<?xml version="1.0"?>
<!DOCTYPE BlastOutput PUBLIC "-//NCBI//NCBI BlastOutput/EN" "http://www.ncbi.nlm.nih.gov/dtd/NCBI_BlastOutput.dtd">
<BlastOutput>
<BlastOutput_program>blastx</BlastOutput_program>
<BlastOutput_version>blastx 2.2.22 [Sep-27-2009]</BlastOutput_version>
<BlastOutput_reference>~Reference: Altschul, Stephen F., Thomas L. Madden, Alejandro A. Schaffer, ~Jinghui Zhang, Zheng Zhang, Webb Miller, and David J. Lipman (1997), ~"Gapped BLAST and PSI-BLAST: a new generation of protein database search~programs", Nucleic Acids Res. 25:3389-3402.</BlastOutput_reference>
<BlastOutput_db>/Applications/blast/db/viral1.protein.faa</BlastOutput_db>
<BlastOutput_query-ID>lcl|1_0</BlastOutput_query-ID>
<BlastOutput_query-def>DSAD-090629_plate11A01a.g1 CHROMAT_FILE: DSAD-090629_plate11A01a.g1 PHD_FILE: DSAD-090629_plate11A01a.g1.phd.1 CHEM: term DYE: big TIME: Thu Sep 17 15:33:59 2009 TEMPLATE: DSAD-090629_plate11A01a DIRECTION: rev</BlastOutput_query-def>
<BlastOutput_query-len>1024</BlastOutput_query-len>
<BlastOutput_param>
<Parameters>
<Parameters_matrix>BLOSUM62</Parameters_matrix>
<Parameters_expect>1e-05</Parameters_expect>
<Parameters_gap-open>11</Parameters_gap-open>
<Parameters_gap-extend>1</Parameters_gap-extend>
<Parameters_filter>F</Parameters_filter>
</Parameters>
</BlastOutput_param>
<BlastOutput_iterations>
<Iteration>
<Iteration_iter-num>1</Iteration_iter-num>
<Iteration_query-ID>lcl|1_0</Iteration_query-ID>
<Iteration_query-def>DSAD-090629_plate11A01a.g1 CHROMAT_FILE: DSAD-090629_plate11A01a.g1 PHD_FILE: DSAD-090629_plate11A01a.g1.phd.1 CHEM: term DYE: big TIME: Thu Sep 17 15:33:59 2009 TEMPLATE: DSAD-090629_plate11A01a DIRECTION: rev</Iteration_query-def>
<Iteration_query-len>1024</Iteration_query-len>
<Iteration_stat>
<Statistics>
<Statistics_db-num>68007</Statistics_db-num>
<Statistics_db-len>19518578</Statistics_db-len>
<Statistics_hsp-len>0</Statistics_hsp-len>
<Statistics_eff-space>0</Statistics_eff-space>
<Statistics_kappa>0.041</Statistics_kappa>
<Statistics_lambda>0.267</Statistics_lambda>
<Statistics_entropy>0.14</Statistics_entropy>
</Statistics>
</Iteration_stat>
<Iteration_message>No hits found</Iteration_message>
</Iteration>
<Iteration>
<Iteration>
<Iteration_iter-num>6</Iteration_iter-num>
<Iteration_query-ID>lcl|6_0</Iteration_query-ID>
<Iteration_query-def>DSAD-090629_plate11A05a.g1 CHROMAT_FILE: DSAD-090629_plate11A05a.g1 PHD_FILE: DSAD-090629_plate11A05a.g1.phd.1 CHEM: term DYE: big TIME: Thu Sep 17 15:33:59 2009 TEMPLATE: DSAD-090629_plate11A05a DIRECTION: rev</Iteration_query-def>
<Iteration_query-len>1068</Iteration_query-len>
<Iteration_hits>
<Hit>
<Hit_num>1</Hit_num>
<Hit_id>gnl|BL_ORD_ID|23609</Hit_id>
<Hit_def>gi|38707884|ref|NP_945016.1| Putative ribose-phosphate pyrophosphokinase [Enterobacteria phage Felix 01]</Hit_def>
<Hit_accession>23609</Hit_accession>
<Hit_len>293</Hit_len>
<Hit_hsps>
<Hsp>
<Hsp_num>1</Hsp_num>
<Hsp_bit-score>49.2914</Hsp_bit-score>
<Hsp_score>116</Hsp_score>
<Hsp_evalue>5.15408e-06</Hsp_evalue>
<Hsp_query-from>580</Hsp_query-from>
<Hsp_query-to>792</Hsp_query-to>
<Hsp_hit-from>202</Hsp_hit-from>
<Hsp_hit-to>273</Hsp_hit-to>
<Hsp_query-frame>-1</Hsp_query-frame>
<Hsp_identity>26</Hsp_identity>
<Hsp_positive>45</Hsp_positive>
<Hsp_gaps>2</Hsp_gaps>
<Hsp_align-len>73</Hsp_align-len>
<Hsp_qseq>MHIIGDVE--GRTCILVDDMVDTAGTLCHAAKALKERGAAKVYAYCTHPVLSGRAIENIENSVLDELVVTNTI</Hsp_qseq>
<Hsp_hseq>MRILDDVDLTDKTVMILDDICDGGRTFVEAAKHLREAGAKRVELYVTHGIFS-KDVENLLDNGIDHIYTTNSL</Hsp_hseq>
<Hsp_midline>M I+ DV+ +T +++DD+ D T AAK L+E GA +V Y TH + S + +EN+ ++ +D + TN++</Hsp_midline>
</Hsp>
</Hit_hsps>
</Hit>
<Hit>
<Hit_num>2</Hit_num>
<Hit_id>gnl|BL_ORD_ID|2466</Hit_id>
<Hit_def>gi|51557505|ref|YP_068339.1| large tegument protein [Suid herpesvirus 1]</Hit_def>
<Hit_accession>2466</Hit_accession>
<Hit_len>3084</Hit_len>
<Hit_hsps>
<Hsp>
<Hsp_num>1</Hsp_num>
<Hsp_bit-score>48.9062</Hsp_bit-score>
<Hsp_score>115</Hsp_score>
<Hsp_evalue>6.70494e-06</Hsp_evalue>
<Hsp_query-from>369</Hsp_query-from>
<Hsp_query-to>875</Hsp_query-to>
<Hsp_hit-from>2312</Hsp_hit-from>
<Hsp_hit-to>2465</Hsp_hit-to>
<Hsp_query-frame>-2</Hsp_query-frame>
<Hsp_identity>52</Hsp_identity>
<Hsp_positive>70</Hsp_positive>
<Hsp_gaps>4</Hsp_gaps>
<Hsp_align-len>173</Hsp_align-len>
<Hsp_qseq>APESQEPGASTWRSSTSVVKKGQPSQK*CTSSVTSKAVPASWSTTWSTLPAPCATPPKR*KSAAPPRSTPTAPTRCCPAAPSRTSRIPSWTSWWSPTPSRCPLRRSPARVFASSTSPR-SSPKRSAASATKNRSAP---CSAKRNWPDHTAPPRAGLFALPPEAGRKPQGGLV</Hsp_qseq>
<Hsp_hseq>APPAQKPPAQPATAAATTAPKATPQTQPPTRAQTQTAPPPPSAAT-----AAAQVPPQ------PPSSQPAAKPRGAPPAPPAPP--PPSAQTTLPRPAAPPAPPPPS---AQTTLPRPAPPPPSAPAATPTPPAPGPAPSAKKSDGDRIVEPKAG---APPDVRDAKFGGKV</Hsp_hseq>
<Hsp_midline>AP +Q+P A ++ + K P + T + T A P + T A PP+ PP S P A R P AP P P P+ P P+ A +T PR + P SA +AT AP SAK++ D P+AG PP+ GG V</Hsp_midline>
</Hsp>
</Hit_hsps>
</Hit>
</Iteration_hits>
<Iteration_stat>
<Statistics>
<Statistics_db-num>68007</Statistics_db-num>
<Statistics_db-len>19518578</Statistics_db-len>
<Statistics_hsp-len>0</Statistics_hsp-len>
<Statistics_eff-space>0</Statistics_eff-space>
<Statistics_kappa>0.041</Statistics_kappa>
<Statistics_lambda>0.267</Statistics_lambda>
<Statistics_entropy>0.14</Statistics_entropy>
</Statistics>
</Iteration_stat>
</Iteration>
기본적으로 각 쿼리 검색은 반복 요소를 만듭니다. 각 반복에는 여러 번의 히트가있어 여러 HSP를 가질 수 있습니다. 첫 번째 히트 만 받고 싶습니다. 각 반복에서 첫 번째 HSP입니다. 폭발로 인해 맞지 않으면 반복을 무시하고 싶습니다. 이 간단한 코드를 작성했습니다.
#!/usr/bin/env python
from elementtree.ElementTree import parse
from elementtree import ElementTree as ET
file = open("/Applications/blast/blanes_viral_nr_results.xml", "r")
save_file = open("/Applications/blast/Blast_parse_ET.txt", 'w')
tree = parse(file)
elem = tree.getroot()
print elem
Per_ID = ()
save_file.write('>%s\t%s\t%s\t%s\t%s\t%s\t\n\n\n\n' % ("It_Num\t", "It_ID\t", "Hit_Def\t", "Num\t", "ID\t", "ACC\t"))
iteration = tree.findall('BlastOutput_iterations/Iteration')
for iteration in iteration:
for hit in iteration.findall('Iteration_hits/Hit'):
It_Num = iteration.findtext('Iteration_iter-num')
It_ID = iteration.findtext('Iteration_query-def')
Hit_Def = hit.findtext('Hit_def')
Num = hit.findtext('Hit_num')
ID = hit.findtext('Hit_id')
DEF = hit.findtext('Hit_def')
ACC = hit.findtext('Hit_accession')
save_file.write('>%s\t%s\t%s\t%s\t%s\t%s\t' % (It_Num, It_ID[12:26], Hit_Def[1:10], Num, ID, ACC,))
for hsp in hit.findall('Hit_hsps'):
HSPN = hsp.findtext('Hsp/Hsp_num')
identities = hsp.findtext('Hsp/Hsp_identity')
#print 'id: ', identities.rjust(4),
length = hsp.findtext('Hsp/Hsp_align-len')
#print 'len:', length.rjust(4),
Per_ID = int(identities) * 100.0 / int(length)
#print hsp.findtext('Hsp/Hsp_qseq')[:50]
#print hsp.findtext('Hsp/Hsp_midline')[:50]
#print hsp.findtext('Hsp/Hsp_hseq')[:50]
save_file.write('%s\t%s\t%s\%st\n' % ('***', '%', HSPN, Per_ID))
save_file.write('n\n' % ())
모든 도움이 크게 걱정 될 것입니다!
해결책
자신의 파서를 구축하는 것은 "재미있는"일 수 있지만 이미 XML 파일을 파일을 구문 분석 할 수있는 패키지가 있습니다. 원하는 경우 로컬 블래스트 인스턴스의 중간 호출을 수행 할 수도 있습니다.
기본 사이트는 여기에 있습니다.http://biopython.org/wiki/biopython
XML Blast Parser는 여기에 있습니다.http://biopython.org/dist/docs/tutorial/tutorial.html#htoc82
같은 것 :
from Bio.Blast import NCBIXML
with open('xml/results/file') as handle:
all_records = NCBIXML.parse(handle)
first_record = all_records.next()
작동해야합니다. 나는 일반적으로 Biopython Parsers와 Writers를 좋아하지만 계급 구조 조직을 좋아하지 않습니다. 그래서 나는 일반적으로 파서를 사용하고 필요한 정보를 내 구조로 추출합니다. ymmv
도움이되기를 바랍니다.
다른 팁
두 번째로 Judowill이 제안한 것 - Biopython Parser에 대해서는 더 똑똑하지 않습니다. 이것은 당신에게 조금 더 나아가 야합니다 :
from Bio.Blast import NCBIXML
blast = NCBIXML.parse(open('results.xml','rU'))
for record in blast:
if record.alignments:
# to print the "best" matches e-score
print record.alignments[0].hsps[0].expect
# to print the "best" matches bit-score
print record.alignments[0].hsps[0].score
break
첫 번째 쿼리 후에는 중지됩니다 (첫 번째 및 최상의 일치를 반환). 같은 파일 내의 다른 쿼리에 대한 결과를 원할 것 같아요. break 마지막 줄에서.
기본 요구 사항을 이해하면 쿼리 단백질/뉴클레오티드 시퀀스의 최상위 적중/HSP를 얻으려면 NR/NT 데이터베이스를 사용하여 시스템에 독립형 폭발을 설치하지 않습니다. 옵션을 입력하십시오
blastall -p {blast programme blastp for protein,blastn for nucleotide} -d {database} -i {input query} -v 1{for top hit} -b 1{alignment of the top hit with query} -m 7{xml blast output} -o example.xml
MS Excel에서 XML 출력 파일을 엽니 다.
제휴하지 않습니다 StackOverflow
