Come estrarre i primi elementi di successo da un file XML NCBI BLAST?
https://stackoverflow.com/questions/1935385
-
20-09-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianDomanda
Sto cercando di estrarre solo il primo hit da un file BLAST xml NCBI.successivamente vorrei ottenere solo il primo HSP.nella fase finale vorrei ottenerli in base al miglior punteggio.per chiarire le cose ecco un esempio del file 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>
fondamentalmente ogni ricerca di query crea un elemento Iterazione.ogni iterazione può avere più colpi che a loro volta possono avere più HSP.Vorrei ottenere solo il primo successo ed è il primo HSP da ogni iterazione.se BLAST non trovasse risultati, vorrei ignorare l'iterazione.Ho elaborato questo semplice codice:
#!/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' % ())
Qualsiasi aiuto sarebbe molto apprezzato!
Soluzione
Mentre costruire il proprio parser può essere "divertente" C'è già un pacchetto là fuori che può analizzare i file XML BLAST ... può anche fare la chiamata intermedia di un'istanza BLAST locale, per voi, se lo desiderano.
Il sito principale è qui: http://biopython.org/wiki/Biopython
e il BLAST parser XML è qui: http://biopython.org/DIST/docs/tutorial/Tutorial.html# htoc82
Qualcosa di simile:
from Bio.Blast import NCBIXML
with open('xml/results/file') as handle:
all_records = NCBIXML.parse(handle)
first_record = all_records.next()
Dovrebbe funzionare. Io in genere amo i parser Biopython e scrittori, ma non mi piace l'organizzazione di classe-struttura. Così ho in genere basta usare il parser ed estrarre le informazioni che ho bisogno nella mia propria struttura. YMMV
La speranza che aiuta.
Altri suggerimenti
Ti secondo quanto suggerito JudoWill - lavoro più intelligente, non di più con il parser Biopython. Questo dovrebbe farti un po 'oltre:
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
Questo fermerà dopo la prima query (ritorno il primo e miglior match). Da quando ho il sospetto che si potrebbe voler risultati per altre query all'interno dello stesso file, è sufficiente rimuovere la break dall'ultima riga.
Se ho capito i tuoi requisiti di base, allora vuoi ottenere il miglior risultato/HSP della sequenza proteina/nucleotide di query. Perché non installi l'esplosione autonoma sul tuo sistema con il database nr/nt formattato.Digita le opzioni
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
apri il file di output xml in MS Excel dove puoi vedere una forma tabulata del blastoutput con il singolo hit principale per ciascuna sequenza di query
