
Think there has been some discussion about trimming specifically on
Bioc-sig-seq (recently about 08/02/10) ago so you might want to look at
that. 

IN GENERAL you want to use a "iterate" fast  , use XStringViews,  then
it goes very fast!

For example for motif searching in many regions:

myb<-"YAACKG"         
myb.dna<-DNAString(myb)
myb.dna.comp.rev<-reverseComplement(myb.dna)

all.genomic<-getSeq(Mmusculus, the.chrom, starts, ends) # get some sequences can read fasta etc
length(all.genomic)

> [1] 81639
> class(all.genomic)
[1] "character"
> starts[1:5]
[1] 3187620 3487470 3777371 4144203 4274129
> ends[1:5]
[1] 3187780 3487630 3777531 4144363 4274289
> the.chrom[1:5]
[1] "chr1" "chr1" "chr1" "chr1" "chr1"


system.time(x<- XStringViews(all.genomic, "DNAString"))

x.labels<-paste(the.chrom,starts,ends,sep=":") # give the sequences/reads names if you want
names(x)<-x.labels

###################### forward counts #################
all.matches<-matchPattern(myb.dna,x,max.mismatch=0, fixed=FALSE) # needs a stringView to vectorize
the.cov<-coverage(all.matches)
counts<-aggregate(the.cov,start=start(x),end=end(x),FUN=sum)/length(myb.dna)
######################################################

###################### reverse counts #################
all.matches.comp<-matchPattern(myb.dna.comp.rev,x,max.mismatch=0, fixed=FALSE) # needs a stringView to vectorize
the.cov.comp<-coverage(all.matches.comp)
counts.comp<-aggregate(the.cov.comp,start=start(x),end=end(x),FUN=sum)/length(myb.dna)
######################################################

counts.t<-counts+counts.comp # all sequences/reads with a match

####################### forward posns #################
loc.f=rep(NA,times=dim(peaks)[1])
loc.hit<-start(all.matches)
region.start<-start(x)
hits<-grep(TRUE,counts>0) # places where a hits is found
k<-1
no.hits<-counts[counts>0]
for(i in 1:length(hits)){
  loc.f[hits[i]]<-toString((loc.hit[k:(k
+no.hits[i]-1)]-region.start[hits[i]])+1)
  k<-k+no.hits[i]
}
loc.f[1:50]
######################################################

###################### revearse posns #################
loc.c=rep(NA,times=dim(peaks)[1])
loc.hit<-start(all.matches.comp)
region.start<-start(x)
hits<-grep(TRUE,counts.comp>0) # places where a hits is found
k<-1
no.hits<-counts.comp[counts.comp>0]
for(i in 1:length(hits)){
  loc.c[hits[i]]<-toString((loc.hit[k:(k
+no.hits[i]-1)]-region.start[hits[i]])+1)
  k<-k+no.hits[i]
}
######################################################

loc.f 
and loc.c
are locations of the pattern w.r.t forward direction

Any If you want to do trimming see the bioc-sig-seq discussion
too....too they probably have a more elegant solution

Cheers
Paul







-----Original Message-----
From: Paul Shannon <pshannon@systemsbiology.org>
To: bioc <bioconductor@stat.math.ethz.ch>
Subject: [BioC] fast iterator over DNAString's?
Date: Wed, 10 Mar 2010 16:30:40 -0800


I wish to trim a variable length sequence from the end of many thousands of DNAStrings in a DNAStringSet.  

The sequence to be trimmed is any recognizable chunk of a solexa short read adapter, which ends up on the end of, for example, 22nt miRNAs.  The adapter chunk might be found in the middle of a 35 base read, or it might be closer to the end.  In every case, I want to delete every base from the start of the adapter chunk to the end of the read.

I imagine there might be a BString operation equivalent to sed.  See could be used ike this:

  echo 'CGAAGCGGGATGATCTATCTCGTATGCCGTCTTCT' | sed s/TCGTATGCCGTC.*$//      --> GAAGCGGGATGATCTATC

(where TCGTATGCCGTC is only part of the 21-base adapter, but is probably a long enough portion to be representative)

Any way to do this with BStrings and friends?

Thanks!

 - Paul
_______________________________________________
Bioconductor mailing list
Bioconductor@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/bioconductor
Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor

	[[alternative HTML version deleted]]

