[Bioc-sig-seq] trimLRPatterns question
Harris A. Jaffee
hj at jhu.edu
Wed Feb 23 07:22:25 CET 2011
You hint that your trimLen is > 50, yet when you show your
max.Rmismatch, it's clearly of length 47, composed
of rep(1,12) and rep(2,N) with N=trimLen-15. So N=35 and
trimLen=50. A max.Rmismatch of length 47 (= 50-3)
will be augmented by rep(-1,3) at the bottom, preventing any matches
of an adapter prefix of length 3 or less
with a read suffix. This is an effect that you want, but to allow
more mismatch tolerance for longer adapter
prefixes, against read suffixes, your max.Rmismatch vector must be
monotone increasing, even if not strictly,
rather than decreasing. So your max.Rmismatch is at least reversed.
Its max (last value) might also be > 2.
Perhaps like this:
> ceiling(1:47 * 1/12)
[1] 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3
3 3 3 3 4 4
[39] 4 4 4 4 4 4 4 4 4
This will be better than any max.Rmismatch=e, with e in (0,1).
-Harris
On Feb 22, 2011, at 11:10 PM, Kunbin Qu wrote:
> Hi,
>
> I have a 50 cycle RNA-seq run, with variable length of adapters in
> some reads. The length of the adapter is longer than 50 bp, and the
> portion in each read could be as high as the full length of 50
> cycles, as well as zero bp. When I trim the adapter, I do not want
> over trim, so the max.Rmismatch parameter I used was shorter than
> (eg, 3 bp shorter) than the full length, also with more mismatch
> tolerance for longer portion, less mismatch for short segment. So I
> wrote like the following:
>
> trimLen<-length(adapter)
> trimCoords<-trimLRPatterns(Rpattern=adapter, subject=seqs,
> max.Rmismatch = c(rep(2,trimLen-15), rep(1, 12)), ranges=T,
> with.Rindels=TRUE)
>
> the value for the max.Rmismatch is basically:
>> c(rep(2,trimLen-15), rep(1, 12))
> [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
> 2 2 2 2 1 1 1
> [39] 1 1 1 1 1 1 1 1 1
>
> So I want to leave the last 3 bp alone (to avoid over-trimming by
> chance), then followed by 12 bp with 1 mismatch, then 35 bp with 2
> mismatch. Somehow when I tried this to the read, it does not seem
> work, the trimmed reads really do not seem what they should be.
> Anyone can help me on this? Thanks
>
> -Kunbin
>
>
>> sessionInfo()
> R version 2.11.0 (2010-04-22)
> x86_64-unknown-linux-gnu
>
> locale:
> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
> [5] LC_MONETARY=C LC_MESSAGES=en_US.UTF-8
> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
> [9] LC_ADDRESS=C LC_TELEPHONE=C
> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
> other attached packages:
> [1] ShortRead_1.6.2 Rsamtools_1.0.1 lattice_0.19-11
> [4] Biostrings_2.16.7 GenomicRanges_1.0.1 IRanges_1.6.8
>
> loaded via a namespace (and not attached):
> [1] Biobase_2.8.0 grid_2.11.0 hwriter_1.3 tools_2.11.0
>>
>
> ______________________________________________________________________
> The contents of this electronic message, including any attachments,
> are intended only for the use of the individual or entity to which
> they are addressed and may contain confidential information. If you
> are not the intended recipient, you are hereby notified that any
> use, dissemination, distribution, or copying of this message or any
> attachment is strictly prohibited. If you have received this
> transmission in error, please send an e-mail to
> postmaster at genomichealth.com and delete this message, along with
> any attachments, from your computer.
>
> _______________________________________________
> Bioc-sig-sequencing mailing list
> Bioc-sig-sequencing at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing
More information about the Bioc-sig-sequencing
mailing list