Hi,
I have some 454 data which is expected to contain a primer at the beginning.

I noticed that using trimLRPatterns to trim an Lpattern at the 5' start of a
sequence with 1 mismatch does not allow matches to the subject sequence up
to the number of mismatches directly to the right of the sequence start (and
also to the left of an Rpattern at the 3' end).

So if my example primer Lpattern is 23 bases in length, trimLRPatterns with
1 mismatch will match all sequences at positions 0 to 22 and 1 to 23, but
not 2 to 24 whereas if you use vmatchPattern, it will find all three of
these positions.

My question is should an option/arg in trimLRPatterns be made available that
allows matches to the right of the Lpattern up to the number of mismatches,
and to the left of the Rpattern up to the number of mismatches?

I think the arg e.g. 'with.Lindels=TRUE' (which appears to have not been
enabled yet in Biostrings_2.14.8) may partially resolve this if the
insertion is somewhere at the beginning of a subject sequence  (within the
number of mismatches allowed) which can then be matched by the Lpattern.


> mismatches <- 1
> pattern    <- "AAGCAGTGGTATCAACGCAGAGT"
> w          <- width(pattern)
> mm         <- rep(mismatches, w)
> base       <- "G"
> n          <- 20
>
> subjectList <- list(
+                     "Primer0+poly(T)" = paste(substring(pattern,2,w),
polyn(base, n), sep=""),
+                     "Primer1+poly(T)"  = paste(pattern, polyn(base, n),
sep=""),
+                     "Primer2+poly(T)" = paste("A", substring(pattern,1,w),
polyn(base,n), sep="")
+                     )
>
> subjectSet <- DNAStringSet(unlist(subjectList))
>
> print(subjectSet)
  A DNAStringSet instance of length 3
    width seq                                               names

[1]    42 AGCAGTGGTATCAACGCAGAGTGGGGGGGGGGGGGGGGGGGG        Primer0+poly(T)
[2]    43 AAGCAGTGGTATCAACGCAGAGTGGGGGGGGGGGGGGGGGGGG       Primer1+poly(T)
[3]    44 AAAGCAGTGGTATCAACGCAGAGTGGGGGGGGGGGGGGGGGGGG      Primer2+poly(T)
> cat("Primer:  ", pattern, "\n")
Primer:  AAGCAGTGGTATCAACGCAGAGT
>
> ## trimLRPatterns -LPattern
> LRcoords     <- trimLRPatterns(Lpattern = pattern, subject = subjectSet,
 max.Lmismatch=mm, ranges=TRUE, with.Lindels=FALSE)
> ## Fails to match subjectSet[3:4] - trimming with mismatches is to the
left of Lpattern (and to the right of Rpattern)
> DNAStringSet(subjectSet,start(LRcoords), end(LRcoords))
  A DNAStringSet instance of length 3
    width seq                                               names

[1]    20 GGGGGGGGGGGGGGGGGGGG                              Primer0+poly(T)
[2]    20 GGGGGGGGGGGGGGGGGGGG                              Primer1+poly(T)
[3]    43 AAGCAGTGGTATCAACGCAGAGTGGGGGGGGGGGGGGGGGGGG       Primer2+poly(T)
>
> ## vmatchPattern
> tmp <- vmatchPattern(pattern, subjectSet, max.mismatch=mismatches)
> matchIndex <- which(as.logical(countIndex(tmp)))
> ## only one match per sequence so indices preserved with unlist
> VMcoords <- unlist(tmp)
> ## Matches to subjectSet[3]
> DNAStringSet(subjectSet[matchIndex],end(VMcoords)+1,
width(subjectSet)[matchIndex])
  A DNAStringSet instance of length 3
    width seq                                               names

[1]    20 GGGGGGGGGGGGGGGGGGGG                              Primer0+poly(T)
[2]    20 GGGGGGGGGGGGGGGGGGGG                              Primer1+poly(T)
[3]    20 GGGGGGGGGGGGGGGGGGGG                              Primer2+poly(T)

cheers,


 Marcus


> sessionInfo()
R version 2.10.0 (2009-10-26)
powerpc-apple-darwin8.11.1

locale:
[1] C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] hgu95av2probe_2.5.0 AnnotationDbi_1.8.1 Biobase_2.6.1
[4] ShortRead_1.4.0     lattice_0.17-26     BSgenome_1.14.2
[7] Biostrings_2.14.8   IRanges_1.4.9

loaded via a namespace (and not attached):
[1] DBI_0.2-4     RSQLite_0.7-3 grid_2.10.0   hwriter_1.1   tools_2.10.0

	[[alternative HTML version deleted]]

