Notice:
If you happen to see a question you know the answer to, please do chime in and help your fellow community members. We encourage our fourm members to be more involved, jump in and help out your fellow researchers with their questions. GATK forum is a community forum and helping each other with using GATK tools and research is the cornerstone of our success as a genomics research community.We appreciate your help!

Test-drive the GATK tools and Best Practices pipelines on Terra


Check out this blog post to learn how you can get started with GATK and try out the pipelines in preconfigured workspaces (with a user-friendly interface!) without having to install anything.

Multisample SNP Calling

qjehanneqjehanne BordeauxMember
edited October 2016 in Ask the GATK team

Hello!

I'm trying to use GATK to call SNPs on multiple bam files (at once) but I encounter some "problems". I tried HaplotypeCaller and UnifiedGenotyper and it works well. The vcf is generated. But not in the way I expected:

CHROM ... FORMAT 20
...
01_contig_678504 ... GT:AD:DP:GQ:PL 1/1:1,249:251:99:11142,712,0
...

I expected to get one genotype per individual, instead of what I only got one sample/genotype named "20". I suppose it's normal and I probably missed something (down below is how I used these tools).

Is it possible to get all different genotypes (for 96 individuals) by snp in one file?

Thanks!

java -jar GenomeAnalysisTK.jar \
-R refs.fa \
-T HaplotypeCaller \
-I bams.list \
-- ERC GVCF \
-o snps.g.vcf

java -jar GenomeAnalysisTK.jar \
-T UnifiedGenotyper \
-R refs.fa -I bams.list \
--genotyping_mode DISCOVERY \
-o snps.raw.vcf -stand_call_conf 25.0 -stand_emit_conf 10.0

Best Answer

Answers

  • valentinvalentin Cambridge, MAMember, Dev ✭✭
    edited October 2016

    I think HaplotypeCaller -ERC GVCF only can handle a single same, so I think that all those bams in bams.list actually make reference to one sample (named "20"?) otherwise the first command would have failed.

    Can you show us an extract of the header of the input bam files? I think the @RG records would be enough. In Unix-like op systems the following command may work assuming that you have samtools installed:

    cat bams.list | xargs -l1 samtools view -H | grep ^@RG

  • qjehanneqjehanne BordeauxMember

    Thanks for your reply!

    So, I looked for headers as you said but I think I misunderstand something. That file, bams.list, is listing (indeed) all my sample like this:
    REALIGN/001.bam
    REALIGN/003.bam
    REALIGN/004.bam
    ...

    In replacement of all these -I in command line no? Else I'm not aware of its other function(s).

    And each bam file starts without any header:
    AS9WV:01464:00111 0 67 12 60 48M1D46M1D65M22S * 0 0 TACACGATGCATTGCCACTCGAAATAAAGATATT
    GGGATTGGACTTATCAATCAATTCAGAACCTTTCGAGCACGACCAATATCTAT
    ...

  • qjehanneqjehanne BordeauxMember

    You're right. Every bam has the same SM tag.

    I'll try to modify that and keep you in touch. Thank you!

  • qjehanneqjehanne BordeauxMember

    Ok, I make it works with different read group tags (obviously...). Thanks a lot!

  • zillurbmb51zillurbmb51 USAMember

    Hello there,
    I am having the same problem. My command: gatk HaplotypeCaller -I bam.list -O snps.hap.g.vcf -R ../../UPR_CANCER_3/all.chr.fa
    The error is: java.lang.RuntimeException: Cannot use index file with textual SAM file
    And my bam headers are :
    (base) [[email protected] analysis]$ samtools view -H ../first_two_plates/BAM/10A_OR_rawlib.bam @HD VN:1.4 SO:coordinate @SQ SN:chr1 LN:249250621 @SQ SN:chr2 LN:243199373 @SQ SN:chr3 LN:198022430 @SQ SN:chr4 LN:191154276 @SQ SN:chr5 LN:180915260 @SQ SN:chr6 LN:171115067 @SQ SN:chr7 LN:159138663 @SQ SN:chr8 LN:146364022 @SQ SN:chr9 LN:141213431 @SQ SN:chr10 LN:135534747 @SQ SN:chr11 LN:135006516 @SQ SN:chr12 LN:133851895 @SQ SN:chr13 LN:115169878 @SQ SN:chr14 LN:107349540 @SQ SN:chr15 LN:102531392 @SQ SN:chr16 LN:90354753 @SQ SN:chr17 LN:81195210 @SQ SN:chr18 LN:78077248 @SQ SN:chr19 LN:59128983 @SQ SN:chr20 LN:63025520 @SQ SN:chr21 LN:48129895 @SQ SN:chr22 LN:51304566 @SQ SN:chrX LN:155270560 @SQ SN:chrY LN:59373566 @SQ SN:chrM LN:16569 @RG ID:10A_OR CN:CGRDeathstarTorrentServer/Chewbacca DS:T0266RD10uprtrainingPlate2pl260ul DT:2013-05-15T03:12:59 FO:TACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACG KS:TCAGTTCCACTTCGCGAT LB:hg19/10A_OR PG:tmap PL:IONTORRENT PU:PGM/316D/10A_OR SM:none @PG ID:bc PN:BaseCaller CL:BaseCaller --trim-qual-cutoff 15 --trim-qual-window-size 30 --trim-adapter-cutoff 16 --calibration-file basecaller_results/recalibration/flowQVtable.txt --phase-estimation-file basecaller_results/recalibration/BaseCaller.json --model-file basecaller_results/recalibration/hpModel.txt --input-dir=sigproc_results --librarykey=TCAG --tfkey=ATCG --run-id=2N1HB --output-dir=basecaller_results --block-col-offset 0 --block-row-offset 0 --datasets=basecaller_results/datasets_pipeline.json --trim-adapter ATCACCGACTGCCCATAGAGAGGCTGAGAC --barcode-filter 0.00 VN:3.6-30/58509 @PG ID:tmap CL:mapall -n 12 -f /results/referenceLibrary/tmap-f3/hg19/hg19.fasta -r basecaller_results/10A_OR_rawlib.basecaller.bam -v -Y -u --prefix-exclude 5 -o 2 stage1 map4 VN:3.6.30 (58509) (201305031917)
    Any help?
    Best Regards
    Zillur

  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    @zillurbmb51

    The program is saying your bam file is actually a sam file.

  • zillurbmb51zillurbmb51 USAMember
    edited June 3

    Thank you very much. It appears that these are indexed bam file (I am not sure why the program is identifying them as sam file). However, after trying several times, I got a new error message and have no idea what to do next. Any help?

    The command:java -jar -Xmx121g /gondor/zillur/tools/GenomeAnalysisTK.jar -R all.chr2.fa -T UnifiedGenotyper -I bam5.list -o snps2.u.vcf

    The error:
    INFO 12:20:35,509 SAMDataSource$SAMReaders - Done initializing BAM readers: total time 0.51 INFO 12:20:36,172 GenomeAnalysisEngine - Preparing for traversal over 316 BAM files INFO 12:20:37,919 GenomeAnalysisEngine - Done preparing for traversal INFO 12:20:37,919 ProgressMeter - [INITIALIZATION COMPLETE; STARTING PROCESSING] INFO 12:20:37,920 ProgressMeter - | processed | time | per 1M | | total | remaining INFO 12:20:37,920 ProgressMeter - Location | sites | elapsed | sites | completed | runtime | runtime INFO 12:20:37,961 StrandBiasTest - SAM/BAM data was found. Attempting to use read data to calculate strand bias annotations values. INFO 12:20:37,962 StrandBiasTest - SAM/BAM data was found. Attempting to use read data to calculate strand bias annotations values. INFO 12:21:07,925 ProgressMeter - chr1:14334217 1.4319616E7 30.0 s 2.0 s 0.5% 108.0 m 107.5 m INFO 12:21:37,926 ProgressMeter - chr1:34490921 3.448832E7 60.0 s 1.0 s 1.1% 89.8 m 88.8 m INFO 12:22:07,928 ProgressMeter - chr1:57275181 5.726208E7 90.0 s 1.0 s 1.9% 81.1 m 79.6 m INFO 12:22:37,930 ProgressMeter - chr1:85581833 8.5573632E7 120.0 s 1.0 s 2.8% 72.3 m 70.3 m INFO 12:23:07,931 ProgressMeter - chr1:115252157 1.15245056E8 2.5 m 1.0 s 3.7% 67.2 m 64.7 m INFO 12:23:37,933 ProgressMeter - chr1:151486381 1.5147008E8 3.0 m 1.0 s 4.9% 61.3 m 58.3 m INFO 12:24:07,934 ProgressMeter - chr1:180585849 1.80584448E8 3.5 m 1.0 s 5.8% 60.0 m 56.5 m INFO 12:24:37,936 ProgressMeter - chr1:210401345 2.10386944E8 4.0 m 1.0 s 6.8% 58.9 m 54.9 m INFO 12:25:07,937 ProgressMeter - chr1:239549581 2.3953408E8 4.5 m 1.0 s 7.7% 58.2 m 53.7 m ##### ERROR -- ##### ERROR stack trace java.lang.ArrayIndexOutOfBoundsException: -2147478782 at java.util.ArrayList.elementData(ArrayList.java:418) at java.util.ArrayList.set(ArrayList.java:446) at org.broadinstitute.gatk.engine.datasources.reads.GATKBAMIndexFromFile.readReferenceSequence(GATKBAMIndexFromFile.java:134) at org.broadinstitute.gatk.engine.datasources.reads.BAMSchedule.<init>(BAMSchedule.java:105) at org.broadinstitute.gatk.engine.datasources.reads.BAMScheduler.getNextOverlappingBAMScheduleEntry(BAMScheduler.java:296) at org.broadinstitute.gatk.engine.datasources.reads.BAMScheduler.advance(BAMScheduler.java:185) at org.broadinstitute.gatk.engine.datasources.reads.BAMScheduler.next(BAMScheduler.java:156) at org.broadinstitute.gatk.engine.datasources.reads.BAMScheduler.next(BAMScheduler.java:46) at htsjdk.samtools.util.PeekableIterator.advance(PeekableIterator.java:68) at htsjdk.samtools.util.PeekableIterator.next(PeekableIterator.java:54) at org.broadinstitute.gatk.engine.datasources.reads.IntervalSharder.next(IntervalSharder.java:79) at org.broadinstitute.gatk.engine.datasources.reads.IntervalSharder.next(IntervalSharder.java:39) at htsjdk.samtools.util.PeekableIterator.advance(PeekableIterator.java:68) at htsjdk.samtools.util.PeekableIterator.next(PeekableIterator.java:54) at org.broadinstitute.gatk.engine.datasources.reads.LocusShardBalancer$1.next(LocusShardBalancer.java:45) at org.broadinstitute.gatk.engine.datasources.reads.LocusShardBalancer$1.next(LocusShardBalancer.java:39) at org.broadinstitute.gatk.engine.executive.LinearMicroScheduler.execute(LinearMicroScheduler.java:89) at org.broadinstitute.gatk.engine.GenomeAnalysisEngine.execute(GenomeAnalysisEngine.java:316) at org.broadinstitute.gatk.engine.CommandLineExecutable.execute(CommandLineExecutable.java:123) at org.broadinstitute.gatk.utils.commandline.CommandLineProgram.start(CommandLineProgram.java:256) at org.broadinstitute.gatk.utils.commandline.CommandLineProgram.start(CommandLineProgram.java:158) at org.broadinstitute.gatk.engine.CommandLineGATK.main(CommandLineGATK.java:108) ##### ERROR ------------------------------------------------------------------------------------------ ##### ERROR A GATK RUNTIME ERROR has occurred (version 3.7-0-gcfedb67): ##### ERROR ##### ERROR This might be a bug. Please check the documentation guide to see if this is a known problem. ##### ERROR If not, please post the error message, with stack trace, to the GATK forum. ##### ERROR Visit our website and forum for extensive documentation and answers to ##### ERROR commonly asked questions https://software.broadinstitute.org/gatk ##### ERROR ##### ERROR MESSAGE: -2147478782 ##### ERROR ------------------------------------------------------------------------------------------

  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    HI @zillurbmb51

    We do not support GATK3 anymore. UnifiedGenotyper is part of GATK3.

    For your purposes you should used Haplotypecaller.

Sign In or Register to comment.