This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Problem in using VarScan fpfilter

Hi, I ran bam-readcount to get coverage per nucleotide. I got the result after getting this warning `

“WARNING: In read DHDC08P1_0264:2:2204:2494:140080#CTTGTA: Couldn't find single-end mapping quality. Check to see if the SM tag is in BAM”

This was the command:

bam-readcount -f ../index/hg38.fa -q 1 -b 20 -l all.Germline.hc.var normal.sorted.bam > all.Germline.hc.readcount

` and this is a line of the output file:

chr1    931131  C   14  =:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00  A:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00  C:14:37.00:34.86:7.00:2:12:0.38:0.03:9.93:2:0.19:123.86:0.72    G:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00  T:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00  N:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00  +CCCT:7:42.00:0.00:0.00:1:6:0.65:0.04:7.86:1:0.26:122.00:0.60

When I used the readcount file to filter out false positives from snp calling, I encountered this error:

Loading readcounts from all.Germline.hc.readcount...
Parsing variants from all.Germline.hc.var...
Exception thrown while filtering at chr1    931131: 3
java.lang.ArrayIndexOutOfBoundsException: 3
    at net.sf.varscan.FpFilter.<init>(FpFilter.java:365)
    at net.sf.varscan.VarScan.fpfilter(VarScan.java:352)
    at net.sf.varscan.VarScan.main(VarScan.java:174)

This was the command: java -jar VarScan.v2.3.9.jar fpfilter all.Germline.var all.Germline.readcount --min-ref-basequal 20 --min-var-basequal 20 --output-file all.Germline.hc.filtered.vcf --dream3-settings

Where is the problem?

Thanks very much for any help!

bam-readcount varscan fpfilter snp calling

all.Germline.hc.var must be a BED file. Also, you are not using the latest release of Varscan, which is I think 2.4.3. Better get the latest one from Github.

Thank you for your reply. As you suggested I got the latest release of VarScan( 2.4.3) but still get the same error. I think the problem is not related to my all.Germline.hc.var. As I found a BED file is a file with three columns (Chromosome Stop Start) and this is my all.Germline.hc.var file:

chr1    1046480 1046480
chr1    1182186 1182186
chr1    1217483 1217483

This is the command that I used to get the file:

awk 'BEGIN {OFS="\t"} {if (!/^#/) { isDel=(length($4) > length($5)) ? 1 : 0; print $1,($2+isDel),($2+isDel); }}' *.Germline.hc.vcf > all.Germline.hc.var

I am following this instruction link.

.

0 answers

No answers yet.

Log in to answer this question.