Hello everyone, I'm new to variant calling and I'm having trouble running VarScan2 v2.4.4 (Support Protocol 1) which can be found in the following online document: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4278659 /
The authors advise using the fpfilter.pl accessory script (https://sourceforge.net/projects/varscan/files/scripts/). When I run the following command:
perl script_FPfilter sample.varScan.snp.filter.bed varScan.variants.readcounts –output-basename varScan.variants.filter
I get the following output:
27841 variants
27841 failed to get readcounts for variant allele
0 had read position < 0.1
0 had strandedness < 0.01
0 had var_count < 4
0 had var_freq < 0.05
0 had mismatch qualsum difference > 50
0 had mapping quality difference > 30
0 had read length difference > 25
0 had var_distance_to_3' < 0.2
0 passed the strand filter
Apparently, the problem is with obtaining the readcounts.
The bam-readcounts files were obtained as follows:
create indexed bam files
samtools index bqsr_output.bam
create bed file for sample.varScan.snp.filter
vcf2bed --snvs sample.varScan.snp.filter > sample.varScan.snp.filter.bed
run bam-read-count
bam-readcount -q 1 -b 20 -f hg38_genome.fa -l sample.varScan.snp.filter.bed bqsr_output.bam > snp_varScan.variants.readcounts
what am I doing wrong?
Can you show a
head
of the BED and bam-rc output?output Bed
output bam-readcount