ClipReads

Hi,

I am getting an error when trying to soft clip bases with ClipReads --clipSequence that the sequences I am trying to clip are not at the end of a read. I can see this is mentioned in the documentation as being the correct behaviour. However, I would would still like to get a sofclipped BAM as output

Would it be possible to add a --validation_strictness LENIENT option which would allow this error to be overridden? i.e. ignore any matched sequences not at the ends of a read but still clip matches at the beginning of reads?

Thanks

Chris

Best Answer

Answers

  • CarneiroCarneiro Charlestown, MAMember admin

    you can only softclip bases at the tails of a read. What exactly are you trying to do with your reads?

  • cboustredcboustred Member

    Hi, thanks for getting back to me.

    I am analysing FASTQ files generated by a MiSeq TruSeq Custom Amplion experiment (http://www.illumina.com/products/truseq_custom_amplicon.ilmn).

    I would like to align the reads to the human reference genome, create a BAM, then clip the primer sequences so that no variant calls are identified within them.

    Any advice on this would be much appreciated

    Thanks again

    Chris

  • cboustredcboustred Member

    OK, thanks for the advice

    BW

    Chris

  • wchenwchen Member

    Why? The primer sequences should match the human reference genome and be helpful in alignment too, am I right?

    I am also having FASTQ files generated by MiSeq TruSeq Custon Amplicon experiment. I aligned, realigned and recalibrated the BAM files without marking duplication, and tried to soft-clip the primer sequences (first 22 nt here) using:

    -T ClipReads -R ...../ftp.broadinstitute.org/bundle/2.8/hg19/ucsc.hg19.fasta -I recal_merged3.3.bam -o clipped.bam -CT 1-22 -CR SOFTCLIP_BASES -os clip.stat --validation_strictness LENIENT

    But got errors:

    ERROR MESSAGE: Read Clipper cannot soft clip unmapped reads

    Bu the reads were mapped indeed, am I right? Thanks!

    @Carneiro said:
    You should clip the primer sequences before alignment, otherwise your alignment would be screwed up by them.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @wchen The primer sequences include adapters and other additional sequence that won't match the reference. All of that should be clipped off before alignment.

    The error may be happening if you're running on all contigs in the reference build, which includes unmapped reads. Are you using Queue by any chance?

  • wchenwchen Member

    @Geraldine_VdAuwera said:
    The primer sequences include adapters and other additional sequence that won't match the reference. All of that should be clipped off before alignment.

    The error may be happening if you're running on all contigs in the reference build, which includes unmapped reads. Are you using Queue by any chance?

    The adapter sequence should be trimmed definitely, but the primer sequence will help with alignment.

    I am not using Queue, just GATK3.3 HaplotypeCaller. Thanks!

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    I'm not sure why the tool is operating on unmapped reads, but you can resolve that by providing a list of target contigs with the -L argument. That should solve your problem.

  • wchenwchen Member

    @Geraldine_VdAuwera said:
    I'm not sure why the tool is operating on unmapped reads, but you can resolve that by providing a list of target contigs with the -L argument. That should solve your problem.

    Thanks. But I still got errors after including a bed file (with columns: chr start end):

    INFO 15:40:27,219 HelpFormatter - Program Args: -T ClipReads -R
    ucsc.hg19.fasta -I recal_3.3.bam -o recal_3.3_clipped.bam -L manifest.bed -CT 1-22 -CR SOFTCLIP_BASES -os clip.stat --validation_strictness LENIENT
    INFO 15:40:27,223 HelpFormatter - Executing as ... on Linux 2.6.18-194.el5 amd64; Java HotSpot(TM) 64-Bit Server VM 1.7.0-b147.
    ...
    INFO 15:40:28,233 SAMDataSource$SAMReaders - Initializing SAMRecords in serial
    INFO 15:40:28,294 SAMDataSource$SAMReaders - Done initializing BAM readers: total time 0.06
    INFO 15:40:28,459 IntervalUtils - Processing 154829 bp from intervals
    INFO 15:40:28,543 GenomeAnalysisEngine - Preparing for traversal over 1 BAM files
    INFO 15:40:28,601 GenomeAnalysisEngine - Done preparing for traversal
    INFO 15:40:28,602 ProgressMeter - [INITIALIZATION COMPLETE; STARTING PROCESSING]
    INFO 15:40:28,602 ProgressMeter - | processed | time | per 1M | | total | remaining
    INFO 15:40:28,603 ProgressMeter - Location | reads | elapsed | reads | completed | runtime | runtime
    INFO 15:40:28,603 ClipReads - Creating cycle clipper 0-21
    INFO 15:40:28,605 ReadShardBalancer$1 - Loading BAM index data
    INFO 15:40:28,625 ReadShardBalancer$1 - Done loading BAM index data
    INFO 15:40:30,129 ReadShardBalancer$1 - Loading BAM index data
    INFO 15:40:30,152 ReadShardBalancer$1 - Done loading BAM index data
    INFO 15:40:58,619 ProgressMeter - chr1:27534037 0.0 30.0 s 49.6 w 0.2% 3.5 h 3.5 h

    ERROR ------------------------------------------------------------------------------------------
    ERROR A USER ERROR has occurred (version 3.3-0-g37228af):

    ...

    ERROR MESSAGE: Read Clipper cannot soft clip unmapped reads
  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin
    edited February 2015

    @wchen
    Hi,

    Is there a reason you are setting --validation_strictness to LENIENT? There may be issues with your input files. We cannot help you until you have validated your files. https://www.broadinstitute.org/gatk/guide/tooldocs/org_broadinstitute_gatk_engine_CommandLineGATK.php#--validation_strictness

    -Sheila

  • wchenwchen Member

    I used the same input file and ran HaplotypeCaller without any problem. Changed to -validation_strictness to STRICT, and still got the same error.

    Also tried to include primer sequence in fasta format, it ran through initially but also stopped due to error:
    -T ClipReads -R ucsc.hg19.fasta -I ../recal3.3.bam -XF seqsToClip.fasta -o recal3.3_clipped2.bam -L manifest.bed -CR SOFTCLIP_BASES -os clip2.stat --validation_strictness STRICT

    INFO 17:24:50,463 ProgressMeter - [INITIALIZATION COMPLETE; STARTING PROCESSING]
    INFO 17:24:50,463 ProgressMeter - | processed | time | per 1M | | total | remaining
    INFO 17:24:50,464 ProgressMeter - Location | reads | elapsed | reads | completed | runtime | runtime
    ...

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

    java.lang.RuntimeException: Cannot apply soft clipping operator to the middle of a read: M00334:21:000000000-A36BV:1:1112:16743:24319 to be clipped at 27-53
    at org.broadinstitute.gatk.utils.clipping.ClippingOp.apply(ClippingOp.java:136)
    at org.broadinstitute.gatk.utils.clipping.ReadClipper.clipRead(ReadClipper.java:157)
    at org.broadinstitute.gatk.tools.walkers.readutils.ClipReads.reduce(ClipReads.java:471)
    at org.broadinstitute.gatk.tools.walkers.readutils.ClipReads.reduce(ClipReads.java:162)
    at org.broadinstitute.gatk.engine.traversals.TraverseReadsNano$TraverseReadsReduce.apply(TraverseReadsNano.java:251)
    at org.broadinstitute.gatk.engine.traversals.TraverseReadsNano$TraverseReadsReduce.apply(TraverseReadsNano.java:240)
    at org.broadinstitute.gatk.utils.nanoScheduler.NanoScheduler.executeSingleThreaded(NanoScheduler.java:279)
    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: Cannot apply soft clipping operator to the middle of a read: M00334:21:000000000-A36BV:1:1112:16743:24319 to be clipped at 27-53

    Thanks!

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @wchen

    Hi,

    It looks like the adapter sequence was found in the middle of the read instead of the very beginning. The program cannot deal with that; there should be nothing before the sequence provided in the -XF argument.

    -Sheila

  • wchenwchen Member

    @Sheila

    Thanks, that makes sense. But the sequence 27-53 AAAGTAAAACTTGCTGCTGATGAA in this sequence causing error doesn't not matching any primer sequences to be clipped from:
    M00334:21:000000000-A36BV:1:1112:16743:24319 161 chr1 27533857 60 175M = 27533986 304 GAGGTGGTAGCAAGGTTCCACAGGTAAAAGTAAAACTTGCTGCTGATGAAGATGATGATGAT....

    Does it cause error if there is a partial match? Thanks!

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @wchen No, it does not recognize partial matches. But the program checks both forward and reverse strand so there may be a match that you're not seeing because it's in the rev-complement sequence.

  • LMKellyLMKelly PittsburghMember

    @wchen @cboustred Where you able to get this to clip off the primer sequences only at the ends? We are also trying to do the same thing for the MiSeq data from an Illumina TruSight panel (same chemistry as the TruSeq custom panel I believe).

    Illumina uses the primer (oligo) sequences to align the reads before clipping them off as well. Although the Oligo sequences match the reference sequence, they are 100% synthetic and should not be included in the BAM file for sample reads.The sequencing adapters and barcodes have already been trimmed at this point so leaving the 'primer' or oligo (as Illumina calls them) sequences would not interfere with the alignment as suggested above.

Sign In or Register to comment.