We've moved!
This site is now read-only. You can find our new documentation site and support forum for posting questions here.
Be sure to read our welcome blog!

ReadAdaptorTrimmer failed because of the unexcepted read order

Question 1:
ReadAdaptorTrimmer failed to process sorted bam file.

Version:GATK 3.30 and GATK 2.8.1

In the map() process
public List map( final ReferenceContext ref, final GATKSAMRecord readIn, final RefMetaDataTracker metaDataTracker )

Current process except that the read order will be _1 read , then _2 read, then another _1 read.

but In fact, the read oder is NOT as expected, so the process will failed.

Question 2:
We need keep the getAlignmentStart() info in cache.
Without this info, the SAM write will fail.

     firstReadInPair =GATKSAMRecord.emptyRead(readIn);
        firstReadInPair.setReadString(readIn.getReadString());
        firstReadInPair.setReadName(readIn.getReadName());
        firstReadInPair.setBaseQualities(readIn.getBaseQualities());

    //bugfix
    firstReadInPair.setAlignmentStart(readIn.getAlignmentStart());

Best Regards
Wang Yugui

Answers

  • SheilaSheila Broad InstituteMember, Broadie ✭✭✭✭✭

    @wang_yugui
    Hi Wang,

    Can you clarify what you are trying to ask? Did you run ReadAdapterTrimmer and get an error? Please post the commands and the errors you got.

    -Sheila

  • wang_yuguiwang_yugui Member

    Hi,Sheila

    ReadAdapterTrimmer failed with the following error message.

    java -XX:-UseCompressedOops -Xms100g -jar /usr/hpc-bio/gatk/GenomeAnalysisTK.jar -T ReadAdaptorTrimmer --minMatches 41 -R /usr/bio-ref/GRCh38.79/sequences.full.fa --bam_compression 0 --removeUnpairedReads -I md.bam -o trim.bam
    INFO 23:08:39,819 HelpFormatter - --------------------------------------------------------------------------------
    INFO 23:08:39,821 HelpFormatter - The Genome Analysis Toolkit (GATK) v3.3-0-g37228af, Compiled 2014/10/24 01:07:22
    INFO 23:08:39,821 HelpFormatter - Copyright (c) 2010 The Broad Institute
    INFO 23:08:39,821 HelpFormatter - For support and documentation go to http://www.broadinstitute.org/gatk
    INFO 23:08:39,825 HelpFormatter - Program Args: -T ReadAdaptorTrimmer --minMatches 41 -R /usr/bio-ref/GRCh38.79/sequences.full.fa --bam_compression 0 --removeUnpairedReads -I /biowrk/bam.md/Project_14686/Sample_100T/md.bam -o trim.bam
    INFO 23:08:39,830 HelpFormatter - Executing as [email protected] on Linux 3.10.0-229.1.2.el7.x86_64 amd64; OpenJDK 64-Bit Server VM 1.7.0_79-mockbuild_2015_04_15_00_14-b00.
    INFO 23:08:39,830 HelpFormatter - Date/Time: 2015/04/28 23:08:39
    INFO 23:08:39,831 HelpFormatter - --------------------------------------------------------------------------------
    INFO 23:08:39,831 HelpFormatter - --------------------------------------------------------------------------------
    INFO 23:08:39,921 GenomeAnalysisEngine - Strictness is SILENT
    INFO 23:08:40,031 GenomeAnalysisEngine - Downsampling Settings: No downsampling
    INFO 23:08:40,039 SAMDataSource$SAMReaders - Initializing SAMRecords in serial
    INFO 23:08:40,063 SAMDataSource$SAMReaders - Done initializing BAM readers: total time 0.02
    INFO 23:08:40,145 GenomeAnalysisEngine - Preparing for traversal over 1 BAM files
    INFO 23:08:40,149 GenomeAnalysisEngine - Done preparing for traversal
    INFO 23:08:40,150 ProgressMeter - [INITIALIZATION COMPLETE; STARTING PROCESSING]
    INFO 23:08:40,151 ProgressMeter - | processed | time | per 1M | | total | remaining
    INFO 23:08:40,151 ProgressMeter - Location | reads | elapsed | reads | completed | runtime | runtime
    INFO 23:08:40,153 ReadShardBalancer$1 - Loading BAM index data
    INFO 23:08:40,154 ReadShardBalancer$1 - Done loading BAM index data
    INFO 23:08:43,948 GATKRunReport - Uploaded run statistics report to AWS S3

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

    java.lang.NullPointerException
    at org.broadinstitute.gatk.tools.walkers.readutils.ReadAdaptorTrimmer.map(ReadAdaptorTrimmer.java:176)
    at org.broadinstitute.gatk.tools.walkers.readutils.ReadAdaptorTrimmer.map(ReadAdaptorTrimmer.java:94)
    at org.broadinstitute.gatk.engine.traversals.TraverseReadsNano$TraverseReadsMap.apply(TraverseReadsNano.java:228)
    at org.broadinstitute.gatk.engine.traversals.TraverseReadsNano$TraverseReadsMap.apply(TraverseReadsNano.java:216)
    at org.broadinstitute.gatk.utils.nanoScheduler.NanoScheduler.executeSingleThreaded(NanoScheduler.java:274)
    at org.broadinstitute.gatk.utils.nanoScheduler.NanoScheduler.execute(NanoScheduler.java:245)
    at org.broadinstitute.gatk.engine.traversals.TraverseReadsNano.traverse(TraverseReadsNano.java:102)
    at org.broadinstitute.gatk.engine.traversals.TraverseReadsNano.traverse(TraverseReadsNano.java:56)
    at org.broadinstitute.gatk.engine.executive.LinearMicroScheduler.execute(LinearMicroScheduler.java:108)
    at org.broadinstitute.gatk.engine.GenomeAnalysisEngine.execute(GenomeAnalysisEngine.java:319)
    at org.broadinstitute.gatk.engine.CommandLineExecutable.execute(CommandLineExecutable.java:121)
    at org.broadinstitute.gatk.utils.commandline.CommandLineProgram.start(CommandLineProgram.java:248)
    at org.broadinstitute.gatk.utils.commandline.CommandLineProgram.start(CommandLineProgram.java:155)
    at org.broadinstitute.gatk.engine.CommandLineGATK.main(CommandLineGATK.java:107)

    ERROR ------------------------------------------------------------------------------------------
    ERROR A GATK RUNTIME ERROR has occurred (version 3.3-0-g37228af):
    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: Code exception (see stack trace for error itself)
    ERROR ------------------------------------------------------------------------------------------

    When I checked the public source, I found that firstReadInPair(ReadAdaptorTrimmer.java:176) was NULL.

    By checking the value of readIn.getReadName()/firstReadInPair.getReadName(), I found that the read order is NOT _1 read , then _2 read, then another _1 read. But the source except read order is _1 read , then _2 read, then another _1 read.

    I thought it is a problem of ReadAdaptorTrimmer.java.

    When I tried to fix the problem myself, I found another exception
    [Tue Apr 28 20:44:57.752 2015] ##### ERROR MESSAGE: Alignments added out of order in SAMFileWriterImpl.addAlignment for /bio/hpc-bio/gatk/trim.bam. Sort order is coordinate. Offending records are at [1:258879] and [1:0]

    I found that this exception may be caused by the lack of setAlignmentStart()
    firstReadInPair.setAlignmentStart(readIn.getAlignmentStart());

    I'm sorry I failed to fix the problem finally.
    Please release a hofix of ReadAdaptorTrimmer.java if this is a problem of ReadAdaptorTrimmer.java

    Best Regards
    Wang Yugui

  • SheilaSheila Broad InstituteMember, Broadie ✭✭✭✭✭

    @wang_yugui
    Hi Wang,

    Can you run Picard's Validate Sam File and let me know what errors you get? http://broadinstitute.github.io/picard/command-line-overview.html#ValidateSamFile

    Thanks,
    Sheila

  • wang_yuguiwang_yugui Member

    Hi Sheila,

    result of picard ValidateSamFile:

    [[email protected] gatk]# java -jar /usr/hpc-bio/picard-tools/picard.jar ValidateSamFile I=/biowrk/bam.md/Project_14686/Sample_100T/md.bam
    [Wed Apr 29 00:10:21 CST 2015] picard.sam.ValidateSamFile INPUT=md.bam MODE=VERBOSE MAX_OUTPUT=100 IGNORE_WARNINGS=false VALIDATE_INDEX=true IS_BISULFITE_SEQUENCED=false MAX_OPEN_TEMP_FILES=8000 VERBOSITY=INFO QUIET=false VALIDATION_STRINGENCY=STRICT COMPRESSION_LEVEL=5 MAX_RECORDS_IN_RAM=500000 CREATE_INDEX=false CREATE_MD5_FILE=false
    [Wed Apr 29 00:10:21 CST 2015] Executing as [email protected] on Linux 3.10.0-229.1.2.el7.x86_64 amd64; OpenJDK 64-Bit Server VM 1.7.0_79-mockbuild_2015_04_15_00_14-b00; Picard version: 1.130(8b3e8abe25f920f5aa569db482bb999f29cc447b_1427207353) IntelDeflater
    No errors found
    [Wed Apr 29 00:10:22 CST 2015] picard.sam.ValidateSamFile done. Elapsed time: 0.03 minutes.
    Runtime.totalMemory()=2058354688
    [[email protected] gatk]#

    by the way, the md.sam is created by:
    1, head -n 10000 to the part of r1.fastq/r2.fastq
    2, bwa mem 0.7.12/samtools 1.2
    3, picard SortSam/MarkDuplicates 1.130

    I also tested GATK 2.6-4, the same error happened.

    java -XX:-UseCompressedOops -Xms100g -jar /usr/hpc-bio/gatk/GATK-2.6.jar -T ReadAdaptorTrimmer --minMatches 41 -R /usr/bio-ref/GRCh38.79/sequences.full.fa --bam_compression 0 --removeUnpairedReads -I md.bam -o trim.bam
    INFO 00:19:28,328 HelpFormatter - --------------------------------------------------------------------------------
    INFO 00:19:28,330 HelpFormatter - The Genome Analysis Toolkit (GATK) v2.6-4-g3e5ff60, Compiled 2013/06/24 14:48:56
    INFO 00:19:28,330 HelpFormatter - Copyright (c) 2010 The Broad Institute
    INFO 00:19:28,330 HelpFormatter - For support and documentation go to http://www.broadinstitute.org/gatk
    INFO 00:19:28,334 HelpFormatter - Program Args: -T ReadAdaptorTrimmer --minMatches 41 -R /usr/bio-ref/GRCh38.79/sequences.full.fa --bam_compression 0 --removeUnpairedReads -I /biowrk/bam.md/Project_14686/Sample_100T/md.bam -o trim.bam
    INFO 00:19:28,334 HelpFormatter - Date/Time: 2015/04/29 00:19:28
    INFO 00:19:28,334 HelpFormatter - --------------------------------------------------------------------------------
    INFO 00:19:28,334 HelpFormatter - --------------------------------------------------------------------------------
    INFO 00:19:28,429 GenomeAnalysisEngine - Strictness is SILENT
    INFO 00:19:28,532 GenomeAnalysisEngine - Downsampling Settings: No downsampling
    INFO 00:19:28,540 SAMDataSource$SAMReaders - Initializing SAMRecords in serial
    INFO 00:19:28,562 SAMDataSource$SAMReaders - Done initializing BAM readers: total time 0.02
    INFO 00:19:28,638 GenomeAnalysisEngine - Preparing for traversal over 1 BAM files
    INFO 00:19:28,642 GenomeAnalysisEngine - Done preparing for traversal
    INFO 00:19:28,643 ProgressMeter - [INITIALIZATION COMPLETE; STARTING PROCESSING]
    INFO 00:19:28,643 ProgressMeter - Location processed.reads runtime per.1M.reads completed total.runtime remaining
    INFO 00:19:28,645 ReadShardBalancer$1 - Loading BAM index data
    INFO 00:19:28,646 ReadShardBalancer$1 - Done loading BAM index data
    WARN 00:19:32,271 RestStorageService - Error Response: PUT '/GATK_Run_Reports/EKpuM5DQGgoF3nZPHND3PpW8p10UCOjK.report.xml.gz' -- ResponseCode: 403, ResponseStatus: Forbidden, Request Headers: [Content-Length: 756, Content-MD5: D6AS2VqYp05OnifQQITgAg==, Content-Type: application/octet-stream, x-amz-meta-md5-hash: 0fa012d95a98a74e4e9e27d04084e002, Date: Tue, 28 Apr 2015 16:19:29 GMT, Authorization: AWS AKIAIMHBU7X642TCHQ2A:245K0NUZ3J+EqGpHqgZA7ozjehs=, User-Agent: JetS3t/0.8.1 (Linux/3.10.0-229.1.2.el7.x86_64; amd64; en; JVM 1.7.0_79), Host: s3.amazonaws.com, Expect: 100-continue], Response Headers: [x-amz-request-id: 7C7EF783759AB9B7, x-amz-id-2: epZk3XmLnb21h0brjIk2ot0o+wYR0a8upW9kQcy3lXsJjci8B5uRe5ucNDIqEyaiaEajXMjzlFE=, Content-Type: application/xml, Transfer-Encoding: chunked, Date: Tue, 28 Apr 2015 16:21:48 GMT, Connection: close, Server: AmazonS3]

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

    java.lang.NullPointerException
    at org.broadinstitute.sting.gatk.walkers.readutils.ReadAdaptorTrimmer.map(ReadAdaptorTrimmer.java:173)
    at org.broadinstitute.sting.gatk.walkers.readutils.ReadAdaptorTrimmer.map(ReadAdaptorTrimmer.java:91)
    at org.broadinstitute.sting.gatk.traversals.TraverseReadsNano$TraverseReadsMap.apply(TraverseReadsNano.java:228)
    at org.broadinstitute.sting.gatk.traversals.TraverseReadsNano$TraverseReadsMap.apply(TraverseReadsNano.java:216)
    at org.broadinstitute.sting.utils.nanoScheduler.NanoScheduler.executeSingleThreaded(NanoScheduler.java:274)
    at org.broadinstitute.sting.utils.nanoScheduler.NanoScheduler.execute(NanoScheduler.java:245)
    at org.broadinstitute.sting.gatk.traversals.TraverseReadsNano.traverse(TraverseReadsNano.java:102)
    at org.broadinstitute.sting.gatk.traversals.TraverseReadsNano.traverse(TraverseReadsNano.java:56)
    at org.broadinstitute.sting.gatk.executive.LinearMicroScheduler.execute(LinearMicroScheduler.java:108)
    at org.broadinstitute.sting.gatk.GenomeAnalysisEngine.execute(GenomeAnalysisEngine.java:311)
    at org.broadinstitute.sting.gatk.CommandLineExecutable.execute(CommandLineExecutable.java:113)
    at org.broadinstitute.sting.commandline.CommandLineProgram.start(CommandLineProgram.java:245)
    at org.broadinstitute.sting.commandline.CommandLineProgram.start(CommandLineProgram.java:152)
    at org.broadinstitute.sting.gatk.CommandLineGATK.main(CommandLineGATK.java:91)

    ERROR ------------------------------------------------------------------------------------------
    ERROR A GATK RUNTIME ERROR has occurred (version 2.6-4-g3e5ff60):
    ERROR
    ERROR Please check the documentation guide 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: Code exception (see stack trace for error itself)
    ERROR ------------------------------------------------------------------------------------------

    Best Regards
    Wang Yugui

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @wang_yugui Would you be able to share a test case so that we can reproduce this locally? Instructions are here: https://www.broadinstitute.org/gatk/guide/article?id=1894

  • wang_yuguiwang_yugui Member

    Hi Sheila,Geraldine_VdAuwera

    I have upload ReadAdaptorTrimmer.Exception.by.wang.zip(md.bam/md.bai) to ftp://[email protected]/

    The reference genome is too big to upload, maybe you need picard SamToFastq to extract the seq, and then rebuild the bam locally.

    command usage:
    java -XX:-UseCompressedOops -Xms100g -jar /usr/hpc-bio/gatk/GenomeAnalysisTK.jar -T ReadAdaptorTrimmer --minMatches 41 -R /usr/bio-ref/GRCh38.79/sequences.full.fa --bam_compression 0 --removeUnpairedReads -I md.bam -o trim.bam

    Best Regards
    Wang Yugui

  • wang_yuguiwang_yugui Member

    Hi Sheila,

    for that md.bam, I added the following process to trace read info.
    System.out.println("readIn.getReadName():"+readIn.getReadName()+readIn.getReadName()+
    ":"+readIn.getReferenceName()+"/"+readIn.getAlignmentStart());

    the first exception will happen for:
    readIn.getReadName():HWI-7001456:548:C61LCANXX:2:1101:9201:2753HWI-7001456:548:C61LCANXX:2:1101:9201:2753:1/2181982
    readIn.getReadName():HWI-7001456:548:C61LCANXX:2:1101:3350:3582HWI-7001456:548:C61LCANXX:2:1101:3350:3582:1/2182144
    readIn.getReadName():HWI-7001456:548:C61LCANXX:2:1101:3350:3582HWI-7001456:548:C61LCANXX:2:1101:3350:3582:1/2182217
    readIn.getReadName():HWI-7001456:548:C61LCANXX:2:1101:9201:2753HWI-7001456:548:C61LCANXX:2:1101:9201:2753:1/2182248

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @wang_yugui I'm afraid we cannot accept a bug report that does not include an appropriate reference (if using something other than the standard human references). If your reference file is too large to upload, you can subset the reference and the bam file to just the contig where the bug occurs.

  • wang_yuguiwang_yugui Member

    Hi Sheila,Geraldine_VdAuwera

    By studying the public source of ReadAdaptorTrimmer.java , I found that ReadAdaptorTrimmer was designed for unsorted bam files(the raw output of bwa mem|samtools view, and without SortSam).

    With the flowing command usage and source patch, ReadAdaptorTrimmer works now.
    1, command usage
    java -XX:-UseCompressedOops -Xms100g -jar /usr/hpc-bio/gatk/GenomeAnalysisTK.jar -T ReadAdaptorTrimmer --minMatches 41 -R /usr/bio-ref/GRCh38.79/sequences.full.fa --bam_compression 0 --removeUnpairedReads -I PE.raw.bam -o trim.bam --unsafe
    a) --unsafe for unsorted bam without index.
    The read order of unsorted bam matches ReadAdaptorTrimmer's expect.
    b) the ReadAdaptorTrimmer.java seems NOT thread-safe, so no "-nct" option.

    2,source patch
    a) set the output bam as unsorted
    b) add the missing AlignmentStart info of read.

    --- org/broadinstitute/gatk/tools/walkers/readutils/ReadAdaptorTrimmer.java Wed Apr 29 12:41:12 2015
    +++ org/broadinstitute/gatk/tools/walkers/readutils/ReadAdaptorTrimmer.java Wed Apr 29 20:43:46 2015

    import com.google.java.contract.Ensures;
    import com.google.java.contract.Requires;
    import htsjdk.samtools.SAMFileWriter;
    +import htsjdk.samtools.SAMFileHeader;
    import org.apache.log4j.Logger;
    import org.broadinstitute.gatk.engine.walkers.NanoSchedulable;
    import org.broadinstitute.gatk.engine.walkers.PartitionBy;

    * @return SAMFileWriter, set to the BAM output file if the command line option was set, null otherwise
    */
    public SAMFileWriter reduceInit() {

    • out.getFileHeader().setSortOrder(SAMFileHeader.SortOrder.unsorted);
      return out;
      }


    firstReadInPair.setReadString(readIn.getReadString());
    firstReadInPair.setReadName(readIn.getReadName());
    firstReadInPair.setBaseQualities(readIn.getBaseQualities());

    • firstReadInPair.setAlignmentStart(readIn.getAlignmentStart());
      }
      else {
      if (!readIn.getReadName().equals(firstReadInPair.getReadName())) {

    I hope these info maybe useful.

    Best Regards
    Wang Yugui

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @wang_yugui Thanks for the additional details. Please feel free to submit a patch against the public repository. Unfortunately we cannot do anything without a full test case, and we cannot yet work with the Grch38 reference.

Sign In or Register to comment.