HaplotypeCaller doesn't output the VCF

Question

[Flexogore]

Hello everyone!

Still struggling to implement variant calling with my BAM file with the following command:

java -Xmx10g -jar gatk-package-4.1.9.0-local.jar HaplotypeCaller -R illumina.fasta -I final.bam -O illumina.vcf

BAM has been previously indexed and identified:

@SQ SN:GL000222.1   LN:186861   AS:human_g1k_v37_modified.ndx
@SQ SN:GL000200.1   LN:187035   AS:human_g1k_v37_modified.ndx
@SQ SN:GL000193.1   LN:189789   AS:human_g1k_v37_modified.ndx
@SQ SN:GL000194.1   LN:191469   AS:human_g1k_v37_modified.ndx
@SQ SN:GL000225.1   LN:211173   AS:human_g1k_v37_modified.ndx
@SQ SN:GL000192.1   LN:547496   AS:human_g1k_v37_modified.ndx
@RG ID:Sample1  LB:lib1M    PL:illumina SM:1M   PU:unit1
@PG ID:samtools PN:samtools VN:1.10 CL:samtools view -h illumina.initial.bam
@PG ID:samtools.1   PN:samtools PP:samtools VN:1.10 CL:samtools view -bo your_fixed.bam -
@PG ID:samtools.2   PN:samtools PP:samtools.1   VN:1.10 CL:samtools view -H final.bam

The command above seems to be right, however instead of VCF I keep getting such an output listed below (a fragment):

A00128:303:HW5CTDSXX:1:1272:4743:8234/2, A00128:303:HW5CTDSXX:2:2130:26304:21746/2, 
A00128:303:HW5CTDSXX:3:2569:18394:18223/2, A00128:303:HW5CTDSXX:3:1532:10890:10081/1, 
A00128:303:HW5CTDSXX:2:2139:4869:15499/2, A00128:303:HW5CTDSXX:4:2641:19226:19445/2, 
A00128:303:HW5CTDSXX:3:1631:3088:28510/1, A00128:303:HW5CTDSXX:2:2628:13801:23296/2, 
A00128:303:HW5CTDSXX:3:1111:23357:28134/2, A00128:303:HW5CTDSXX:2:1546:11930:7905/2, 
A00128:303:HW5CTDSXX:2:2545:32479:36636/2, A00128:303:HW5CTDSXX:3:2368:16749:30874/2, 
A00128:303:HW5CTDSXX:1:1578:31991:1720/1, A00128:303:HW5CTDSXX:1:1609:25446:6136/1, 
A00128:303:HW5CTDSXX:2:1434:12029:24173/1, A00128:303:HW5CTDSXX:3:2109:28782:36464/2, 
A00128:303:HW5CTDSXX:4:2649:12671:18239/2, A00128:303:HW5CTDSXX:4:2106:27344:22482/1, 
A00128:303:HW5CTDSXX:3:1121:18945:16924/2, A00128:303:HW5CTDSXX:2:2348:4255:23234/2, 
A00128:303:HW5CTDSXX:1:2461:25066:3912/2, A00128:303:HW5CTDSXX:4:2455:13105:34867/2, 
A00128:303:HW5CTDSXX:3:2368:32289:1391/1, A00128:303:HW5CTDSXX:1:1255:14606:10942/2, 
A00128:303:HW5CTDSXX:4:2358:14895:3771/1, A00128:303:HW5CTDSXX:2:1232:21377:21856/1, 
A00128:303:HW5CTDSXX:3:2202:3875:20541/1, A00128:303:HW5CTDSXX:1:1173:16568:12493/2, 
A00128:303:HW5CTDSXX:3:2347:13349:1939/1, A00128:303:HW5CTDSXX:2:1633:20835:5979/2, 
A00128:303:HW5CTDSXX:1:1210:29496:2440/2, A00128:303:HW5CTDSXX:2:1678:6497:10144/2, 
A00128:303:HW5CTDSXX:2:2169:21377:6136/1, A00128:303:HW5CTDSXX:4:2367:31467:29152/1, 
A00128:303:HW5CTDSXX:1:1478:11460:4022/1, A00128:303:HW5CTDSXX:4:2242:6994:16736/2, 
A00128:303:HW5CTDSXX:1:2220:9426:34601/1, A00128:303:HW5CTDSXX:1:2209:12292:6277/1, 
A00128:303:HW5CTDSXX:4:2377:28248:3599/2, A00128:303:HW5CTDSXX:2:2667:18285:36605/1, 
A00128:303:HW5CTDSXX:3:1543:15591:28244/1, A00128:303:HW5CTDSXX:2:1367:31783:14481/1, 
A00128:303:HW5CTDSXX:2:1367:31783:14481/2, A00128:303:HW5CTDSXX:4:2320:27019:28995/1, 
A00128:303:HW5CTDSXX:4:2320:27019:28995/2]
reads contigs = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, X, Y, MT, GL000207.1, GL000226.1, GL000229.1, 
GL000231.1, GL000210.1, GL000239.1, GL000235.1, GL000201.1, GL000247.1, GL000245.1, GL000197.1, GL000203.1, GL000246.1, 
GL000249.1, GL000196.1, GL000248.1, GL000244.1, GL000238.1, GL000202.1, GL000234.1, GL000232.1, GL000206.1, GL000240.1, 
GL000236.1, GL000241.1, GL000243.1, GL000242.1, GL000230.1, GL000237.1, GL000233.1, GL000204.1, GL000198.1, GL000208.1, 
GL000191.1, GL000227.1, GL000228.1, GL000214.1, GL000221.1, GL000209.1, GL000218.1, GL000220.1, GL000213.1, GL000211.1, 
GL000199.1, GL000217.1, GL000216.1, GL000215.1, GL000205.1, GL000219.1, GL000224.1, GL000223.1, GL000195.1, GL000212.1, 
GL000222.1, GL000200.1, GL000193.1, GL000194.1, GL000225.1, GL000192.1]

***********************************************************************
Set the system property GATK_STACKTRACE_ON_USER_EXCEPTION (--java-options '- 
DGATK_STACKTRACE_ON_USER_EXCEPTION=true') to print the stack trace.

I don't get any error messages so I am wondering what's the possible reason. Do you have any ideas?

Answer 1

[Pierre Lindenbaum]

The reference used for mapping the reads is not the same than the one you're using for haplotype caller.

e.g: Incompatible contigs No overlapping contigs found.; GATK BaseRecalibrator error: input files reference and features have incompatible contigs ; error using bqsr and incompatible contigs ; etc... etc...

[Flexogore]

I've used FASTA extracted from BAM via samtools but seems it didn't work. Are there any other options to get the reference without having it itself?

Answer 2

[guyho]

I do the following steps for variant calling (Tool used):

Starting with fastq files

Align to genome with (BWA-mem)

Sort (samtools sort)

Add read groups (AddOrReplaceReadGroups)

mark duplicates (MarkDuplicatesSpark )

base recalibration (BaseRecalibrator)

base recalibration (ApplyBQSR )

Variant calling (HaplotypeCaller)

Integrate vcfs (GenotypeGVCFs )

[Pierre Lindenbaum]

this doesn't answer the question.

[guyho]

HaplotypeCaller assumes sorted bam with read-group tag. If these requirements are not fulfilled it might cause the strange behavior. The question did not specify how the bam was created.