HaplotypeCaller doesn't output the VCF
Entering edit mode
13 months ago
Flexogore ▴ 10

Hello everyone!

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

java -Xmx10g -jar gatk-package- 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, 
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?

BAM VCF GATK HaplotypeCaller • 595 views
Entering edit mode
13 months ago

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

Entering edit mode

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?

Entering edit mode
13 months ago
guyho • 0

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 )

Entering edit mode

this doesn't answer the question.

Entering edit mode

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.


Login before adding your answer.

Traffic: 1629 users visited in the last hour
Help About
Access RSS

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6