<p>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:</p>

<pre><code>21066679 read(s) filtered by: MappingQualityReadFilter 
79045475 read(s) filtered by: MappingQualityAvailableReadFilter 
</code></pre>

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

<p>Specifically, you're missing <code>SplitNCigarReads</code>, 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 <a href="https://github.com/gatk-workflows/gatk3-4-rnaseq-germline-snps-indels/blob/master/rna-germline-variant-calling.wdl" rel="nofollow noreferrer">https://github.com/gatk-workflows/gatk3-4-rnaseq-germline-snps-indels/blob/master/rna-germline-variant-calling.wdl</a>, line 440-477. </p>

<p>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.</p>

CC BY-SA 4.0
2018-07-05T10:18:38.837Z
What is meant by the option --limitOutSJcollapsed in the STAR command given in the link you suggested https://github.com/gatk-workflows/gatk3-4-rnaseq-germline-snps-indels/blob/a7a890b38af1309edc040ac89814df37a57f77b7/rna-germline-variant-calling.wdl#L345
Ammar Sabir Cheema
CC BY-SA 4.0
2018-07-05T10:56:35.930Z
A 'collapsed' junction is what gets written to the `SJ.out.tab` file during a run of STAR (the collapsing here refers to the removal of duplicate junctions). This option merely limits the maximum number of these junctions that are written. It defaults to 1,000,000.
sjcockell
CC BY-SA 4.0
2018-07-04T10:55:54.470Z
2018-07-04T10:55:54.470Z
sjcockell
<p>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 <code>-ploidy 6</code> as an argument.  Since the ploidy is a key term in the underlying Bayesian stats the HC employs (see this <a href="https://software.broadinstitute.org/gatk/documentation/article.php?id=4442" rel="nofollow noreferrer">link</a> 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. </p>

CC BY-SA 4.0
2018-07-16T10:28:24.580Z
You might also want to consider the GATK support forum too, as they are very helpful on GATK matters: https://gatkforums.broadinstitute.org/gatk/categories/gatk-support-forum
Matt Bashton
CC BY-SA 4.0
2018-07-13T15:04:11.323Z
2018-07-13T15:04:11.323Z
Matt Bashton
<p>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:</p>

<p>Alignment through <strong>STAR 2.6.0c</strong> :</p>

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

<p>Then <strong>AddOrReplaceReadGroups</strong> from <strong>picard</strong> was used with default settings:</p>

<pre><code> AddOrReplaceReadGroups \
   I=input.bam \
   O=Read_groups_added.bam \
   RGID=4 \
   RGLB=lib1 \
   RGPL=illumina \
   RGPU=unit1 \
   RGSM=20
</code></pre>

<p>Finally <strong>GATK 4.0 haplotypeCaller</strong> was applied:</p>

<pre><code>gatk HaplotypeCaller --native-pair-hmm-threads 28  \
    -I Read_groups_added.bam -R wheat_ref.fa -O final.vcf
</code></pre>

<p>The command ran successfully as shown below:</p>

<pre><code>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
</code></pre>

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

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

<p>Is there any mistake in the haplotypeCaller command or any argument is missing ?</p>

CC BY-SA 4.0
2018-07-03T23:21:44.710Z
All those compute times of 0.0 seem a little suspicious.
swbarnes2
CC BY-SA 4.0
2018-07-04T05:01:39.813Z
I have edited the question to add more details.
Ammar Sabir Cheema
CC BY-SA 4.0
2018-07-04T10:31:01.173Z
@swbarnes could you tell about the type of suspicions?
Ammar Sabir Cheema
<vcf><variant-calling><gatk>
No variant found using GATK 4.0 HaplotypeCaller
# No variant found using GATK 4.0 HaplotypeCaller

(Asked by: username_0 on 2018-07-03)

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=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##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 ?

**Tags:** vcf, variant-calling, gatk
> **Comments:**
> **username_1**: All those compute times of 0.0 seem a little suspicious.  
> **username_0**: I have edited the question to add more details.  
> **username_0**: @swbarnes could you tell about the type of suspicions?  

---

## 2 answers

### Answer by username_2 on 2018-07-04. Score: 3 | Accepted answer

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.

> **Comments:**
> **username_0**: What is meant by the option --limitOutSJcollapsed in the STAR command given in the link you suggested https://github.com/gatk-workflows/gatk3-4-rnaseq-germline-snps-indels/blob/a7a890b38af1309edc040ac89814df37a57f77b7/rna-germline-variant-calling.wdl#L345  
> **username_2**: A 'collapsed' junction is what gets written to the `SJ.out.tab` file during a run of STAR (the collapsing here refers to the removal of duplicate junctions). This option merely limits the maximum number of these junctions that are written. It defaults to 1,000,000.  

---

### Answer by username_3 on 2018-07-13. Score: 4

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](https://software.broadinstitute.org/gatk/documentation/article.php?id=4442) 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.

> **Comments:**
> **username_3**: You might also want to consider the GATK support forum too, as they are very helpful on GATK matters: https://gatkforums.broadinstitute.org/gatk/categories/gatk-support-forum  

---

