Questions about BQSR

Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie
edited November 2014 in Ask the GATK team
This discussion was created from comments split from: Base Quality Score Recalibration (BQSR).

To ask new questions, please post on the article linked above.

Comments

  • Hi,

    Could anyone provide a mathematical notes on what exactly BaseRecalibrator is doing?

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie
  • Hi,

    I've generated the plots during recalibration, but I'm only getting plots for the original data (not pairs of plots like above). How can I get the plots for the recalibrated scores in order to see what's going on with my data?

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    This has been addressed previously on this forum, please do a search (box in the top right corner) and you will find an explanation.

  • priesgopriesgo PanamáMember

    About the BaseRecalibrator performance. I am using GATKLite-2.3.4 and in a "best-practices"-like pipeline the base recalibration is the most time consuming task.

    The base recalibration (BaseRecalibrator + PrintReads) took about 5 hours for an exome sample, while local realignment took 1 hour end 15 minutes and variant discovery little more than 1 hour. I am using multi-threading whenever it is possible.

    I've read that there were some issues with BaseRecalibrator multi-threading and performance in previous releases, but when running BaseRecalibrator in a single thread performance does not seem to increase. Is there any alternative to increase this performance? Maybe disabling some functionality or preparing the data somehow?

    Thanks,

    Pablo.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Hi Pablo,

    Have you tried using scatter-gather (e.g. using Queue)? There was a recent forum discussion on the topic, you may find it interesting to look up.

  • Hello,
    Having trouble with the PrintReads aspect of the GATK 2+ recalibration scheme. I am using v2.3-9-ge5ebf34 with reducereads on illumina generated exomes. The BaseRecalibrator & PrintReads have been complaining about the cycle covariate length being too long

    ERROR MESSAGE: The maximum allowed value for the cycle is 500, but a larger cycle was detected in read C47. Please use the --maximum_cycle_value argument to increase this value (at the expense of requiring more memory to run)

    So I implemented a filtering step to clip the cycles > 500 as below. BaseRecalibrator now runs without errors, however PrintReads still complains maximum cycle length. Additionally it suggests I use the --maximum_cycle_value argument, which neither works in the command line nor appears in the documentation for PrintReads.

    Is this an error or omission in PrintReads? I may try HardClipping the bam file though this is not supported.

    thoughts?

    java -Xmx8g -jar /home/jpriest/priest_apps/GATKv2.3.9/GenomeAnalysisTK.jar \
    -T ClipReads \
    -I ./CHD001_novomask_sorted_realigned.rr.bam \
    -o ./CHD001_novomask.clip.sort.realign.rr.bam \
    -R /home/jpriest/reference_genomes/hg19/ucsc.hg19.fasta \
    -CT 500-10000 \
    -CR SOFTCLIP_BASES \
    -os testfilter.txt

    java -Xmx8g -jar /home/jpriest/priest_apps/GATKv2.3.9/GenomeAnalysisTK.jar \
    -T BaseRecalibrator \
    -I ./CHD001_novomask.clip.sort.realign.rr.bam \
    -R /home/jpriest/reference_genomes/hg19/ucsc.hg19.fasta \
    -knownSites /home/jpriest/reference_genomes/hg19/dbsnp_135.hg19.vcf \
    -knownSites /home/jpriest/reference_genomes/hg19/Mills_and_1000G_gold_standard.indels.hg19.vcf \
    -o CHD001.recal_data.grp \
    -nct 4 \
    -filterMBQ

    java -jar /home/jpriest/priest_apps/GATKv2.3.9/GenomeAnalysisTK.jar \
    -T PrintReads \
    -R /home/jpriest/reference_genomes/hg19/ucsc.hg19.fasta \
    -I ./CHD001_novomask.clip.sort.realign.rr.bam \
    -BQSR CHD001.recal_data.grp \
    -o ./CHD001_novomask.final.rr.bam \

  • update: using the -CR HARDCLIP_BASES \ flag with ClipReads seems to solve the problem, though this is a little unnerving as it is not directly supported

  • ebanksebanks Broad InstituteMember, Broadie, Dev

    Hi @jamesrpriest, I just re-read your original post and it dawned on me what you are doing: you are trying to run BQSR on a reduced BAM file. This is very, very bad - and does not at all follow our best practices recommendations for a good workflow. BQSR absolutely must be done before BAM reduction (which should be the very last step before variant calling) because reduction squashes all of the base qualities into a condensed form.

    That being said, it is true that reads that are too long will fail in the PrintReads step so I will add a patch for that case (to be available in the upcoming 2.4 release).

  • Can your input be multiple BAM files?

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    @hmk123, are you asking if BQSR accepts multiple bam files as input? If so the answer is yes.

  • Great thank you for the update

  • pmintpmint Member

    (# mismatches + 1) / (# bases + 2) --> why +1 & +2?

  • @pmint said:
    (# mismatches + 1) / (# bases + 2) --> why +1 & +2?

    We use this correction factor to smooth out low occupancy bins in the calculations. It is saying that each bin starts with one error observation and one true observation.

  • pmintpmint Member

    Oh! I understand. Thank you!!

  • I was wondering where I can find the QualityScoreCovariate.java source code. I was browsing through the source code on github but difficult to find as there are so many folders and classes, so would you be able to point me to the full path? Thanks!

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    It is located in package org.broadinstitute.sting.utils.recalibration.covariates

  • Is Base Quality Score Recalibration (BQSR) supported for Ion Torrent bams? It is stated that it is supported for: Illumina, SOLiD, 454, Complete Genomics, Pacific Biosciences, etc. and I was not sure if "etc." included Ion Torrent. The error I receive is "##### ERROR MESSAGE: The platform (PL) associated with read group [email protected] is not a recognized platform. Implemented options are e.g. illumina, 454, and solid" so obviously the algorithm is checking for a keyword for PL.

    Thanks.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Hi @Mutagenic, sorry about the vague error message. Ion data should work; we'll check and update the list of compatible platforms.

  • how many memorys need for GATK PrintReads ?

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    @murphy, it can depend on what you're doing with it, but you should start with 2 Gb and increase it if the program complains that it's not enough.

  • AshuAshu Member

    I want to set the quality scores in a SAM file to 20 as mentioned here -> http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3443046/. What program in GATK should I use to do this, I could not find that in the paper. Please get back as soon as you could. Thank you so much.

  • AshuAshu Member

    Thank you so very much.

  • AshuAshu Member

    The link you just gave me - reassigns Read Mapping Qualities. I want to re-assign the Base Qualities - (I beleive the BSQR recalibrates base qualities and not the read mapping qualities). Could you please let me know of a tool to reassign base quality?

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Oh, I didn't realize you wanted to change base qualities, sorry. I don't know of any tool that does that. Maybe in Samtools...

  • AshuAshu Member

    But, the paper I asked you about was published by Broad Institute, it says - "Quality scores for the bases were generated by using the GATK’s base quality score recalibration pipeline, starting with a default original base quality of Q20 for every base in the read. " This means that base qual scores in the BAM files were changed to a default value of 20 for every base. I just wanted to know of the tool you used to do that. Mauricio O Carneiro is the first author in the paper. - http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3443046/

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Hmm, I'm not familiar with the methodology used in that paper. I'll ask @Carneiro to clarify.

  • AshuAshu Member

    great thank you so much. i'd appreciate that.

  • CarneiroCarneiro Charlestown, MAMember

    That was before PacBio had quality scores in their BAM files. Now that they do, you can safely use their Base Qualities as priors for BQSR as you do with illumina reads.

  • AshuAshu Member

    Dear Carneiro,
    Thank you so much for replying.
    Ok, Just to confirm, you are saying that in case/If quality scores are absent in the pacbio's bam files, we need to update them to 20; but there was a time when pacbio's BAM files did not have quality scores?? and that is when you updated them to a default value of 20?? And I do not need to do this anymore. I just need this confirmation as I have to answer this to my leads :)
    Also, while waiting for your reply, I used the following command to change the qual scores to 20: Do you think I should get rid of the options ddq/idq/mdq??
    java -jar new_gatk/GenomeAnalysisTK.jar -T BaseRecalibrator -I test/A2.bam -R test/reference.fasta -knownSites test/knownSites.vcf -ddq 20 -idq 20 -mdq 20 --maximum_cycle_value 10000 -o test/A2_recal_data.grp

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Hi Ashu,

    That's right, what Mauricio is saying is that he had to artificially add quality scores to the PacBio bams because at the time PacBio didn't output quality scores. So if you have to work with older bam files, you need to add quality scores, but if you're working with more recent files that already have quality scores, you don't need to modify them.

  • AshuAshu Member

    awesome thanks a lot. :)
    I dont get any email updates when you reply to my comments. I dont know if that is an issue or is it supposed to work that way? Just letting you guys know, I have to search this page explicitly to read these comments regularly.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Well, it's not really a technical issue -- there is an option you can turn on to receive notifications, but the forum software doesn't let us turn it on by default, so people have to turn it on themselves. I forget exactly where (and am on a smartphone right now so I'm limited) but the setting is in your user profile somewhere. If you can't find it, remind me and I will look it up for you.

  • AshuAshu Member

    found it... thanks again :)

  • Dear all, did you find that the value in "Errors" column contained the decimal fraction, in the recalibration report generated by the the new version GATK (2.4-7), as shown in the picture . But, why? I think it should be a integer?
    Looking forward to your reply, many thanks!

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Hi there,

    That value is a count on the number of errors that are in that bin, but it is now a probabilistic count based on the likelihood of actually being at that location so it is a decimal value instead of an integer.

  • JackJack Member

    Hello there, I have been studying on a non-human organism (chicken). Recently, I have a problems confusing me for almost half a month: when I use the BaseRecalibrator to generate a grp file, it always complains there is something wrong with my input file.
    I downloaded the reference genome from UCSC and the snp files from NCBI, since all snp files are provided according to chromosomes, I merge them (linux split command) into a dbSNP.vcf file which I used as the input file. However, the GATK complained "The provided VCF file is malformed at approximately line number 626572: The VCF specification does not allow for whitespace in the INFO field." well, when I used a perl script to solve this problem and rerun, it again gives another complaints "Input vcf and reference have incompatible contigs: Relative ordering of overlapping contigs differs, which is unsafe" or "The provided VCF file is malformed at approximately line number 15: 4554 is not a valid start position in the VCF format" and so on.
    Now my work completely stop. I can not move on. Does anyone have encountered similar problems? any help is appreciated!

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Hi @Jack,

    I would recommend that you start again from scratch and try merging the original VCF files again, but this time instead of using linux commands, you should use a tool specialized for merging this type of file, such as the GATK's CombineVariants tool, or VCFtools.

  • JackJack Member

    Thanks for your quick response,Geraldine. Now I have fixed it by combing those vcf files using GATK CombineVariants, I can move on now!

  • We are now trying with the Queue.jar where we wrote a script to do scatter gather with the BQSR. However, whenever we run this script, we always end up not getting the pdf file and got the error: org.broadinstitute.sting.queue.QException: Unable to find gather inputs
    We have checked the director stated by the command and we couldn't find the intermediate file with the pdf. However, in the scala script, we did included the plot_pdf argument:

        trait GATKCommon extends CommandLineGATK{
                this.memoryLimit = pipeline.memory
                this.intervals = pipeline.intervals
                this.excludeIntervals = pipeline.excludeIntervals
                this.reference_sequence = pipeline.referenceFasta
    
            }
    def script(){    
        val baseRecalibrate = new BaseRecalibrator with GATKCommon
        baseRecalibrate.input_file:+= inputFile
        baseRecalibrate.scatterCount = pipeline.thread
        baseRecalibrate.knownSites = realignTarget.known
        baseRecalibrate.out = swapExt(indelRealigner.out, "bam", "bqsr.grp")
        baseRecalibrate.plot_pdf_file = swapExt(indelRealigner.out, "bam", "bqsr.pdf")
        add(baseRecalibrate)
    }       
    

    Is it that if we do scatter gather with the BQSR we will not get the pdf? If so, is there anyway we can quickly generate the pdf file? If not, is there any problem with our command?

    Thank you for your help

  • CarneiroCarneiro Charlestown, MAMember
    edited May 2013

    There's nothing obviously wrong about your code, and plots should run with scatter/gather (as far as I know -- this may have changed). But I have a few questions first:

    1. What version of the GATK are you using?
    2. What version of R are you
    3. Can you show us the output of any scattered BQSR job here to see if the plot was generated?

    Unfortunately there is no quick way to generate the plots other than through the program.

  • First, I use the latest GATK and R 2.15. When I tried to re-run the script with additional

    this.S_=(ValidationStringency.LENIENT)

    in the trait and were able to pull out more information. The run time error code was as follow:

    org.broadinstitute.sting.queue.QException: Unable to find gather inputs: /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_01_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_02_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_03_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_04_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_05_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_06_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_07_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_08_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_09_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf
            at org.broadinstitute.sting.queue.function.scattergather.GatherFunction$class.waitForGatherParts(GatherFunction.scala:70)
            at org.broadinstitute.sting.queue.function.scattergather.SimpleTextGatherFunction.waitForGatherParts(SimpleTextGatherFunction.scala:36)
            at org.broadinstitute.sting.queue.function.scattergather.SimpleTextGatherFunction.run(SimpleTextGatherFunction.scala:38)
            at org.broadinstitute.sting.queue.engine.InProcessRunner.start(InProcessRunner.scala:53)
            at org.broadinstitute.sting.queue.engine.FunctionEdge.start(FunctionEdge.scala:84)
            at org.broadinstitute.sting.queue.engine.QGraph.runJobs(QGraph.scala:434)
            at org.broadinstitute.sting.queue.engine.QGraph.run(QGraph.scala:156)
            at org.broadinstitute.sting.queue.QCommandLine.execute(QCommandLine.scala:171)
            at org.broadinstitute.sting.commandline.CommandLineProgram.start(CommandLineProgram.java:245)
            at org.broadinstitute.sting.commandline.CommandLineProgram.start(CommandLineProgram.java:152)
            at org.broadinstitute.sting.queue.QCommandLine$.main(QCommandLine.scala:62)
            at org.broadinstitute.sting.queue.QCommandLine.main(QCommandLine.scala)
    

    Typing ls /.pdf within the ./queue/scatterGather/Pipeline-5-sg folder I found only one of the pdf file (temp_10_of_10)

    The .out file within those 9 folders were relatively long but I noticed the following error:

    WARN 11:35:23,960 RScriptExecutor - RScript exited with 1. Run with -l DEBUG for more info.

    As a result of that, I also enable DEBUG and look into the RScript errors:

    Loading required package: methods
    Loading required package: gtools
    Loading required package: gdata
    gdata: read.xls support for 'XLS' (Excel 97-2004) files ENABLED.
    
    gdata: read.xls support for 'XLSX' (Excel 2007+) files ENABLED.
    
    Attaching package: âgdataâ
    
    The following object(s) are masked from âpackage:statsâ:
    
        nobs
    
    The following object(s) are masked from âpackage:utilsâ:
    
        object.size
    
    Loading required package: caTools
    Loading required package: grid
    Loading required package: KernSmooth
    KernSmooth 2.23 loaded
    Copyright M. P. Wand 1997-2009
    Loading required package: MASS
    
    Attaching package: âgplotsâ
    
    The following object(s) are masked from âpackage:statsâ:
    
        lowess
    
    Loading required package: plyr
    
    Attaching package: âreshapeâ
    
    The following object(s) are masked from âpackage:plyrâ:
    
        rename, round_any
    
    Error in `colnames<-`(`*tmp*`, value = c("ReadGroup", "EventType", "EmpiricalQuality",  :
      'names' attribute [6] must be the same length as the vector [1]
    Calls: source ... finishTable -> .gsa.assignGATKTableToEnvironment -> colnames<-
    Execution halted
    

    It seems like there is some error with the RScript? I am not sure why such error occurs considering that I am using the Queue directly without manipulating the data....

  • CarneiroCarneiro Charlestown, MAMember
    edited May 2013

    the error is before the Rscript.

    Unable to find gather inputs

    something happened in your Queue run that your scattered jobs did not produce the right outputs or your filesystem may have failed there? It seems like a network failure that prevented the scattered files to be reached by the gather function. Without all the files, the Rscript will fail.

    Can you check what happened in the .queue scatter directories, see what output files are there? Maybe you can figure out what happened on your side. If not, just rerun with -startFromScratch and hopefully this won't happen again.

  • Thank you Carneiro, I will try and see if that is the case. I am not sure why only one folder contain the pdf file though. Maybe startFromScratch will help.

  • I have just finished the startFromScratch yet I still got the following error:

    SimpleTextGatherFunction: List(/nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_01_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_02_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_03_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_04_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_05_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_06_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_07_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_08_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_09_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_10_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf) > List(/nfs/home/sam/Test_Zone/ScalaTest/105NT.chrY.dedup.indelRe.bqsr.pdf)
    org.broadinstitute.sting.queue.QException: Unable to find gather inputs: /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_01_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_02_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_03_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_04_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_05_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_06_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_07_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_08_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf, /nfs/home/sam/Test_Zone/ScalaTest/.queue/scatterGather/Pipeline-5-sg/temp_09_of_10/105NT.chrY.dedup.indelRe.bqsr.pdf
        at org.broadinstitute.sting.queue.function.scattergather.GatherFunction$class.waitForGatherParts(GatherFunction.scala:70)
        at org.broadinstitute.sting.queue.function.scattergather.SimpleTextGatherFunction.waitForGatherParts(SimpleTextGatherFunction.scala:36)
        at org.broadinstitute.sting.queue.function.scattergather.SimpleTextGatherFunction.run(SimpleTextGatherFunction.scala:38)
        at org.broadinstitute.sting.queue.engine.InProcessRunner.start(InProcessRunner.scala:53)
        at org.broadinstitute.sting.queue.engine.FunctionEdge.start(FunctionEdge.scala:84)
        at org.broadinstitute.sting.queue.engine.QGraph.runJobs(QGraph.scala:434)
        at org.broadinstitute.sting.queue.engine.QGraph.run(QGraph.scala:156)
        at org.broadinstitute.sting.queue.QCommandLine.execute(QCommandLine.scala:171)
        at org.broadinstitute.sting.commandline.CommandLineProgram.start(CommandLineProgram.java:245)
        at org.broadinstitute.sting.commandline.CommandLineProgram.start(CommandLineProgram.java:152)
        at org.broadinstitute.sting.queue.QCommandLine$.main(QCommandLine.scala:62)
        at org.broadinstitute.sting.queue.QCommandLine.main(QCommandLine.scala
    

    Again, I only got a pdf on temp_10_of_10, however, I have checked that PDF and it is empty. To make sure it is a real problem, I have re-run the programme 3 times and I still got the same error message. So I am not sure if that is the

  • CarneiroCarneiro Charlestown, MAMember

    this looks like a bug, let me take a look.

  • Thank you Carneiro

  • CarneiroCarneiro Charlestown, MAMember

    it is a bug, we are fixing it internally and will be in the nightly build in the next few days. In the meantime, I recommend you don't use the -plots option when scatter/gathering.

  • Dear Carneiro,

    Thank you very much, it is nice to know that it is being solved. Thank you for your help!

    Sam

  • does GATK support Ion_Torrent data for teh base quality recalibration process?
    (i know taht it was asked here before, but i wonder if anything changed since then.)
    i get this message error:
    RROR MESSAGE: The platform (ION_TORRENT) associated with read group [email protected] is not a recognized platform. Implemented options are e.g. illumina, 454, and solid

    Thanks :-)

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Hey @galmoran, can you tell me what version of GATK you are using? That looks like an older version of the platform lineup. We added Ion Torrent to the allowable platforms a little while ago IIRC.

  • Hello, I was wondering if BQSR was using depth as an indication of erroneous base. Does it look at every single reads or only at a consensus sequence with degenerate base?
    Thanks

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Hi Yannick,

    BQSR looks at all reads, but it does so independently for each read, using its CIGAR information (= it doesn't look at pileups). So depth is not among the covariates used for recalibration.

  • JackJack Member

    Hi,there. I'm a little confused about BQSR, does this step rely on a large dataset of dbSNP? I mean my dbSNP is not very large, does the results will still be better if I run with my small dbSNP set than that got without recalibrating the realn.bam file ?

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Yes, it is better to do BQSR with a small database rather than not do it at all.

  • Hi there. I'm wondering if you can help with an issue I've been stumbling over with BQSR? I'm working with a non-model species; we have done exon-capture and sequencing on about 600 individuals. As part of our SNP calling pipeline, we are planning to use BQSR, but we need to build a new dbSNP. To do this, we used BWA-mem to align, Picard to mark duplicates, and samtools mpileup to call SNPs. The mpileup output vcf files are very large, and I was wondering how I can reformat the vcf file to omit the genotypes and just have the first 9 columns that are necessary to use a vcf file as dbSNP in BQSR? I've done some testing but I keep getting error messages like this one:

    The provided VCF file is malformed at approximately line number 1036: there are 1 genotypes while the header requires that 9 genotypes be present for all records at Scaffold_1:6720

    I feel I must be doing something obviously wrong, but I can't seem to find any documentation telling me how to fix it. Thanks for all your hard work on this useful tool!

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Hi @yeaman,

    Simplest way to do that is to feed your variants to CombineVariants with the -minimalVCF argument. It will produce a sites-only file without any genotypes.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie
    edited September 2013

    Oh wait, just realized -minimalVCF is the wrong argument; it will strip out some info but not the genotypes. What you want is --sites_only. It's not currently documented anywhere, but you can use it with any GATK tool that writes a VCF.

  • Thanks for the quick reply! I tried running CombineVariants with -sitesOnly but got an error message saying:

    "Argument with name 'sitesOnly' isn't defined"

    I tried upgrading to the newest version but still got the same message. I also tried -minimalVCF but as you mentioned, that won't do it. I forgot to mention above that I had previously just used AWK to strip the genotypes off, leaving the first 9 columns. Any thoughts?

  • Just for reference, this is the command I was using:

    java -Xmx32G -XX:MaxPermSize=16g -XX:PermSize=16g -jar /data/programs/GenomeAnalysisTK-2.7-2/GenomeAnalysisTK.jar -T CombineVariants -R /data/seqcap/spruce/pseudo.sequences.all3_reduced_genome.fa --sitesOnly --variant var.flt.vcf.1 --variant var.flt.vcf.2 -o test_combine_variants.vcf --assumeIdenticalSamples -genotypeMergeOptions UNSORTED

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Yeah sorry that was my fault, I misremembered the name (that's the problem with undocumented arguments), the actual argument name is --sites_only. This is useful if you want to create a set of knowns with one of our callers and make GATK output a sites-only file directly, and not have to strip it in a second step.

    If you awk the samples info out, remember to edit the header accordingly, otherwise the parser will expect additional fields and flip out (which might be what you're seeing already).

  • Thanks! That worked like a charm and seems to be running fine with the BQSR tools. I fiddled around with awk and the headers for a bit but had trouble getting a format that would work (an assortment of other errors). But that's unimportant now that CombineVariants is working.

  • For future reference, I found a quicker way to use awk to chop off the genotypes and create a vcf formatted file that BQSR will accept. I printed out the first 8 columns (#CHROM POS ID REF ALT QUAL FILTER INFO), excluding the 'format' column (rather than the 6 or 9 I had tried previously) and removing the '##FORMAT' lines from the header. Using grep and cat, I then combined these together and added a header back on. If you use GNU parallel with multiple threads for the awk chopping, this can be quite fast. It doesn't seem to matter that the order of the contents of the 'INFO' column differs between the mpileup output and the format that is outputted by 'CombineVariants'.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Good to know, thanks for reporting your solution, @yeaman.

  • ajc8ajc8 Member
    edited October 2013

    Hello,
    I am running BQSR and I'm getting the following error message:

    ERROR MESSAGE: SAM/BAM file SAMFileReader{/Volumes/Passport/IA_084/IA084C.realigned.fixed.dedup.bam} is malformed: BAM file has a read with mismatching number of bases and base qualities. Offender: HWI-ST1122:264:C2LCWACXX:1:2113:17084:65824 [101 bases] [0 quals]

    So I found the place in the bam using grep and it looks like this:

    HWI-ST1122:264:C2LCWACXX:1:2113:17084:65824 99  1   13405   22  101M    =   13495   186 CCTCCACCACCCCGAGATCACATTTCTCACTGCCTTTTGTCTGCCCAGTTTCACCAGAAGTAGGCCTCTTCCTGACAGGCAGCTGCACCACTGCCTGCCGC   *   PG:Z:MarkDuplicates RG:Z:IA084C NM:i:1  MQ:i:22 AS:i:97 XS:i:97
    HWI-ST1122:264:C2LCWACXX:1:2113:17084:65824 147 1   13495   22  5S96M   =   13405   -186    CTCCTCTGCCTGGCGATGTGCCCGTCCTTTGCTCTGACCGCTGGAGACAGTGTTTGTCATTGGCATGGTCTGCAGGGATCCTGCTACAAAGGTGAAACCCA   " 
    
    
                                                              "
    
                                                               "!
    
    
    
     !
        PG:Z:MarkDuplicates RG:Z:IA084C NM:i:6  MQ:i:22 AS:i:66 XS:i:66
    

    I am guessing that the large space with the ", "! and ! is not supposed to be there. How can I remove this from the bam file? This is not the only exome I'm having trouble with at this step - we just got 8 exome sequences back from a new company and half of them are giving me strange errors like this. Thanks for your help.

    Post edited by Geraldine_VdAuwera on
  • Also, I tried adding "-filterMBQ" to the command and got a different error:

    ERROR stack trace

    java.lang.ArrayIndexOutOfBoundsException: -5
    at org.broadinstitute.sting.utils.baq.BAQ.calcEpsilon(BAQ.java:158)
    at org.broadinstitute.sting.utils.baq.BAQ.hmm_glocal(BAQ.java:246)
    at org.broadinstitute.sting.utils.baq.BAQ.calcBAQFromHMM(BAQ.java:542)
    at org.broadinstitute.sting.utils.baq.BAQ.calcBAQFromHMM(BAQ.java:595)
    at org.broadinstitute.sting.utils.baq.BAQ.calcBAQFromHMM(BAQ.java:530)
    at org.broadinstitute.sting.utils.baq.BAQ.baqRead(BAQ.java:663)
    at org.broadinstitute.sting.gatk.walkers.bqsr.BaseRecalibrator.calculateBAQArray(BaseRecalibrator.java:428)
    at org.broadinstitute.sting.gatk.walkers.bqsr.BaseRecalibrator.map(BaseRecalibrator.java:243)
    at org.broadinstitute.sting.gatk.walkers.bqsr.BaseRecalibrator.map(BaseRecalibrator.java:112)
    at org.broadinstitute.sting.gatk.traversals.TraverseReadsNano$TraverseReadsMap.apply(TraverseReadsNano.java:203)
    at org.broadinstitute.sting.gatk.traversals.TraverseReadsNano$TraverseReadsMap.apply(TraverseReadsNano.java:191)
    at org.broadinstitute.sting.utils.nanoScheduler.NanoScheduler$MapReduceJob.run(NanoScheduler.java:468)
    at java.util.concurrent.Executors$RunnableAdapter.call(Executors.java:439)
    at java.util.concurrent.FutureTask$Sync.innerRun(FutureTask.java:303)
    at java.util.concurrent.FutureTask.run(FutureTask.java:138)
    at java.util.concurrent.ThreadPoolExecutor$Worker.runTask(ThreadPoolExecutor.java:895)
    at java.util.concurrent.ThreadPoolExecutor$Worker.run(ThreadPoolExecutor.java:918)
    at java.lang.Thread.run(Thread.java:680)

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    If you got the exome bams like that straight from the company, you should contact their customer support and let them know you have this problem. They should not be giving you malformed bams. If you did some kind of processing before getting this error, you should check if the originals were okay.

  • For large datasets, it is often reasonable to split the base quality recalibration process across chromosomes. Given that I have one bam file and one grp file per chromosome, is it possible to merge the reads with the printReads command, e.g:

    java -jar GenomeAnalysisTK.jar \
    -T PrintReads \
    -I input_1.bam -I input_2.bam \
    -BQSR rep_1.grp -BQSR rep_1.grp \
    -o output.bam

  • Hi,
    Thanks for your help. I only got the raw fastq files from the company - could these errors be something in the raw data?

  • I have been using GATK 2.3. So I downloaded the latest version (2.7) and updated java. I'm still getting the error message: Exception in thread "main" java.lang.UnsupportedClassVersionError: org/broadinstitute/sting/gatk/CommandLineGATK : Unsupported major.minor version 51.0

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Hi @julian,

    No, you have to run PrintReads individually per chromosome bam because there is no way to pass multiple grp files. But the good news is that in the next step, if you're doing ReduceReads, you can just pass in all the chromosome as a list, and GATK will produce a single merged and reduced output bam.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    @acj8, check the FASTQ records and see if the same reads are okay or not in there. Maybe it's the aligner that screwed them up -- what program (and what version) did you use?

    Re: the java issue, if you're on a Mac, you have to install the java 7 JDK from Oracle, not just the regular JRE.

  • ajc8ajc8 Member
    edited October 2013

    Hi again,
    I had not installed the correct Java, but doing so did not make a difference. Here are the two fastq reads:

    @HWI-ST1122:264:C2LCWACXX:1:2113:17084:65824 1:N:0:AAACAT
    CCTCCACCACCCCGAGATCACATTTCTCACTGCCTTTTGTCTGCCCAGTTTCACCAGAAGTAGGCCTCTTCCTGACAGGCAGCTGCACCACTGCCTGCCGC
    +
    [email protected][email protected]@GGE><?9EHBH<F<[email protected]:?D<[email protected];@C)[email protected]?B7?;;A?;=??##################
    
    @HWI-ST1122:264:C2LCWACXX:1:2113:17084:65824 2:N:0:AAACAT
    TGGGTTTCACCTTTGTAGCAGGATCCCTGCAGACCATGCCAATGACAAACACTGTCTCCAGCGGTCAGAGCAAAGGACGGGCACATCGCCAGGCAGAGGAG
    +
    :[email protected]+22?22<+++22,@A7++3A7=<=+1**1?A#################################################################
    

    And here is the original bam:

    HWI-ST1122:264:C2LCWACXX:1:2113:17084:65824 99  1   13405   22  101M    =   13495   186 CCTCCACCACCCCGAGATCACATTTCTCACTGCCTTTTGTCTGCCCAGTTTCACCAGAAGTAGGCCTCTTCCTGACAGGCAGCTGCACCACTGCCTGCCGC   [email protected][email protected]@GGE><?9EHBH<F<[email protected]:?D<[email protected];@C)[email protected]?B7?;;A?;=??##################   NM:i:1  AS:i:97 XS:i:97 RG:Z:IA084C
    HWI-ST1122:264:C2LCWACXX:1:2113:17084:65824 147 1   13495   22  5S96M   =   13405   -186    CTCCTCTGCCTGGCGATGTGCCCGTCCTTTGCTCTGACCGCTGGAGACAGTGTTTGTCATTGGCATGGTCTGCAGGGATCCTGCTACAAAGGTGAAACCCA   #################################################################A?1**1+=<[email protected],22+++<[email protected]+==:   NM:i:6  AS:i:66 XS:i:66 RG:Z:IA084C
    

    I am using bwa-mem to do the alignment and then making the bam in samtools. Does any of this look strange to you.
    I'm currently pulling the ID out of the bam files created at each step (e.g. indel realignment (GATK), mark duplicates (Picard).

    Allison

    Post edited by Geraldine_VdAuwera on
  • ajc8ajc8 Member
    edited October 2013

    update: it is happening during the indel realignment - this is the entry from the realigned bam:

    MDF869:samtools-0.1.18 Allison$ ./samtools view /Volumes/Passport/IA_084/IA084C.realigned.bam |grep HWI-ST1122:264:C2LCWACXX:1:2113:17084:65824
    HWI-ST1122:264:C2LCWACXX:1:2113:17084:65824 99  1   13405   22  101M    =   13495   186 CCTCCACCACCCCGAGATCACATTTCTCACTGCCTTTTGTCTGCCCAGTTTCACCAGAAGTAGGCCTCTTCCTGACAGGCAGCTGCACCACTGCCTGCCGC   *   RG:Z:IA084C NM:i:1  MQ:i:22 AS:i:97 XS:i:97
    HWI-ST1122:264:C2LCWACXX:1:2113:17084:65824 147 1   13495   22  5S96M   =   13405   -186    CTCCTCTGCCTGGCGATGTGCCCGTCCTTTGCTCTGACCGCTGGAGACAGTGTTTGTCATTGGCATGGTCTGCAGGGATCCTGCTACAAAGGTGAAACCCA   " 
    
    
                                                                                              "
    
                                                                                               "!
    
    
    
     !
        RG:Z:IA084C NM:i:6  MQ:i:22 AS:i:66 XS:i:66
    
    Post edited by Geraldine_VdAuwera on
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Hi @ajc8, that looks like the realigner is messing up the read. Could you make a couple of BAM snippets (before/after realignment) with the read (plus about 100 bp of padding before/after it), and upload it to our FTP so we can debug locally? Instructions are here: http://www.broadinstitute.org/gatk/guide/article?id=1894

  • Yes, I will do that now. I have been troubleshooting with one of the files since my last post. I normally fix mates and mark duplicates in Picard after realigning around indels, but I switched the order - first Picard steps and then indel realignment. After doing so, I was able to use -filterMBQ (was getting an error before) - which removed 5% of the reads. Then, running BQST still gives me an error but not a malformed bam error, but instead saying that there is an array out of index error.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Hmm. Well, we typically recommend doing MarkDuplicates then Indel Realignment, so you're good there. Note that Indel Realigner should fix mates for you automatically, so you don't need to run FixMate separately. For the BQSR error, what's your command line and exact error message?

  • hi,
    Well I ran everything through again and I got the malformed bam error message again.

    Here is the command line:java -Xmx8g -jar GenomeAnalysisTK.jar -T BaseRecalibrator -nct 4 -I /Volumes/Thunderbolt/IA_084/IA084C.realigned.bam -R /Volumes/Thunderbolt/ref_genome/human_g1k_v37.fasta -knownSites /Volumes/Thunderbolt/2.3/b37/1000G_phase1.indels.b37.vcf -knownSites /Volumes/Thunderbolt/2.3/b37/Mills_and_1000G_gold_standard.indels.b37.vcf -knownSites /Volumes/Thunderbolt/2.3/b37/dbsnp_137.b37.vcf -o /Volumes/Thunderbolt/IA_084/IA084C.recal_data.grp

    And the error message is:

    SAM/BAM file SAMFileReader{/Volumes/Thunderbolt/IA_084/IA084C.realigned.bam} is malformed: BAM file has a read with mismatching number of bases and base qualities. Offender: HWI-ST1122:264:C2LCWACXX:1:2113:17084:65824 [101 bases] [0 quals]

    If I add -filterMBQ to the command, I get this error message:

    ERROR stack trace

    java.lang.ArrayIndexOutOfBoundsException: -5

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Oh well that's still the same problem -- that bad read is causing different issues depending on which tool you run, because they use different parts of the read info. You have to find out exactly what is causing this bad read to occur before you can continue with your work. That's why I asked for the before/after snippets.

  • I made the snippet from the pre-aligned bam no problem. When I try to do it with the aligned bam, I get the same malformed read error message. Is it okay if I just send you the prealigned snippet? If not, what should I do?

    Here is the command and error message:

    java -Xmx8g -jar GenomeAnalysisTK.jar -T PrintReads -R /Volumes/Thunderbolt/ref_genome/human_g1k_v37.fasta -I /Volumes/Thunderbolt/IA_084/IA084C.realigned.bam -o /Volumes/Thunderbolt/IA_084/IA084C.aligned_snippet.bam -L 1:13300-13600

    ERROR MESSAGE: SAM/BAM file SAMFileReader{/Volumes/Thunderbolt/IA_084/IA084C.realigned.bam} is malformed: BAM file has a read with mismatching number of bases and base qualities. Offender: HWI-ST1122:264:C2LCWACXX:1:2113:17084:65824 [101 bases] [0 quals]

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Sure, that's fine. Let me know when you have uploaded it.

  • I just uploaded it. The filename is ajc8_GATK_help. I also unintentionally uploaded it to a file on the site - Breakendbug. Can you delete it for me? Sorry about that.
    Thanks again for all your help.

  • Is there any difference between
    1: Running BQSR on all samples simultaneously generating a single .grp file, then running printreads for each sample, generating multiple recalibrated bams.
    2: Running BQSR on all samples simultaneously generating a single.grp file, then running printreads for all samples simultenaously, generating a single recalibrated bam.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Hi @chongm,

    No, there is no functional difference in terms of how the data will be processed, since recalibration is done at the read group level. However, be aware that there is some computational cost to running all the samples through simultaneously. We typically process lane-level or sample-level bams individually for performance reasons. Also, keep in mind that if you want to compress your files using ReduceReads, you need to process them individually (you cannot do multisample compression).

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Hi @ajc8, sorry to get back to you so late -- your bug slipped off my radar somehow. So, I've looked at your files and ran through everything, but I can't reproduce your bug -- at least not with the current version of GATK, which is 2.7. I realize reading your logs that you were using 2.3, so it's very possible that the bug was fixed by subsequent updates. But of course if you run into this again with the latest GATK, please do let me know. Good luck with your work!

  • My apologies for posting an error in this forum, but given the fact that "this (error) should never happen," I decided I needed to post it. If Mauricio, or any other members of the Broad/GATK team, can provide insight as to why one of my BAMs consistently gets hung up with the BQSR step of the GATK "Best Practices" workflow, I would be greatly appreciative. I am processing additional BAMs using the same pipeline (all other which past the BQSR step fine) and have previously used GATK for local realignment of my this problematic BAM file. Is it possible this bug could have been fixed if I run a later version of GATK? If I updated my version of GATK to a more recent version, will I need to re-run all my samples using the same version? Look forward to hearing back from you

    Thanks!

    ERROR ------------------------------------------------------------------------------------------
    ERROR A GATK RUNTIME ERROR has occurred (version 2.3-9-ge5ebf34):
    ERROR
    ERROR Please visit the wiki to see if this is a known problem
    ERROR If not, please post the error, 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: START (96) > (95) STOP -- this should never happen -- call Mauricio!
    ERROR ------------------------------------------------------------------------------------------
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Hi @mpvivero,

    I'm happy to inform you (on Mauricio's behalf) that this is a known bug that has been fixed in subsequent versions of GATK.

    Our recommendation is to use the same version of GATK for all processing steps, because that is the safest way to avoid version-change incompatibilities. So ideally you should reprocess all your samples with the newer version when you upgrade, yes. Sorry for the inconvenience!

  • I will try using the latest version then! Thank you!

  • Hello all,

    Can someone explain exactly how covariation works, in simple terms? I looked at the presentations in the dropbox folder from the link at the beginning of the comments, but I still don't understand. How exactly are context, cycle, read quality, and read group used to recalibrate the qualities? What is the dependent variable? And how does one determine that the rate of mismatch is higher than expected? If we determine that a covariate (context, e.g.) has a high positive covariation with the dependent variable (which is quality?), then how does that help us recalibrate the quality? Any help is highly appreciated! Thanks!

    • Nik.
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Hi Nik,

    Have a look at the video of the presentation on BQSR from out latest workshop here, hopefully it will help:

    http://www.broadinstitute.org/gatk/guide/events?id=3391

  • Hello,

    I am trying to analyse a rat genome (I am still a newbie) and I have downloaded the Ensembl reference from igenomes and more recently the .vcf file from Ensembl. The latest version of the Ensembl available ".vcf" file is 4.2 so it seems it is not supported. but when I try to run the previous version (4.1) (I got rid as well of some non standard annotations used in the .vcf file) but I am getting the following error:

    _##### ERROR MESSAGE: Input files /Volumes/Data2/genomes/Rn5/Rattus_norvegicus/Ensembl/Rnor_5.0/Annotation/Rattus_norvegicus_mod2.vcf and reference have incompatible contigs: Relative ordering of overlapping contigs differs, which is unsafe.

    ERROR /Volumes/Data2/genomes/Rn5/Rattus_norvegicus/Ensembl/Rnor_5.0/Annotation/Rattus_norvegicus_mod2.vcf contigs = [1, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 2, 20, 3, 4, 5, 6, 7, 8, 9, MT, X]
    ERROR reference contigs = [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 1, 20, 2, 3, 4, 5, 6, 7, 8, 9, MT, X]_

    Any ideas or suggestions about why might this be happening or how to solve it?
    Thanks in advance

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Hi @ghlopez,

    As the error message states, the order of contigs is not the same in your data and in your reference. You need to re-order the data to match the reference. Have a look at this document for more details: http://www.broadinstitute.org/gatk/guide/article?id=1204

  • Hello,

    I'm working on an exome capture project in a non-model species where we have multiplexed 12 individuals/lane of Illumina HiSeq. We have been running BQSR to give a recalibration table for each individual, but we thought perhaps it would be better to pool across all individuals in a lane before running BQSR, so that we'd have more data for the error profile. Do you have any recommendations about the best practices here?

    Thanks!
    Sam

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Hi @yeaman,

    I think in production we separate multiplexed individuals from each lane as well, but in theory there is nothing that prevents you from doing recalibration on them together in order to have more data for the model. It depends how much data you're getting per individual per lane. May be worth comparing the recalibration results you get in both cases.

  • Thanks for the help! If we get around to trying the pooling, I'll report back.

  • mlyonmlyon UKMember

    Hi Geraldine,

    Could you explain how the cycle number is referenced within the BQSR GATK report?

    Best wishes,
    Matthew

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie
    edited May 2014

    Hi @mlyon‌, if I recall correctly it should be referenced as a number in the CovariateValue column for lines where CovariateName is Cycle.

  • vifehevifehe SpainMember

    Hi,

    I have relaized that my recal_realigned.bam (280 M) file is much smaller than the previous realigned.bam (6.5 G) file. I don't think that is normal, right? However I am not getting any error message.
    Any hints on where may be any possible error?

    Thanks

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    @vihefe That doesn't sound right. Are you sure the job completed successfully? If so, what was your command line?

This discussion has been closed.