GATK AlternateReferenceMaker in non-model organism & SNP flanking regions

TammyTammy Washington, DCMember

Hi everyone,

I have a question about using the GATK AlternateReferenceMaker and how it may affect down stream analysis. I apologize for the long post but I tried to include only relevant info. Any thoughts on analysis are helpful!

Here is basically what I did:

  • Illumina sequencing of 300 libraries with methods similar to genotype-by-sequencing (like RAD-tag libraries without shearing)
  • Created contigs of sequences within individuals (cap3) and compared across 300 individuals
  • Selected contigs found in at least 5 individuals and mapped to a draft reference genome of sister species (bwa-mem)
  • Used GATK AlternateReferenceMaker to insert differences and create a "pseudo"-reference for my species

  • Then I have did SNP calling from raw sequences using GATK best practices (alignment of each individual to pseudo-reference, realignment, g.vcf, Haplotype Caller)

  • I did hard filtering for quality, coverage, and missing data across SNPs
  • I have sub-selected SNPs with minor allele frequencies > 0.25 to use in designing a SNP array for parentage (~800 SNPs)
  • When I try to select the flanking regions around these SNPs using bedtools, it pulls the flanking regions from the pseudo-reference and many of them have a lot of N's... so my questions are:
  1. How can I tell if these were produced from the process I used to create a "pseudo"-reference for my species? Could they be misalignments from sequence contigs to reference or between species? How would this affect the SNP calling proceess?

  2. What is the best way to ground truth the sequencing region? Should I try to pull the alignments from the region around each SNP for each individual and look at them manually? How can I do this?

  3. Finally, the sister species and my pseudo-reference are not indexed in the same way. I'm guessing I should've sorted them before indexing or specified the indexing to use when using the GATK AlternateReferenceMaker. Can I convert the indexing?

Thanks for your thoughts!

Best Answers


  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    Hi Tammy,

    1) Were there N's in your original reference at those positions? Were there N's in the sample contigs at those positions? Can you post the exact command line you ran for FastaAlternateReferenceMaker? GATK will skip over N bases.

    2) Yes, you can do that by checking both the original reference and the newly created contigs against the reads in IGV.

    3) What do you mean "not indexed in the same way"?


  • TammyTammy Washington, DCMember

    Hi Sheila,

    Here are the commands I used:
    java -Xmx64g -jar ~/Programs/GenomeAnalysisTK.jar -R ~/Catharus_ustulatus_draft_assembly_Feb13_2015.fasta -T RealignerTargetCreator -DBQ 20 -I woth-goodcontigs-clean-sort.bam -o woth-goodcontigs-clean-sort.bam.indelrealigner.intervals

    java -Xmx64g -jar ~/Programs/GenomeAnalysisTK.jar -R ~/Catharus_ustulatus_draft_assembly_Feb13_2015.fasta -T IndelRealigner -DBQ 20 -I woth-goodcontigs-clean-sort.bam -targetIntervals woth-goodcontigs-clean-sort.bam.indelrealigner.intervals -o woth-goodcontigs-clean-sort-realigned.bam

    java -Xmx64g -jar ~/Programs/GenomeAnalysisTK.jar -R ~/Catharus_ustulatus_draft_assembly_Feb13_2015.fasta -T HaplotypeCaller -I woth-goodcontigs-clean-sort-realigned.bam -o woth-goodcontigs-snps.indels.raw.vcf

    java -Xmx64g -jar ~/Programs/GenomeAnalysisTK.jar -R ~/Catharus_ustulatus_draft_assembly_Feb13_2015.fasta -T FastaAlternateReferenceMaker -o woth-vcf-reference.fasta -V woth-goodcontigs-snps.indels.raw.vcf

    I was trying to say that the chromosome names and lengths are not the same between the two genomes. I'm guessing they are in the same order. Looking at the fai files:

    Catharus ustulatus draft genome:
    16 8833 4 60 61
    LG5 13109 8990 60 61
    LG2 109741 22323 60 61
    1B_random 142709 133905 60 61
    new woth reference made with GATK alternate reference makers:
    1 8833 3 60 61
    2 13109 8987 60 61
    3 109742 22318 60 61
    4 142709 133893 60 61

    So if I have genomic intervals from the alternate genome I created with GATK, how do I get the appropriate intervals from the original draft genome?

    Therefore I have only pulled out the intervals from the pseduo-reference I created. I need to figure this out so I can look at the contigs and original genome to see where the N's came from.

    You said I could pull out the aligned sequences from all the individuals. I have the vcf file, the bed file, and the 600 bam files used to create them... is there an easy way to extract all the alignments for all individuals for one SNP? (I have been having trouble running IGV on my server which is running Debian. But I could examine each SNP region individually if I can pull out all of the correct alignments)

    I'm sorry if these questions are simple; I am new to genomic datasets!


    Issue · Github
    by Sheila

    Issue Number
    Last Updated
    Closed By
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi Tammy,

    Sorry for the delayed response; we had to discuss this internally. Unfortunately this kind of question is largely out of scope of the support we can provide, since it's a fairly complex experimental design that we don't have experience with, so the amount of help we can give you is going to be limited.

    I can tell you however that one way to pull out the regions of interest is to use the PrintReads tool with the -L argument to specify intervals of interest.

    Also, I recommend you make sure to use the latest version of GATK, because some improvements were made recently in how the alternate reference making tools name output intervals. It used to be that we simply named them 1, 2, 3 etc in the order in which they were created; now the name includes a reference to the original contig names.

    Good luck!

  • TammyTammy Washington, DCMember

    Hi Geraldine,

    Thanks for your time and thoughts on this question. I will try sub setting my data with the PrintReads tool. I am still working on comparing the genomes to find out what files included these ambiguous bases.

    I think it would be difficult at this point to redo the alternate reference with the latest version of GATK, since all of the downstream analysis would need to be redone. When you say that the output intervals were named in the order they were created, does that mean the same order as the original reference?


  • TammyTammy Washington, DCMember

    Great! Thank you.

Sign In or Register to comment.