Virfinder AWK Filtering Maybe not Correct

[benyoung93]

Hi there :)

First of all, a wonderful collection of tools and a pipeline !! I have had trouble with some if the installation but I am just chopping and cutting the pipeline to fit my needs.

A quick query for the Virfinder step. So I know you want to do q-value < 0.01 and length > 1000bp. I have run the virfinder rscript file and got the output folder. I just want to double check that the filtering awk script is correct.

head P10_virfinder.tsv

"name"  "length"    "score" "pvalue"    "qvalue"
"1" "k141_270856 flag=1 multi=2.0000 len=320"   320 0.138513157806992   0.588997020854022   0.56810885275558
"2" "k141_33860 flag=1 multi=2.0000 len=361"    361 0.248927602838187   0.42085402184707    0.525288351341074
"3" "k141_237000 flag=1 multi=6.0000 len=478"   478 0.392398453460999   0.271737835153923   0.470800520085438
"4" "k141_93112 flag=1 multi=5.0000 len=419"    419 0.316555183433318   0.343932472691162   0.499862326324076
"5" "k141_177752 flag=1 multi=2.0000 len=398"   398 0.461420194784674   0.219920556107249   0.450421369675625
"6" "k141_321640 flag=1 multi=4.0000 len=383"   383 0.456895550455814   0.222641509433962   0.450905940743768
"7" "k141_211608 flag=1 multi=2.0000 len=477"   477 0.152712635250371   0.562542204568024   0.561554065996342
"8" "k141_16930 flag=1 multi=2.0000 len=426"    426 0.165548747858572   0.539602780536246   0.555654697140864
"9" "k141_8465 flag=1 multi=3.0000 len=431" 431 0.855015816724245   0.0269116186693148  0.268291393791584

The subsequent awk script is then this obviously (ignore all mu '"'"'" this is for a loop script that generates jobs for all of my samples).

cat '"$PALPAL"'_virfinder.tsv | awk -F'\t' '"'"'{ if ($4 <= 0.01) print }'"'"' | awk -F'_' '"'"'{ if ($4 >= 1000) print }'"'"' | cut -f2 | sed "s/\"//g" > '"$PALPAL"'_virfinder_filtered_data.txt

So this awk script is currently firstly using awk on the fourth \t column which is p-value ($4), should this not be $5 as the <0.01 filtering is for q-value ?

Secondly, second awk is extracting the fourth instance after underscore, looking at the output this doesn't seem to be correct. If you are wanting length >1000 would it not be a better idea to do awk -F'\t' '"'"'{ if ($2 >= 1000) print }'"'"'

Looking at the virfinder github it seems the original results format would work with your awk script, but I think a new update may of switched this ?

I used the mamba install in the mudoger install scripts but edited a wee bit (mamba create -n virfinder_env -c bioconda r-virfinder) so I believe both my local and the mudogoer install versions would be the same.

Happy to provide more info if needed, and I apologise in advance if I am barking up the wrong tree, but I was doing some QC and testing and noticed that the awk was not approproate for the output file generated by virfinder :).

Ben

[benyoung93]

Also, additional thing I have just realised, you need to identify your column and + 1 as the awk takes into account the rowname column in the output from write.table. It may be prudent therefore to put a row.names=F in the r script, or for the awk it would need to be $6 =<0.01 and $3 =>100 for q-value and length respectively