Base recalibration & realignment without known sites

gilgigilgi Posts: 19Member
edited October 2012 in Ask the GATK team

Dear GATK team,

Thanks a lot for the new GATK version and GATK forum!

I am trying to use GATK for yeast strains. I do not have files of known sites of SNPs/indels.
I understand that the BaseRecalibrator must get such a file.
Do you suggest to skip calibration and realignment, or is there another way to go here?

Post edited by Geraldine_VdAuwera on

Best Answer

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 8,169Administrator, GATK Dev admin
    edited September 2012 Answer ✓

    Hi Gilgi,

    Glad you like the new digs!

    Skipping Base recalibration and realignment are never a good option. Theere is an article on Base Quality Score Recalibration in the Methods & Workflows section of the Guide that you need to read. Towards the end, in the "Troubleshooting" section, there is a discussion of what to do if you don't have a set of known variants for your organism.

    Good luck!

    Post edited by Geraldine_VdAuwera on

    Geraldine Van der Auwera, PhD

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 8,169Administrator, GATK Dev admin
    edited September 2012 Answer ✓

    Hi Gilgi,

    Glad you like the new digs!

    Skipping Base recalibration and realignment are never a good option. Theere is an article on Base Quality Score Recalibration in the Methods & Workflows section of the Guide that you need to read. Towards the end, in the "Troubleshooting" section, there is a discussion of what to do if you don't have a set of known variants for your organism.

    Good luck!

    Post edited by Geraldine_VdAuwera on

    Geraldine Van der Auwera, PhD

  • gilgigilgi Posts: 19Member

    Thanks a lot!!!

  • FrankFrank Posts: 2Member

    The link is broken right now.

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 8,169Administrator, GATK Dev admin

    Link fixed, sorry about that.

    Geraldine Van der Auwera, PhD

  • mhtmht Posts: 8Member

    Hi, do you have any solution for the IndelRealigner as well? I have a non-model organism as reference and there are no known indel sites/document which I can use.

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 8,169Administrator, GATK Dev admin

    If you don't have any known indels you can run without.

    Geraldine Van der Auwera, PhD

  • mhtmht Posts: 8Member

    Hi, I ran it without any known indels, but some regions didn't seem to realign as expected. As a result, GATK is calling both a SNP and single base insertion at the same location (the true variant is the insert). Any idea why this is happening?

  • ebanksebanks Broad InstitutePosts: 689Member, Administrator, GATK Dev, Broadie, Moderator, DSDE Dev, GP Member admin

    What's your command line?

    Eric Banks, PhD -- Senior Group Leader, MPG Analysis, Broad Institute of Harvard and MIT

  • mhtmht Posts: 8Member

    Hi Eric,

    to identify intervals:
    java -Xmx1g -jar GenomeAnalysisTKLite.jar -T RealignerTargetCreator -R reference.fa -I sampleA_clean_paired.sorted.dedup.bam -o sampleA_clean_paired.sorted.dedup.intervals

    to realign:
    java -Xmx2g -jar GenomeAnalysisTKLite.jar -T IndelRealigner -R reference.fa -I sampleA_clean_paired.sorted.dedup.bam -targetIntervals sampleA_clean_paired.sorted.dedup.intervals -o sampleA_clean_paired.sorted.dedup.realigned.bam

    Thanks for your help

  • ebanksebanks Broad InstitutePosts: 689Member, Administrator, GATK Dev, Broadie, Moderator, DSDE Dev, GP Member admin

    Okay, good. Those look reasonable.
    Now why do you think some regions didn't realign properly? Could you please post an example IGV screenshot with an explanation? Thanks.

    Eric Banks, PhD -- Senior Group Leader, MPG Analysis, Broad Institute of Harvard and MIT

  • mhtmht Posts: 8Member

    Hi Eric,

    I've attached a screenshot of a region as an example. (see: http://imageshack.us/photo/my-images/803/realigned.png)
    Also, I noticed some messages during running GATK that goes like this:

    INFO 15:56:32,707 IndelRealigner - Not attempting realignment in interval contig_31331:6142-6275 because there are too many reads.

    I'm wondering if this is the reason why... and if so, any suggestions on how to deal with this?

    Thanks so much for your help.

  • ebanksebanks Broad InstitutePosts: 689Member, Administrator, GATK Dev, Broadie, Moderator, DSDE Dev, GP Member admin

    That is certainly a possible reason. Please take a look at the documentation for command-line arguments that allow you to increase the maximum number of reads allowed in a given realignment region.

    Eric Banks, PhD -- Senior Group Leader, MPG Analysis, Broad Institute of Harvard and MIT

  • mhtmht Posts: 8Member

    Hi Eric,

    I found the --maxReadsForRealignment option and I set it at a higher number (~ 100000) and am getting less of those messages. There are, however, still a few of these cases. There is a note in the documentation which says "For expert users only!". Is there other harms of increasing this number, other than the toll it takes on memory? How do I decide what to set this number at?

    Much thanks.

  • ebanksebanks Broad InstitutePosts: 689Member, Administrator, GATK Dev, Broadie, Moderator, DSDE Dev, GP Member admin

    Just memory. The value to use really depends on your average coverage and whether it's important to realign every single region that might need it.

    Eric Banks, PhD -- Senior Group Leader, MPG Analysis, Broad Institute of Harvard and MIT

  • nanacnanac Posts: 9Member

    Hi, I would like to use the BaseRecalibrator but I do not have a database of known SNPs for the organism I am studying. I have a VCF files with potential real SNPs and I want to use the best ones as known sites. How many SNPS at least do we need? Are 500 enough?

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 8,169Administrator, GATK Dev admin

    There is no hard number, it depends how divergent you expect your organism samples to be. Do you have any way to estimate this?

    Geraldine Van der Auwera, PhD

  • nanacnanac Posts: 9Member
    edited November 2012

    We will study 40 genotypes. So far, the genetic differentiation between 2 genotypes has been assessed from microsatellite markers (Fsc), and varies from 50% to 70%. We currently have GATK results on 7 genotypes from a previous study reporting at least 14500 polymorphic sites present on each genotype with a quality between 1329 and 4867. Maybe we can use these ones or a part of them?

    Post edited by nanac on
  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 8,169Administrator, GATK Dev admin

    Yes, you can use a high-confidence subset of those.

    Geraldine Van der Auwera, PhD

  • joanajoana Posts: 2Member

    Hi
    I also work on nonmodel organisms and don't have a set of known polymorphisms. I have 1-4 individuals each from several quite divergent species. Would you use high-confidence SNPs found with Unified Genotyper run over all individuals of all species together or would you define the set of high-confidence SNPs for each species separately, even though I only have few individuals each?
    Thanks a lot.

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 8,169Administrator, GATK Dev admin

    Hi Joana,

    If the species are quite divergent I would probably generate separate high-confidence sets for each. It may be interesting to then look at how the separate sets intersect, how much overlap there is.

    Geraldine Van der Auwera, PhD

Sign In or Register to comment.