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.

BQSR Warnings: NAs introduced by coercion

micknudsenmicknudsen DenmarkMember ✭✭

Hi,

When running the BQSR script, I get the following warnings

Warning messages:
1: NAs introduced by coercion
2: NAs introduced by coercion

I have managed to track down exactly where they come from:

for(cov in levels(data$CovariateName)) {
  d = data[data$CovariateName==cov,]
  if( cov == "Context" ) {
    d$CovariateValue = as.character(d$CovariateValue)
    d$CovariateValue = substring(d$CovariateValue,nchar(d$CovariateValue)-2,nchar(d$CovariateValue))
  } else {
    d$CovariateValue = as.numeric(levels(d$CovariateValue))[as.integer(d$CovariateValue)]
  }

Here the problem is that levels(d$CovariateValue) contains both integers and strings (short DNA sequences), and the latter causes as.numeric to introduce NAs.

Is this something to be worried about? I am using GATK 3.4-46, but the error also occurs in 3.3-0.

Thanks,
Michael Knudsen

Comments

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @micknudsen
    Hi Michael Knudsen,

    Can you tell me a little more about the data you are working with? Are you following our Best Practices? Can you post the exact commands you ran for Base Recalibrator step?

    Thanks,
    Sheila

  • micknudsenmicknudsen DenmarkMember ✭✭
    edited August 2015

    Hi Sheila,

    I run gatkMem using the following settings:

    gatkMem="java -Djava.awt.headless=true  -Xmx60g  -jar $gatkdir/../jar-bin/GenomeAnalysisTK.jar -et NO_ET -K [PATH_TO_MY_KEY_FILE]
    

    The analysis is the performed as follows (names in curly braces representing actual files):

    $gatkMem -T AnalyzeCovariates -R {reference_genome} -before {recalibrated_data_table} -after {post_recalibrated_data_table} -plots {recalibrated_plots} -l DEBUG
    

    I do get an output plot, but I am curious about the coercion error. It comes from the fact that the csv file (excerpt below) contains non-numeric entries in the CovariateName column. This causes as.numeric(levels(d$CovariateValue)) in the R script to return NA values.

    ReadGroup,CovariateValue,CovariateName,EventType,Observations,Errors,EmpiricalQuality,AverageReportedQuality,Accuracy,Recalibration
    H129LBGXX.1,AA,Context,Base Substitution,152573,9413.58,12.00,14.13,-2.13,Before
    H129LBGXX.1,AAA,Context,Base Insertion,56627,11.17,44.00,45.00,-1.00,Before
    H129LBGXX.1,AAA,Context,Base Deletion,56627,17.48,44.00,45.00,-1.00,Before
    H129LBGXX.1,CA,Context,Base Substitution,110720,10661.52,10.00,13.87,-3.87,Before
    H129LBGXX.1,CAA,Context,Base Insertion,27630,2.01,45.00,45.00,0.00,Before
    (...)
    H129LBGXX.1,1,Cycle,Base Substitution,5813,102.22,18.00,18.94,-0.94,Before
    H129LBGXX.1,-1,Cycle,Base Substitution,5616,114.83,17.00,16.84,0.16,Before
    H129LBGXX.1,2,Cycle,Base Substitution,5836,136.43,17.00,18.72,-1.72,Before
    H129LBGXX.1,-2,Cycle,Base Substitution,5629,133.51,16.00,16.86,-0.86,Before
    H129LBGXX.1,3,Cycle,Base Substitution,5870,186.58,16.00,18.86,-2.86,Before
    H129LBGXX.1,-3,Cycle,Base Substitution,5656,153.00,16.00,17.39,-1.39,Before
    (...)

    /Michael

  • micknudsenmicknudsen DenmarkMember ✭✭
    edited August 2015

    For some reason, the system will not render my reply. Below is a screenshot. Sorry!

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin
    edited August 2015

    @micknudsen
    Hi Michael,

    What is GATKMem? I have not heard of this!
    Also, what was your command for BaseRecalibrator?

    Thanks,
    Sheila

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @micknudsen
    Okay. Your csv file is missing the QualityScore column. It looks like the CovariateName column looks fine in your csv file. There is an example csv file here: http://gatkforums.broadinstitute.org/discussion/44/base-quality-score-recalibration-bqsr

    I'm not sure why your csv is missing the QualityScores. What kind of data are you working with?

  • micknudsenmicknudsen DenmarkMember ✭✭

    @Sheila said:

    What is GATKMem? I have not heard of this!

    GATKMem is just an alias forthe command:

    java -Djava.awt.headless=true -Xmx60g -jar $gatkdir/../jar-bin/GenomeAnalysisTK.jar -et NO_ET -K $gatkdir/../foo.key

    Also, what was your command for BaseRecalibrator?

    I have run three iterations of BaseRecalibrator. Things in curly braces represent files.

    $gatkMem -T BaseRecalibrator -L {target_intervals} --interval_padding 100 -I {realigned_bam_file} -R {reference_genome} -knownSites {dbsnp} -knownSites {indel_calls} -knownSites {g1k} -o {recalibrated_data_table} -nct {cores}
            $gatkMem -T BaseRecalibrator -I {realigned_bam_file} -R {reference_genome} -knownSites {dbsnp} -knownSites {indel_calls} -knownSites {g1k} -BQSR {recalibrated_data_table} -o {post_recalibrated_data_table} -nct {cores}
            $gatkMem -T BaseRecalibrator -L {target_intervals} --interval_padding 100 -I {realigned_bam_file} -R {reference_genome} -knownSites {dbsnp} -knownSites {indel_calls} -knownSites {g1k} -BQSR {recalibrated_data_table} -o {post_recalibrated_data_table} -nct {cores}
    

    /Michael

  • micknudsenmicknudsen DenmarkMember ✭✭

    @Sheila said:
    What is GATKMem? I have not heard of this!

    GATKMem is just an alias for:

    java -Djava.awt.headless=true -Xmx60g -jar $gatkdir/../jar-bin/GenomeAnalysisTK.jar -et NO_ET -K $gatkdir/../foo.key

    Also, what was your command for BaseRecalibrator?

    I have run BaseRecalibrator three times. The names in curly braces refer to files.

    $gatkMem -T BaseRecalibrator -L {target_intervals} --interval_padding 100 -I {realigned_bam_file} -R {reference_genome} -knownSites {dbsnp} -knownSites {indel_calls} -knownSites {g1k} -o {recalibrated_data_table} -nct {cores}
    
    $gatkMem -T BaseRecalibrator -I {realigned_bam_file} -R {reference_genome} -knownSites {dbsnp} -knownSites {indel_calls} -knownSites {g1k} -BQSR {recalibrated_data_table} -o {post_recalibrated_data_table} -nct {cores}
    
    $gatkMem -T BaseRecalibrator -L {target_intervals} --interval_padding 100 -I {realigned_bam_file} -R {reference_genome} -knownSites {dbsnp} -knownSites {indel_calls} -knownSites {g1k} -BQSR {recalibrated_data_table} -o {post_recalibrated_data_table} -nct {cores}
    

    /Michael

  • micknudsenmicknudsen DenmarkMember ✭✭

    This is really odd. Each time I post a reply, the content disappears. I don't know what else to do but post a screenshot of my comment before hitting the 'Post Comment' button...

  • micknudsenmicknudsen DenmarkMember ✭✭

    @Sheila said:
    I'm not sure why your csv is missing the QualityScores. What kind of data are you working with?

    I am working with paired-end reads from MiSeq whole-exom sequencing.

    /Michael

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @micknudsen
    Hi Michael,

    I just verified your account, so you should be able to post properly here now. As for the issue, I have not seen this before. Is it possible that you have too few bases to run BaseRecalibrator? We recommend having a minimum of 100M bases. http://gatkforums.broadinstitute.org/discussion/4272/targeted-sequencing-appropriate-to-use-baserecalibrator-bqsr-on-150m-bases-over-small-intervals

    -Sheila

Sign In or Register to comment.