No variant found using GATK 4.0 HaplotypeCaller

Question

I am doing variant calling on RNA-seq datasets from wheat which is hexaploid,the binary alignment (BAM) files were created using STAR version 2.6.0c and variant calling was done using GATK 4.0 HaplotypeCaller.The whole pipeline is as follows:

Alignment through STAR 2.6.0c :

STAR --runMode alignReads --runThreadN 20  --outSAMtype BAM  SortedByCoordinate --sjdbOverhang 100  --genomeDir wheat --readFilesIn G1_cleaned_R1.fastq G1_cleaned_R2.fastq --sjdbGTFtagExonParentTranscript wheat.gff3

Then AddOrReplaceReadGroups from picard was used with default settings:

 AddOrReplaceReadGroups \
   I=input.bam \
   O=Read_groups_added.bam \
   RGID=4 \
   RGLB=lib1 \
   RGPL=illumina \
   RGPU=unit1 \
   RGSM=20

Finally GATK 4.0 haplotypeCaller was applied:

gatk HaplotypeCaller --native-pair-hmm-threads 28  \
    -I Read_groups_added.bam -R wheat_ref.fa -O final.vcf

The command ran successfully as shown below:

20:33:01.018 INFO  ProgressMeter - TGACv1_scaffold_470693_5DS:301            102.2              45108670         441320.0
20:33:11.039 INFO  ProgressMeter - TGACv1_scaffold_724195_U:301            102.4              45111270         440625.5
20:33:21.072 INFO  ProgressMeter - TGACv1_scaffold_724579_U:301            102.5              45113880         439932.4
20:33:31.108 INFO  ProgressMeter - TGACv1_scaffold_248869_3B:301            102.7              45116440         439240.9
20:33:41.143 INFO  ProgressMeter - TGACv1_scaffold_248918_3B:301            102.9              45118990         438551.6
20:33:51.150 INFO  ProgressMeter - TGACv1_scaffold_306137_4AL:301            103.0              45121550         437866.7
20:34:01.176 INFO  ProgressMeter - TGACv1_scaffold_193559_2DS:301            103.2              45124100         437182.5
20:34:03.134 INFO  HaplotypeCaller - 100112154 read(s) filtered by: ((((((((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter) AND NotDuplicateReadFilter) AND PassesVendorQualityCheckReadFilter) AND NonZeroReferenceLengthAlignmentReadFilter) AND GoodCigarReadFilter) AND WellformedReadFilter)
100112154 read(s) filtered by: (((((((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter) AND NotDuplicateReadFilter) AND PassesVendorQualityCheckReadFilter) AND NonZeroReferenceLengthAlignmentReadFilter) AND GoodCigarReadFilter)
100112154 read(s) filtered by: ((((((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter) AND NotDuplicateReadFilter) AND PassesVendorQualityCheckReadFilter) AND NonZeroReferenceLengthAlignmentReadFilter)
100112154 read(s) filtered by: (((((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter) AND NotDuplicateReadFilter) AND PassesVendorQualityCheckReadFilter)
100112154 read(s) filtered by: ((((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter) AND NotDuplicateReadFilter)
100112154 read(s) filtered by: (((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter)
100112154 read(s) filtered by: ((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter)
100112154 read(s) filtered by: (MappingQualityReadFilter AND MappingQualityAvailableReadFilter)
21066679 read(s) filtered by: MappingQualityReadFilter 
79045475 read(s) filtered by: MappingQualityAvailableReadFilter 

20:34:03.134 INFO  ProgressMeter - TGACv1_scaffold_725922_U:301            103.2              45124592         437049.1
20:34:03.134 INFO  ProgressMeter - Traversal complete. Processed 45124592 total regions in 103.2 minutes.
20:34:04.749 INFO  VectorLoglessPairHMM - Time spent in setup for JNI call : 0.0
20:34:04.749 INFO  PairHMM - Total compute time in PairHMM computeLogLikelihoods() : 0.0
20:34:04.749 INFO  SmithWatermanAligner - Total compute time in java Smith-Waterman : 0.00 sec
20:34:04.749 INFO  HaplotypeCaller - Shutting down engine
[July 3, 2018 8:34:04 PM PKT]   org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller done. Elapsed time: 103.99 minutes.
Runtime.totalMemory()=4789370880

but still the generated vcf file i.e final.vcf does not contain any variants as shown below:

[dellemc@localhost bilal]$ tail final.vcf
##contig=<ID=TGACv1_scaffold_725916_U,length=500>
##contig=<ID=TGACv1_scaffold_725917_U,length=500>
##contig=<ID=TGACv1_scaffold_725918_U,length=500>
##contig=<ID=TGACv1_scaffold_725919_U,length=500>
##contig=<ID=TGACv1_scaffold_725920_U,length=500>
##contig=<ID=TGACv1_scaffold_725921_U,length=500>
##contig=<ID=TGACv1_scaffold_725922_U,length=500>
##contig=<ID=TGACv1_scaffold_731344_7BL,length=500>
##source=HaplotypeCaller
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  20

Is there any mistake in the haplotypeCaller command or any argument is missing ?

Answer 1

There are steps of the 'best practices' workflow missing, which I think is leading to the majority of reads being filtered by the haplotype caller, notice these lines in the output:

21066679 read(s) filtered by: MappingQualityReadFilter 
79045475 read(s) filtered by: MappingQualityAvailableReadFilter 

That's >100m reads being filtered by the HC.

Specifically, you're missing SplitNCigarReads, which separates reads that bridge a splice junction (padded by Ns in the STAR alignment) and reassigns the mapping quality. At the moment this element of the workflow is not validated or available in GATK4 (RNA-Seq variant calling is largely abandonware in newer versions of GATK, AFAIK), but in the best practices WDL utilises GATK3.5 - see https://github.com/gatk-workflows/gatk3-4-rnaseq-germline-snps-indels/blob/master/rna-germline-variant-calling.wdl, line 440-477.

There are other elements of the recommended workflow missing too, but I think this is the critical bit. I won't harp on here about the unsuitability of RNA-Seq data for variant calling, as I'm all too aware there often isn't much option.

Answer 2

As well as the issues @sjcockell mentions, you have another potential problem you've got a hexaploid organism, by default the HaplotypeCaller (HC) is hardwired to a ploidy of 2, you need to set this to 6 to match your organisms ploidy using -ploidy 6 as an argument. Since the ploidy is a key term in the underlying Bayesian stats the HC employs (see this link for an explanation) it's likely that some of the lower allelic frequencies that hexaploidy would generate will mean quite a few variants won't be called as they'd be statistically unlikely for a diploid organism.