Question: Hi,

I am attempting to reanalyze some old data. Recently, an updated genome was released for the organism I work with. I had to map the reads to the new genome and what not. So what I have now is a bam file but now I want to do some variant calling so I can generate a vcf.

I know there is a lot of programs out there to do this but I am attempting to find the best one. I have tried using Platypus but the vcf file that it generates never contains any genome data. That is probably my fault but I can't seem to figure it out and there isn't much support. I tried VarDict but when I go to look at the vcf file, I get an error of "failure to parse tbx_vcf". I can really seem to find why I am getting that error. I have tried Freebayes and that has worked for me, but I have an issue with indels. Freebayes seems to call SNPs instead of acknowledging indels. I know there is GATK, but it seem like such a tedious piece of software. It looks like you have to go through 15 steps to just get anything worthwhile. I could be looking at the wrong examples. What I have notice about GATK is they try to make it so simple that it ends up just making things complicated.

Could anyone maybe give some suggestions or maybe some insight? Any help is very appreciated

Thanks!!

Answer: **Edit September 10, 2019:**
----
Although I have used samtools / bcftools mpileup, I have doubts surrounding its false positive rate; although, this can be controlled by setting a minimum *total* position read depth for each variant call (e.g. 30 reads). mpileup is neither fine-tuned to identify indels. Were I to implement a clinical pipeline in the future for calling germline variants, I would likely run both GATK and the method described below with mpileup, and then cross compare.

-----------------
------------

**Original answer**
----
While the GATK is popular, from my experience from analysing both clinical and research NGS data, GATK is *not* the most sensitive for detecting germline single nucleotide variants (SNVs). Sensitivity here is gauged by comparing NGS results to those of Sanger sequencing.

As you've mentioned, there are many variant callers out there and each has advantages, including the GATK. A fool would be the person who would come here and say that this or that particular program / pipeline is the best, as they are each good in different scenarios.

My 'go to' pipeline is the 'old faithful':

1. *BWA mem* for alignment (if read lengths > 70bp)
2. Remove low MAPQ reads
3. Mark and remove PCR / optical duplicates
4. Produce random BAM subsets per sample (usually 3 subsets at 25%,
50%, and 75% 'randomly' selected reads)
5. Call variants with *BCFtools mpileup* (main BAM and 3 subsets
together) piped into *BCFtools call*

The final command would be something like:

bcftools mpileup \
--redo-BAQ \
--min-BQ 30 \
--per-sample-mF \
--output-tags DP,AD \
-f "${Ref_FASTA}" \
--BCF Full.bam 75pcReads.bam 50pcReads.bam 25pcReads.bam | \
bcftools call \
--multiallelic-caller \
--variants-only \
-Ob > Variants.bcf

This has greater sensitivity for ***SNVs*** when compared to other pipelines that I've tested, which includes GATK and Illumina's BaseSpace pipeline. Granted, there are many programs and pipelines out there, and one would ideally have to test all of them, and on an expanded number of samples.

For somatic variant calling, I use Lancet.

Kevin