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!

Sign In or Register to comment.