How can I calculate query coverage using perl If I have a file like this. I want to calculate it on my own using perl.
Query id, q. start, q. end, gene_id Q.length Calculated_q_coverage
lcl|NZ_CP018664.1_gene_1 1 1149 NZ_CP018664.1_gene_1 1149
lcl|NZ_CP018664.1_gene_2 1 837 NZ_CP018664.1_gene_2 837
lcl|NZ_CP018664.1_gene_3 1 606 NZ_CP018664.1_gene_3 606
lcl|NZ_CP018664.1_gene_4 1 327 NZ_CP018664.1_gene_4 327
lcl|NZ_CP018664.1_gene_5 1 471 NZ_CP018664.1_gene_5 471
lcl|NZ_CP018664.1_gene_6 1 363 NZ_CP018664.1_gene_6 363
lcl|NZ_CP018664.1_gene_7 1 645 NZ_CP018664.1_gene_7 645
•
written
by
sharmatina189059 •
1 answer
Surprising answer maybe, you cannot calculate the correct query coverage from tabular output of blast.
Query coverage is not equal to: abs(qend-qstart)/qlength. While often weekly defined, I normally have seen query coverage as the proportion of sequence covered by the alignment. Why is this != to the difference of start and end? Because of the existence of HSPs and gaps. It could well be the case that there are multiple HSPs or large gaps in the alignment, such that parts of the sequence between qstart and qend are not covered by the alignment and therefore do not count for the coverage. If I am not mistaken, blast tabular report contains 1 line per hit, not per HSP.
Think of it in a similar way as with mRNA length vs gene length: the mRNA length is not equal to the difference of genomic start and stop of the transcript but the combined length of all exons in the transcript.
what do you mean by relative coverage?
Log in to answer this question.
I have these 4 fields. Query id, q. start, q. end, gene_id and want to calculate query coverage using perl.
Please use
ADD COMMENT/ADD REPLYwhen responding to existing posts to keep threads logically organized. Additional information should be added to the original post by editing it.Also for reference, you should validate (up-vote or accept answers) comments/posts that help with your questions. That shows appreciation for those who help and also provides closure to open threads.
I added markup to your post for increased readability. You can do this by selecting the text and clicking the 101010 button. When you compose or edit a post that button is in your toolbar, see image below:
want to calculate query coverage using perl <=== what that mean ? you already have coverage field 5 or you want the same result as field 5 ?
How is coverage defined in your case? Is it q.end-q.start? If so, you can do
extending a little more this;
Why did not you divide it by query length. I could not get the hang of it. Please elaborate.
Because in your example your coverage field = query length ....