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 ?