Problem with qcovs in blast outfmt 6

Question

[Bara'a]

Hi all ...

I'm using Biopython 3.4.1 with Blast 2.2.29 and I'm trying to get the query coverage of some blast search beside the default report info in tabular format using -outfmt 6 as following :

comline=NcbiblastxCommandline(cmd='blastn',query="A.fasta",db="B.db",evalue=0.001,outfmt='6 qcovs',out="result.txt")

But it throws the following error:

Error: Too many positional arguments (1), the offending value: qcovs

I know it seems a very basic question ,but it's seriously driving me crazy ,since I spent two days trying to figure out the solution for this error and nothing seems to fix it !!

Some threads have addressed this issue and suggested to rewrite the command by hand instead of copying and pasting ... I did rewrite it by hand like thousand times and the error pertains .

Also, I have tried to run it from the Command Prompt cmd and again it threw the same error!!

Is there anyway to get around this problem?

Many thanks in advance

[Bara'a]

It's really funny how deep searches sometimes leads you to nowhere , while only chance can help you to find the solution !!

I found this by accident Does Ncbitblastncommandline Reads The Outfmt (Options 6, 7, And 10 ) Correctly? and I'm surprisingly mad for wasting 2 days of my life chasing a double quote criminal !!

I hate when the error is that tiny to waste that much time ... it's a terrible feeling.

The only thing that made me laugh is that the one who solved this obstacle is the same one who always answers my questions in this forum ... it's Peter from Scotland ,UK :D

He is truly a life saver ... Is there away I can thank him properly beside Upvoting and accepting his answers?!

[Peter]

I guess just pass it forward - by which I mean help other people when you can. :)

[nikulina]

Thank you Peter, you saved me from insanity!

Answer 1

[Michael]

Did you try outfmt='"6 qcovs"'? The command string needs to literally contain these quotes, such that the command line looks like:

blastn -outfmt "6 qcovs" ...

[Michael]

Well, I'm not Peter from Scotland, but I submitted in the same instant seemingly as you submitted your comment, thanks for the upvotes anyway :D

[Bara'a]

Yes ... I found a thread 10 minutes ago with your same solution.

Thanks anyways ... I do appreciate your help .

[HG]

I was having a try as you suggested. But my result coming like this:

Query id                                  Database id                               % Identity     Alignment length     Mismatches     Gap     QSS*       QSE**      DSS***     DSE****     e-value     Bit score
NODE_1_length_612857_cov_36.6605_ID_1     (Flq)OqxBgb:EU370913:47851-51003:3153     28.88          1049                 695            16      558271     561309     4          1037        1e-90       321
NODE_1_length_612857_cov_36.6605_ID_1     (Flq)OqxBgb:EU370913:47851-51003:3153     30             1030                 665            17      561364     564330     3          1017        1e-78       284
NODE_1_length_612857_cov_36.6605_ID_1     (Flq)OqxA:EU370913:46652-47827:1176       28.87          336                  230            3       557153     558133     35         370         8e-31       126
NODE_1_length_612857_cov_36.6605_ID_1     (Flq)OqxA:EU370913:46652-47827:1176       31.3           246                  155            6       103747     104445     59         303         5e-11       64.3
NODE_1_length_612857_cov_36.6605_ID_1     (Flq)QepA2:EU847537:1672-3107:1536        29.59          392                  264            5       291389     292555     35         417         7e-18       87
NODE_1_length_612857_cov_36.6605_ID_1     (Flq)QepA2:EU847537:1672-3107:1536        26.01          419                  265            15      564421     565608     4          400         2e-04       43.9
NODE_1_length_612857_cov_36.6605_ID_1     (Flq)QepA2:EU847537:1672-3107:1536        29.55          132                  92             1       52523      52915      44         175         0.028       36.6
NODE_1_length_612857_cov_36.6605_ID_1     (Flq)QepA2:EU847537:1672-3107:1536        34.62          104                  65             2       122466     122771     77         179         0.029       36.6

Could you please suggest me did I missed something. I cant find query coverage in result.

[Bara'a]

@HG Could you please write the code snippet you used to have those results?!

We can't guess whether the fault is in your code or something else!

Answer 2

[Bara'a]

Hey @HG :)

The parameter qcovs is not of the defaults!

You should write it explicitly in your code ... something like:

comline=NcbiblastxCommandline(cmd='blastn', query="Test.fasta", db="Test2.db", evalue=0.001, outfmt='"6 std qlen slen qseq sseq qcovs"', out="Result.txt" )

Note that std represents the default parameters reported to the out file (Those you've got in your trial) whereas qlen , slen , qseq , sseq and qcovs are NOT!

Good luck :)

[HG]

Hi

I used like in such away:

makeblastdb -in Db_Flq.fasta -dbtype prot;
for R1 in *.fasta
do blastx -db Db_Flq.fasta -query $R1 -out $R1.blasted_Flq -outfmt "6 qcovs"...;
done

Is it proper way to do or else I am missing something?

[Bara'a]

@HG What programming language are you using? Is that pseudo code?

Anyways ... You're missing the single quote wrapper!

outfmt= '"6 std qlen slen qseq sseq qcovs"'

[HG]

Thanks for reply. I am going to have a try

[HG]

It's working. Could you please let me know how can we print specific field and comment few field: For example I don't want to print last field and 396559 in final result? Is there any command line option? or else I have to do by other process.

ODE_1_length_396559_cov_64.0552_ID_1    (Flq)QepA:AB263754:7052-8587:1536    24.91    285    184    10    293103    293879    76    356    5e-08    54.7    396559    511    LGPLSDRIGRRPVLI

[Bara'a]

@HG ... Hello once again :)

You can always have the fields you want in any order just by writing them explicitly without using std .

It will be something like this :

'"6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore qlen slen qseq sseq qcovs"'

Note that the first 12 parameters here represent what std usually reports in the exact same order , but following this approach ( write fields explicitly ) allows you to print them in any order you wish to have .

To add or eliminate other fields , refer to their names in your code according to this: Blast - Formatting Output

Please, don't forget to upvote my answer if you find it correct ;)

[HG]

Sorry to disturb you again: my out put is coming like this

My bash command:

blastx -db DbFlq.fasta -query H1.fasta -out H1.fasta.blasted_test1 -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore qlen slen qcovs"

Output: I did not understand what is the last filed 0 is here is it "qcovs" bit confusing, could you please help me out

NODE_1_length_255599_cov_23.9893_ID_1     (Flq)norA:D90119:478-1644:1167            30.88     68       46      1      245612     245409     2       68       0.001     39.3     255599     388      0
NODE_1_length_255599_cov_23.9893_ID_1     (Flq)norA:D90119:478-1644:1167            24.16     149      105     2      102981     102547     5       149      0.02      35.4     255599     388      0
NODE_1_length_255599_cov_23.9893_ID_1     (Flq)norA:D90119:478-1644:1167            27.27     77       44      3      166288     166082     5       77       9.8       26.9     255599     388      0
NODE_1_length_255599_cov_23.9893_ID_1     (Flq)OqxA:EU370913:46652-47827:1176       24.18     91       63      1      251849     251595     175     265      2.4       28.9     255599     391      0
NODE_1_length_255599_cov_23.9893_ID_1     (Flq)QnrC:EU917444:1717-2382:666          36.84     38       21      1      222270     222374     131     168      9         26.6     255599     221      0
NODE_2_length_229923_cov_24.0916_ID_3     (Flq)QnrC:EU917444:1717-2382:666          37.21     43       27      0      227934     227806     87      129      7.3       26.9     229923     221      0
NODE_2_length_229923_cov_24.0916_ID_3     (Flq)QnrB46:JQ349154:1-680:680            24.73     93       55      3      57220      57462      128     217      9.1       26.6     229923     226      0
NODE_3_length_192885_cov_24.5171_ID_5     (Flq)OqxBgb:EU370913:47851-51003:3153     40.29     1050     600     11     5778       2677       4       1042     0         648      192885     1050     2
NODE_3_length_192885_cov_24.5171_ID_5     (Flq)OqxBgb:EU370913:47851-51003:3153     27.63     76       54      1      16032      16259      740     814      2.1       28.9     192885     1050     2

[Bara'a]

@HG ... Although your sample is disturbing ant not comfortable to look at , but it follows the exact order of parameter's appearance in your bash command, so yes: 0 is qcovs.

If you find it confusing you can just eliminate it from the command and it will no longer appears in the output.

You're welcome anytime :)

Answer 3

[Sam]

Hy everyone! I'm trying BLASTn, with the command line:

"blastn -db ensembl34.fa -query chr0.fna -word_size 7 -dust yes -num_threads 5 -evalue 0.0001 -outfmt "6 std qcovs qcovhsp" > res_teste0.csv"

But the "qcovs" and the "qcovshsp" parameters are both "0" in the results... Anyone knows why?

[GenoMax]

I don't think you are able to redirect output like you are doing. You should write the result file out with a -out filename option. What did you get in res_teste0.csv file out of curiosity?

[Sam]

Hi, I got all the results correctly, except for the "qcovs"

[GenoMax]

Interesting.

I think the problem may be the nested double quotes that are in your command because of the output redirection you are using. Can you try following?

blastn -db ensembl34.fa -query chr0.fna -word_size 7 -dust yes -num_threads 5 -evalue 0.0001 -out res_teste0.csv -outfmt "6 std qcovs qcovhsp"

or if you want to keep the redirection then try this (ref: Does Ncbitblastncommandline Reads The Outfmt (Options 6, 7, And 10 ) Correctly? )

"blastn -db ensembl34.fa -query chr0.fna -word_size 7 -dust yes -num_threads 5 -evalue 0.0001 -outfmt ' "6 std qcovs qcovhsp" ' > res_teste0.csv"

[Sam]

Ok, I'll try. Super thx! :)