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.

AnalyzeCovariates error

estif74estif74 Saint Paul, MN, USAMember

Hi there, I'm trying to get AnalyzeCovariates to work, and I get the following error message...any thoughts as to what is going wrong? Not sure whether this is a program-related issue or not?

Thanks!!

docsmb17:No_Backup sgfriede$ java -jar /Users/sgfriede/GeneApps/GenomeAnalysisTK-3.1.1/GenomeAnalysisTK.jar -T AnalyzeCovariates \ -R canFam3.fa \ -before recal.table \ -after after_recal.table \ -plots Mariah_recal_plots.pdf
INFO 09:00:35,945 HelpFormatter - --------------------------------------------------------------------------------
INFO 09:00:35,947 HelpFormatter - The Genome Analysis Toolkit (GATK) v3.1-1-g07a4bf8, Compiled 2014/03/18 06:09:21
INFO 09:00:35,947 HelpFormatter - Copyright (c) 2010 The Broad Institute
INFO 09:00:35,948 HelpFormatter - For support and documentation go to http://www.broadinstitute.org/gatk
INFO 09:00:35,951 HelpFormatter - Program Args: -T AnalyzeCovariates -R canFam3.fa -before recal.table -after after_recal.table -plots Mariah_recal_plots.pdf
INFO 09:00:36,278 HelpFormatter - Executing as [email protected] on Mac OS X 10.9.3 x86_64; Java HotSpot(TM) 64-Bit Server VM 1.7.0_55-b13.
INFO 09:00:36,279 HelpFormatter - Date/Time: 2014/05/26 09:00:35
INFO 09:00:36,279 HelpFormatter - --------------------------------------------------------------------------------
INFO 09:00:36,279 HelpFormatter - --------------------------------------------------------------------------------
INFO 09:00:36,706 GenomeAnalysisEngine - Strictness is SILENT
INFO 09:00:37,005 GenomeAnalysisEngine - Downsampling Settings: Method: BY_SAMPLE, Target Coverage: 1000
INFO 09:00:37,163 GenomeAnalysisEngine - Preparing for traversal
INFO 09:00:37,179 GenomeAnalysisEngine - Done preparing for traversal
INFO 09:00:37,179 ProgressMeter - [INITIALIZATION COMPLETE; STARTING PROCESSING]
INFO 09:00:37,179 ProgressMeter - Location processed.sites runtime per.1M.sites completed total.runtime remaining
INFO 09:00:37,514 ContextCovariate - Context sizes: base substitution model 2, indel substitution model 3
INFO 09:00:37,802 ContextCovariate - Context sizes: base substitution model 2, indel substitution model 3
INFO 09:00:37,852 AnalyzeCovariates - Generating csv file '/var/folders/jx/3t4ylccn0g32xwl0qs6y4tz80000gp/T/AnalyzeCovariates3417596129729096429.csv'
INFO 09:00:38,256 AnalyzeCovariates - Generating plots file 'Mariah_recal_plots.pdf'
INFO 09:00:39,706 GATKRunReport - Uploaded run statistics report to AWS S3

ERROR ------------------------------------------------------------------------------------------
ERROR stack trace

org.broadinstitute.sting.utils.R.RScriptExecutorException: RScript exited with 1. Run with -l DEBUG for more info.
at org.broadinstitute.sting.utils.R.RScriptExecutor.exec(RScriptExecutor.java:174)
at org.broadinstitute.sting.utils.recalibration.RecalUtils.generatePlots(RecalUtils.java:548)
at org.broadinstitute.sting.gatk.walkers.bqsr.AnalyzeCovariates.generatePlots(AnalyzeCovariates.java:380)
at org.broadinstitute.sting.gatk.walkers.bqsr.AnalyzeCovariates.initialize(AnalyzeCovariates.java:394)
at org.broadinstitute.sting.gatk.executive.LinearMicroScheduler.execute(LinearMicroScheduler.java:83)
at org.broadinstitute.sting.gatk.GenomeAnalysisEngine.execute(GenomeAnalysisEngine.java:313)
at org.broadinstitute.sting.gatk.CommandLineExecutable.execute(CommandLineExecutable.java:121)
at org.broadinstitute.sting.commandline.CommandLineProgram.start(CommandLineProgram.java:248)
at org.broadinstitute.sting.commandline.CommandLineProgram.start(CommandLineProgram.java:155)
at org.broadinstitute.sting.gatk.CommandLineGATK.main(CommandLineGATK.java:107)

ERROR ------------------------------------------------------------------------------------------
ERROR A GATK RUNTIME ERROR has occurred (version 3.1-1-g07a4bf8):
ERROR
ERROR This might be a bug. Please check the documentation guide to see if this is a known problem.
ERROR If not, please post the error message, with stack trace, to the GATK forum.
ERROR Visit our website and forum for extensive documentation and answers to
ERROR commonly asked questions http://www.broadinstitute.org/gatk
ERROR
ERROR MESSAGE: RScript exited with 1. Run with -l DEBUG for more info.
ERROR ------------------------------------------------------------------------------------------

Best Answer

Answers

  • SheilaSheila admin Broad InstituteMember, Broadie, Moderator admin
    edited May 2014

    @estif74‌

    Hi,

    It looks like you might not have the R libraries that the script depends on installed. Please refer to http://gatkforums.broadinstitute.org/discussion/1244/what-is-a-gatkreport. This tells you how to install the libraries you need.

    -Sheila

  • estif74estif74 Saint Paul, MN, USAMember

    Hi Sheila - didn't see anything posted here? Thanks!

  • SheilaSheila admin Broad InstituteMember, Broadie, Moderator admin

    @estif74‌

    Hi,

    Sorry about that; I pressed enter too fast! There is an answer now.

    -Sheila

  • estif74estif74 Saint Paul, MN, USAMember

    @Sheila‌
    Thanks so much for your help. I've done this already, and even changed my .Rprofile to include:

    .libPaths("/Library/Frameworks/R.framework/Resources/library/")

    But I still get the same error message. Any other thoughts? I'm sure it's a problem with my R configuration but I'm a bit of a novice at all of this...

    Thanks a ton.

  • SheilaSheila admin Broad InstituteMember, Broadie, Moderator admin
    edited May 2014

    @estif74‌

    Hi,

    You can run this in R directly instead of running it in GATK. This will tell you exactly what the error is.

    1) You can get the source code (script itself) Find this here: https://github.com/broadgsa/gatk-protected

    2) Identify the parameters needed to run it (run with -l DEBUG to get parameters)

    3) Run in R and see what error message is

    -Sheila

  • estif74estif74 Saint Paul, MN, USAMember

    @Sheila‌
    Hi Sheila, thanks for this input, but being a bit of a newbie - is there any way you can give me some more specifics (or point me to a good resource) to learn how to run the source code directly in R?

    Also, is the source you're suggesting I run found here:
    gatk-protected / protected / gatk-protected / src / main / java / org / broadinstitute / sting / gatk / walkers / bqsr /
    I just want to make sure I'm downloading the right one.

    Thank you again!
    Steve

  • Geraldine_VdAuweraGeraldine_VdAuwera admin Cambridge, MAMember, Administrator, Broadie admin

    Hi Steve (@estif74),

    The script you want is here: https://github.com/broadgsa/gatk-protected/blob/master/public/gatk-framework/src/main/resources/org/broadinstitute/sting/utils/recalibration/BQSR.R

    Here's what you need to do:

    1. Re-run AnalyzeCovariates with these additional parameters: -l DEBUG (that's a lowercase L, not an uppercase i, to be clear) and -csv my-report.csv (where you can call the file anything; this is so the intermediate csv file will be saved).

    2. Identify the line in the log output that says what parameters the RScript is given. If you can't find it, post the log output here and we will show you how to find it.

    3. Run the script manually with those arguments. For newbies the easiest way to do this is to do it from within an IDE program like RStudio. Or you can start up R at the command line and run it that way, whatever you're comfortable with.

  • estif74estif74 Saint Paul, MN, USAMember

    @Geraldine_VdAuwera‌
    Thanks so much for your help. When I look a the log output, this seems to be where the error occurs:


    Loading required package: methods
    KernSmooth 2.23 loaded
    Copyright M. P. Wand 1997-2009

    Attaching package: ‘gplots’

    The following object is masked from ‘package:stats’:

    lowess
    

    Error: Use 'theme' instead. (Defunct; last used in version 0.9.1)

    Execution halted

    Any thoughts on how to fix this?

  • Geraldine_VdAuweraGeraldine_VdAuwera admin Cambridge, MAMember, Administrator, Broadie admin

    Ah, darn it. I've seen this before. This refers to a function that was merely deprecated in previous versions of R (or the library itself, not sure), but it seems you have a newer version that treats it as a blocking error. We need to update the script to use the newer function but this is not a priority so it might take a while longer to get done. In the meantime, you can either downgrade the software you have installed, or you can tweak the script to use the latest function names. Either way it's going to be a bit of hassle, sorry.

  • Geraldine_VdAuweraGeraldine_VdAuwera admin Cambridge, MAMember, Administrator, Broadie admin

    Actually, can you tell me the versions of R and of gplots that you have installed? That might help for debugging.

  • estif74estif74 Saint Paul, MN, USAMember

    Sure, that I think I can handle :-)

    gplots: 2.13.0
    R: 3.0.2

  • Geraldine_VdAuweraGeraldine_VdAuwera admin Cambridge, MAMember, Administrator, Broadie admin
  • estif74estif74 Saint Paul, MN, USAMember

    As an aside for the newbies out there, you can take this CSV file that is created here (see above) and make your own plots in Excel. It's clunky compared to the GATK solution, but you can at least see the results that way. Mine looked great!

  • MutagenicMutagenic Member

    I haven't tried this yet because I'm a newbie with R, but from what I can gather opts() needs to be replaced with theme() for all instances in the R script.

  • SheilaSheila admin Broad InstituteMember, Broadie, Moderator admin
  • estif74estif74 Saint Paul, MN, USAMember

    How exactly does one do this?

  • Geraldine_VdAuweraGeraldine_VdAuwera admin Cambridge, MAMember, Administrator, Broadie admin
  • estif74estif74 Saint Paul, MN, USAMember

    Thanks @Geraldine_VdAuwera‌ for posting this to dropbox. I'm a bit confused as to what to do with this exactly - I took a look at the script and it looks like your script is referencing 2 files that I don't have:

    1. NA12878.6.1.dedup.realign.recal.bqsr.grp.csv
    2. NA12878.6.1.dedup.realign.recal.bqsr.grp

    Do I need to edit the script to changes these values to my files and then re-run the Rscript? I do have a .csv file, but not a .grp file. Sorry to be so clueless!

  • estif74estif74 Saint Paul, MN, USAMember

    I ended up renaming the files RTest.csv, RTest.recal, and placed them in the same directory as your new BQSR.R script. Then I typed:

    RScript BQSR.R RTest.csv RTest.recal RTest.pdf

    at the command prompt. Worked like a charm. THANK YOU!

  • zhunterzhunter United StatesMember

    Thanks Geraldine! I recently went through and tried to debug the RScript files to get queue working on our system. Its all ggplot2 errors that can be easily fixed. You can find all the R scripts doing a search for files ending in .R and then scan them for anything using the depreciated "opts" command. Basically, change "opts" to "theme" and "theme_XXX()" -> "element_XXX()". If opts is being used to add a title (ie opts(title="my title")) then just rewrite it to ggtitle("my title").

    The exception is the ququeJobReport.R. In addition to the above, you also need to add "library(grid)" to the top of script with the other library calls as arrow() now requires grid. Hope this helps!

  • bioinfo_89bioinfo_89 IndiaMember

    Hi!!
    I get the same error which I got using BaseRecalibration step when I use the below command:

    Rscript BQSR.R report.csv recal_data.table recal.pdf

    Attaching package: ‘gplots’

    The following object is masked from ‘package:stats’:

    lowess
    

    Error: Use 'theme' instead. (Defunct; last used in version 0.9.1)
    Execution halted

    Please guide me through ...

  • bioinfo_89bioinfo_89 IndiaMember

    Sorry, I get this same error when I use AnalyzeCovariates.

  • Geraldine_VdAuweraGeraldine_VdAuwera admin Cambridge, MAMember, Administrator, Broadie admin
  • bioinfo_89bioinfo_89 IndiaMember

    Thnx Geraldine!! I did the changes @zhunter mentioned and got the following result. Can you tell me if it worked or not!

  • bioinfo_89bioinfo_89 IndiaMember

    I not able to exactly understand the plots!

  • bioinfo_89bioinfo_89 IndiaMember

    Also I would like to know how to plot Variant recalibration plots in the similar way. Is there any command like what was provided in the link using Rscript and BQSR.R script available?

    Thanks!

  • bioinfo_89bioinfo_89 IndiaMember

    I again ran the recalibration step and this time I got a better plot. Is there any documentation to understand what the plots are trying to say? First 9 graphs are easy to understand just by looking that recalibration worked but the following graphs show the reverse case which is confusing to me.

    Any help would be grateful!!

  • Geraldine_VdAuweraGeraldine_VdAuwera admin Cambridge, MAMember, Administrator, Broadie admin

    @bioinfo_89‌ Please try to group your questions or comments into fewer posts. You are posting too many messages in a short time. Think carefully about what you want to ask or say before posting, and you won't need to post three messages in 5 minutes.

    Yes, there is documentation that explains the plots. Look in the Guide section. I recommend you watch the video presentations in the "Presentations" section of the Guide.

    We just released a new version with many bugfixes; I believe if you use that the R scripts will work fine.

  • bioinfo_89bioinfo_89 IndiaMember

    I did that because I am working on that very issue currently. I will try to group my posts thank you!!

  • ravichavravichav United StatesMember

    I have the same problem and I am unable to access the dropbox R script. Is it still active ?

  • SheilaSheila admin Broad InstituteMember, Broadie, Moderator admin

    @ravichav‌

    Hi,

    I just checked this link https://www.dropbox.com/s/5bplh33rrjzr7v8/BQSR.R and it works.

    Regardless, if you upgrade to version 3.2, you will not need to change anything.

    -Sheila

  • ravichavravichav United StatesMember
  • ravichavravichav United StatesMember
    edited August 2014

    I tried the script. the .csv file generated and the original baserecal file and executed it as

    Rscript BQSR.R Rtest.csv Rtest.baserecal.table Rtest.pdf

    and got this warning

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

    Does it represent any issues?

  • luzluz Member

    @Geraldine_VdAuwera said:

    Yes, there is documentation that explains the plots. Look in the Guide section. I recommend you watch the video presentations in the "Presentations" section of the Guide.

    I don't understand the second part of the graph either. Indeed, why the reverse? I have searched the forum and gone through the all the presentation on the base recalibration but found no explanation. Thank you.

  • Geraldine_VdAuweraGeraldine_VdAuwera admin Cambridge, MAMember, Administrator, Broadie admin

    @luz, can you please quote which part you don't understand? Specifically, what do you mean by "why the reverse"?

  • Geraldine_VdAuweraGeraldine_VdAuwera admin Cambridge, MAMember, Administrator, Broadie admin

    Oh I see, @Sheila explained the problem to me. Can you please post all the command lines you used for the recalibration process?

  • luzluz Member

    Here are the commands I used:

    $ java -Xmx2g -jar ~/bin/software/GATK/GenomeAnalysisTK.jar -T BaseRecalibrator -R Ref3.fa -I Test_Mdup_srt_RG_realigned.bam -L chr21:36099955-42000057 -knownSites dbSNP_vcf_chr_21_6Mbp.vcf -o Test_Mdup_srt_RG_realigned_recal_data.table

    $ java -Xmx2g -jar ~/bin/software/GATK/GenomeAnalysisTK.jar -T BaseRecalibrator -R Ref3.fa -I Test_Mdup_srt_RG_realigned.bam -L chr21:36099955-42000057 -knownSites dbSNP_vcf_chr_21_6Mbp.vcf -BQSR Test_Mdup_srt_RG_realigned_recal_data.table -o Test_6Mbp_Mdup_srt_RG_realigned_post_recal_data.table

    $ java -Xmx2g -jar ~/bin/software/GATK/GenomeAnalysisTK.jar -T AnalyzeCovariates -R Ref3.fa -L chr21:36099955-42000057 -before Test_Mdup_srt_RG_realigned_recal_data.table -after Test_Mdup_srt_RG_realigned_post_recal_data.table --intermediateCsvFile Test_Mdup_srt_RG_realigned_recal_plots.csv

    $ Rscript ~/bin/software/GATK/BQSR.R Test_Mdup_srt_RG_realigned_recal_plots.csv Test_Mdup_srt_RG_realigned_recal_data.table Test_Mdup_srt_RG_realigned_recal_plots.pdf

    I have also attached the plots.

    Thank you.

  • Geraldine_VdAuweraGeraldine_VdAuwera admin Cambridge, MAMember, Administrator, Broadie admin

    @luz OK, so the good news is that you're doing things correctly. But there is a misunderstanding in the nature of the plots which I think causes your confusion.

    The first page of plots shows accuracy, i.e. the distance between what the quality scores are (observed) and what they should be (empirical). Before recalibration, the points are scattered, and after recalibration, they are more aligned (along the 0 mark), which is good.

    The second page shows the actual values (approximated by the mean) of the quality scores in your data, before and after recalibration. For that, there is no expectation of whether points will be more or less aligned after recalibration. In your case, it looks like the distribution of base qualities before and after are roughly the same, but the values are a little lower, which suggests that the machine was overly confident by about the same amount over all the data. But I guess what troubles you are the plots for indels, because the points are all neatly aligned before recalibration and scattered afterward. The reason is because before recalibration, there are no indel quality scores available, so the program gives them all the same default value. Then after recalibration, they have new values that are all different. This is a good sign that everything worked properly.

  • luzluz Member

    Now I understand. Thank you.

    When you said the machine was overly confident, should we be worrying about the performance of the Illumina Hiseq machine we are using? Any idea what people in the lab should be checking on?

    As for the indel "problem", I'm wondering if gaps in the mapping have any effect on this? I forgot to mention that this is sequence (not exome) capturing data with more than 12% gaps in the mapped region.

  • Geraldine_VdAuweraGeraldine_VdAuwera admin Cambridge, MAMember, Administrator, Broadie admin

    @luz,

    I can't really speak to the performance of your particular machine based on just this data. We expect all sequencers to have some issues estimating their own accuracy to some extent, in part due to external factors (handling of the machine, reagents) and in part internal factors (sequencing chemistry, manufacturing quality). In past testing (see figure attached) we have found that Illumina HiSeq machines tended to be overconfident across the board, so your result is not out of the ordinary, and is not necessarily cause for alarm. But it's good to keep an eye on this, e.g. by integrating the BQSR results into QC data tracking. Since QC is typically done on data before recalibration, it can paint an overly optimistic picture of the data quality. By adding a QC step after recalibration, you can get a more accurate evaluation, and if the empirical quality falls below expectations, then you can look into doing some additional performance of the hardware, reagents and personnel. We cannot provide hard numbers on what the quality of your data should be, but that is something that Illumina support should be able to discuss with you.

  • luzluz Member

    Hi Geraldine,
    Thanks for the reassurance. Good to know that.

  • ibseqibseq United KingdomMember

    @Geraldine_VdAuwera said:
    Have a look at this tutorial: http://www.broadinstitute.org/gatk/guide/article?id=4294

    Part of it explains how to run the plotting script manually. The hard-coded file names are only used in our internal testing / debugging; if you run it yourself you can specify your own filenames.

    Hi Geraldine,
    I had a similar problem and ran -l DEBUG but `I cannot see any R script lines. I attach my log output.
    Any advice on how to procede?

    thanks,
    ibseq

  • ibseqibseq United KingdomMember

    @luz said:
    Here are the commands I used:

    $ java -Xmx2g -jar ~/bin/software/GATK/GenomeAnalysisTK.jar -T BaseRecalibrator -R Ref3.fa -I Test_Mdup_srt_RG_realigned.bam -L chr21:36099955-42000057 -knownSites dbSNP_vcf_chr_21_6Mbp.vcf -o Test_Mdup_srt_RG_realigned_recal_data.table

    $ java -Xmx2g -jar ~/bin/software/GATK/GenomeAnalysisTK.jar -T BaseRecalibrator -R Ref3.fa -I Test_Mdup_srt_RG_realigned.bam -L chr21:36099955-42000057 -knownSites dbSNP_vcf_chr_21_6Mbp.vcf -BQSR Test_Mdup_srt_RG_realigned_recal_data.table -o Test_6Mbp_Mdup_srt_RG_realigned_post_recal_data.table

    $ java -Xmx2g -jar ~/bin/software/GATK/GenomeAnalysisTK.jar -T AnalyzeCovariates -R Ref3.fa -L chr21:36099955-42000057 -before Test_Mdup_srt_RG_realigned_recal_data.table -after Test_Mdup_srt_RG_realigned_post_recal_data.table --intermediateCsvFile Test_Mdup_srt_RG_realigned_recal_plots.csv

    $ Rscript ~/bin/software/GATK/BQSR.R Test_Mdup_srt_RG_realigned_recal_plots.csv Test_Mdup_srt_RG_realigned_recal_data.table Test_Mdup_srt_RG_realigned_recal_plots.pdf

    I have also attached the plots.

    Thank you.

    Hi ,
    I have ran similar commands althought the R script ended with an error message saying:

    Loading required package: methods
    Attaching package: ‘gplots’
    The following object is masked from ‘package:stats’:
    lowess
    Warning messages:
    1: NAs introduced by coercion
    2: NAs introduced by coercion
    NULL

    Quick note: I didn't have a dbSNPs to start with, so I have ran Haplotype caller and generated one. Then I used this vcf to re-calibrate and I was then doing the second pass recalibration where I encountered R-related issues...

    I attach my pdf. Does this look like ok?
    thanks for the help
    ibseq

  • SheilaSheila admin Broad InstituteMember, Broadie, Moderator admin

    @ibseq
    Hi ibseq,

    How big is your dataset? The error of NAs introduced by coercion may be due to a small dataset. We recommend running BQSR on at least 150M bases (1B is better).

    That said, your plots look fine, so I am not sure if the error messages are anything to be too concerned about. However, you should really consider upgrading to the latest version of GATK. The version you are using is very old.

    -Sheila

Sign In or Register to comment.