BQSR Warnings: NAs introduced by coercion

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.