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.
Attention:
We will be out of the office on November 11th and 13th 2019, due to the U.S. holiday(Veteran's day) and due to a team event(Nov 13th). We will return to monitoring the GATK forum on November 12th and 14th respectively. Thank you for your patience.

Variant calling of PacBio data using UnifiedGenotyper or HaplotypeCaller

micschmicsch SwitzerlandMember
edited July 2015 in Ask the GATK team

We have PacBio data where we want to do variant calling. I tried both UnifiedGenotyper and Haplotype caller. I was not very successfull doing that. When I used UnifiedGenotyper I got some output, but just SNPs but NO indels... (I skipped the realigment part there). I tried to play arround with the parameters (indelGapContinuationPenalty, indGapOpenPenalty, min_base_quality_score). Setting the "min_base_quality_score" to a lower value is giving at least this output with only SNPs as mentioned above.

The only manual for PacBio data on GATK I got was this: https://www.broadinstitute.org/gatk/guide/topic?name=methods
But this document is pretty old. Are there newer developments regarding GATK for PacBio? Or any more detailed tutorials?

Or would you suggest PBHoney as the right tool to use? Anything else?

Just to mention: For Illumina your toolkit worked like a charm! So basically we are able to work with it...

Thanks!
Michael

Answers

  • SheilaSheila Broad InstituteMember, Broadie admin

    @micsch
    Hi Michael,

    What happens when you use Haplotype Caller? Do you get indel calls? Did you follow the pre-processing steps in the article you posted above?

    Thanks,
    Sheila

  • bayolaubayolau earthMember

    Hello @Sheila,

    I am very interested in gatk with pacbio reads also. I have successfully called SNPs with recalibration+unified genotyper, but no INDELS can be called so far, even with -stand_emit_conf 0.1 -stand_call_conf 0.1 -mbq 1 and various settings of -minIndelCnt, -indelGCP, and -indelGOP.

    I tried using the same recalibrate BAM with haplotype caller, but the out VCF is empty. The -APO output does signify some non-zero activity on the fifth column and the -graph output has something. I wonder if there's a guide to decipher the -graph output, or are there any known settings/examples of haplotype caller with PacBio reads?

  • SheilaSheila Broad InstituteMember, Broadie admin

    @bayolau
    Hi,

    Oh, so you are saying that with Unified Genotyper, SNPs are called, but with Haplotype Caller no SNPs or indels are called?

    Would you mind posting the -APO output? This article may help with bamout files: http://gatkforums.broadinstitute.org/discussion/5484/howto-generate-a-bamout-file-showing-how-haplotypecaller-has-remapped-sequence-reads

    Unfortunately, it looks like we don't have anything else on PacBio data. The testing is done on Illumina reads here. I am quite curious to learn more about your findings.

    -Sheila

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Should check what the qualities are like, both base and mapping. The original post suggests that lowering min_base_quality_score was necessary to call SNPs with HC, which is a bad sign. HC applies very stringent filters.

  • evolvedmicrobeevolvedmicrobe MGHMember

    As @Sheila mentioned, the GATK is really setup for Illumina data and not PacBio data, so using the GATK right out of the box without tweaking it is going to create problems (and will be a sub-optimal use of the PacBio data). Also, are you using raw reads or CCS/ROI reads?

    If you want to try it though, I would suggest the following. First off, for the UnifiedGenotyper, by default SNP is the default option, so to call indels you should use the flag for both, "-glm BOTH"

    Critically, the default error modes for the GATK are primed for Illumina data, so quality score recalibration will be essential. For the Illumina data the GATK typically sees, the QV values represent the probability of a substitution error and there is no direct information on the probability of an insertion or deletion. In contrast, in PacBio data, there are essentially no substitution errors, but there are more insertion and deletion errors (and rather than a single QV they report a rich set of QV values for each error type). The key feature in PacBio data is that these errors are unbiased and so you can just stack reads on top of each other and converge on the right result. However, this consensus generating process will not work very well if you don't have a statistical model that appropriately handles indels and their relative probabilities. It would therefore be essential that the indel probabilities are calculated correctly. It is possible that the current version of the GATK Base Quality Score Recalibration can do this, but it's an experiment you'll have to try. Feel free to get in touch with your local PacBio field application specialist for additional advice on calling variants.

  • bayolaubayolau earthMember
    edited July 2015

    @Sheila thank you. I'll soon have more leisure time to get back to the exercise. If I do cat gatk.activityprofile | awk '$5>0' from HC's APO, I'll get lines like

    Chromosome  Start   End Feature ActivityProfile
    21  15221368    15221369    state   0.91774
    21  15221369    15221370    state   0.91588
    21  15221370    15221371    state   0.91377
    21  15221371    15221372    state   0.91145
    21  15221372    15221373    state   0.90892
    21  15221373    15221374    state   0.90588
    21  15221374    15221375    state   0.90287
    21  15221375    15221376    state   0.89973
    21  15221376    15221377    state   0.89644
    21  15221377    15221378    state   0.89302
    21  15221378    15221379    state   0.88975
    21  15221379    15221380    state   0.88635
    21  15221380    15221381    state   0.88252
    

    @Geraldine_VdAuwera @evolvedmicrobe The mapping qualities are 60. The recalibration gives much higher quality scores

    #:GATKTable:6:3:%s:%s:%.4f:%.4f:%d:%.2f:;
    #:GATKTable:RecalTable0:
    ReadGroup  EventType  EmpiricalQuality  EstimatedQReported  Observations  Errors
    hole       M                   16.0000             10.5585       8841541  228618.23
    hole       I                   26.0000             45.0000       8841541  339355.45
    hole       D                   26.0000             45.0000       8841541  341969.02
    
    #:GATKTable:6:11:%s:%s:%s:%.4f:%d:%.2f:;
    #:GATKTable:RecalTable1:
    ReadGroup  QualityScore  EventType  EmpiricalQuality  Observations  Errors
    hole                  6  M                   14.0000        440103   18474.29
    hole                  7  M                   14.0000        420844   16697.06
    hole                  8  M                   15.0000        441856   15404.89
    hole                  9  M                   15.0000        572300   17710.82
    hole                 10  M                   16.0000        841074   22521.64
    hole                 11  M                   16.0000       1199485   29707.48
    hole                 12  M                   16.0000       1648572   39387.72
    hole                 13  M                   17.0000       2618147   56168.16
    hole                 14  M                   17.0000        659160   12546.18
    hole                 45  I                   26.0000       8841541  339355.45
    hole                 45  D                   26.0000       8841541  341969.02
    

    For the successful SNP-calling UG runs, I've tried for example -stand_emit_conf 0.1 -stand_call_conf 0.1 -mbq 1 -baq CALCULATE_AS_NECESSARY --genotype_likelihoods_model BOTH --defaultBaseQualities 9. I've also tried things in the line of -deletions 0.5 -minIndelCnt 2 -indelGCP 0 -indelGOP 0 to get indels (with no success yet)

    EDIT: @evolvedmicrobe these are raw reads

  • SheilaSheila Broad InstituteMember, Broadie admin

    @bayolau
    Hi,

    Can you post a bam file screenshot and bamout file screenshot of a region you think should have some variants called?
    Also, please do check that the base qualities are good in each region as well.

    Thanks,
    Sheila

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @bayolau where did you get the commands that you are running? Those are using a tool (UG) and arguments that have been deprecated for several years. You should really use our latest recommendations and tools (HaplotypeCaller).

Sign In or Register to comment.