2
$\begingroup$

I need to extract (or to calculate precisely) the query coverage from Blast output. I am using NCBIWWW.qblast to query Blast in XML format.

Is it possible to know the query coverage parsing the XML output? If yes, how?

result_handle_xml = NCBIWWW.qblast("blastp","pdb",my_query.seq,format_type='XML')

here it suggest to use '-6 qcovs', but only for BLAST command line version and only for tabular format

This is the fasta file that I open:

>sp|P02489|CRYAA_HUMAN Alpha-crystallin A chain OS=Homo sapiens 
GN=CRYAA PE=1 SV=2
MDVTIQHPWFKRTLGPFYPSRLFDQFFGEGLFEYDLLPFLSSTISPYYRQSLFRTVLDSG
ISEVRSDRDKFVIFLDVKHFSPEDLTVKVQDDFVEIHGKHNERQDDHGYISREFHRRYRL
PSNVDQSALSCSLSADGMLTFCGPKIQTGLDATHAERAIPVSREEKPTSAPSS
$\endgroup$
3
  • $\begingroup$ Please edit your question and include an example of the file you are trying to parse. We can't really help you parse data you don't show. $\endgroup$
    – terdon
    Dec 28, 2017 at 10:30
  • $\begingroup$ I see you edited, but you still haven't shown us the file you are trying to parse. Is the information even there? Please add an example so we know what you're working with. $\endgroup$
    – terdon
    Dec 28, 2017 at 16:55
  • 2
    $\begingroup$ According to this, the answer is 'no'. $\endgroup$ Dec 28, 2017 at 21:19

2 Answers 2

0
$\begingroup$

I had the same problem. I'm pretty sure you can calculate it yourself using the following:

  1. The length of your sequence (seq_len)

  2. The value from the <Hsp_query-from> tag from the hit (<Hsp_query-from>)

  3. The value from the <Hsp_query-to> tag from the hit (<Hsp_query-to>)

It can be calculated by:

query_cover = (<Hsp_query-to> - <Hsp_query-from> + 1) / seq_len * 100

Please note that I did not get this from any documentation. I figured it out by staring at the xml output returned from the Biopython NCBIWWW.qblast() function. This solution seems to match with the results I get by blasting the same sequence on the interactive web portal.

$\endgroup$
0
$\begingroup$
record = NCBIXML.read(result_handle_xml)
query_len = record.query_length
for alignment in record.alignments:
    _ = alignment.hsps.__next__()
    query_cover = (_.query_end - _.query_start + 1) / query_len * 100
$\endgroup$
2
  • $\begingroup$ Could you please check your code? The use of '_' I always used/saw being used as a null variable not a temporary variable as done here. Please consider your methods against against -outfmt 6 std qcovs which is being suggested as the preferred option. $\endgroup$
    – M__
    Apr 7 at 16:52
  • $\begingroup$ Okay '_' is used in loops. It is a rare usage IMO $\endgroup$
    – M__
    Apr 7 at 17:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.