On Wed, Oct 16, 2013 at 1:38 AM, Lukasz Kielpinski wrote: > It worked, the 'map' function was very helpful! > > In the end I implemented as you said with modifications: > a) reassigned "chrM" is_circular to FALSE for mapping to work (I still use > R 2.15) > > my_exons@unlistData@seqinfo@is_circular <- rep(FALSE, 35) > > Please refrain from direct slot access. Something like isCircular(seqinfo(my_exons)) <- FALSE > b) reassigned start positions of motifsTX to 1 is below 1: > > start(motifsTX)[start(motifsTX)<=1] <- 1 > > or: restrict(motifsTX, 1) > c) extracted sequences without ordering by mapping (motifsTX was already > ordered) and using subseq function (" '[' subsetting by Ranges is defunct ") > > motifSeqs <- subseq(tx[subjectHits(mapping)], start=start(motifsTX), > end=end(motifsTX)) > > Yea, the Ranges support for [ was re-incarnated last devel cycle. You really should update. Thanks a lot for the valuable hint! > > Lukasz > > > > > > On 15 October 2013 19:00, Michael Lawrence wrote: > >> Use map(gr, exonsByResult) to map your GRanges to transcript coordinates. >> Then, do any range transformations like flank(). Use >> GenomicFeatures::extractTranscriptsFromGenome to get the transcript >> sequences. To extract those ranges from each sequence, use [, with some >> repping. Should only be about 4 lines of code. Something like this: >> >> mapping <- map(motifsGR, exonsBy(TxDb.Hsapiens.UCSC.hg19.knownGene)) >> motifsTX <- flank(ranges(mapping), n) >> tx <- extractTranscriptsFromGenome(Hsapiens, >> TxDb.Hsapiens.UCSC.hg19.knownGene) >> motifSeqs <- tx[subjectHits(mapping)][motifsTX[queryHits(mapping)]] >> >> Untested... >> >> Good luck, >> Michael >> >> >> On Tue, Oct 15, 2013 at 9:36 AM, Lukasz [guest] wrote: >> >>> >>> Hi! >>> >>> Problem summary: How to retrieve part of the sequence of mRNA around >>> given location. >>> >>> I have the locations of the binding to mRNA events as GRanges (GRevents) >>> and need to retrieve sequence for motif finding. The problem is that if I >>> use getSeq(flank(GRevents, width=n)) then I get the genomic sequence not >>> transcript sequence, i.e. not accounting for introns or mRNA border. I have >>> tried solving it with exonsBy("transcriptDb object", "tx") function but >>> without success. >>> >>> Question: Is there a bioconductor-supported way of getting resolving the >>> problem? With CLIPseq being more and more popular this will be very >>> demanded function. >>> >>> Thanks, >>> Lukasz >>> >>> -- output of sessionInfo(): >>> >>> > sessionInfo() >>> R version 2.15.0 (2012-03-30) >>> Platform: x86_64-unknown-linux-gnu (64-bit) >>> >>> 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=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 >>> [7] LC_PAPER=C 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] rtracklayer_1.16.0 GenomicFeatures_1.8.1 >>> [3] AnnotationDbi_1.18.4 Biobase_2.16.0 >>> [5] BSgenome.Mmusculus.UCSC.mm9_1.3.17 BSgenome_1.24.0 >>> [7] Biostrings_2.24.1 GenomicRanges_1.8.3 >>> [9] IRanges_1.14.2 BiocGenerics_0.2.0 >>> >>> loaded via a namespace (and not attached): >>> [1] biomaRt_2.12.0 bitops_1.0-4.1 DBI_0.2-5 RCurl_1.91-1 >>> [5] Rsamtools_1.8.0 RSQLite_0.11.1 stats4_2.15.0 tools_2.15.0 >>> [9] XML_3.9-4 zlibbioc_1.2.0 >>> >>> >>> -- >>> Sent via the guest posting facility at bioconductor.org. >>> >>> _______________________________________________ >>> Bioconductor mailing list >>> Bioconductor@r-project.org >>> https://stat.ethz.ch/mailman/listinfo/bioconductor >>> Search the archives: >>> http://news.gmane.org/gmane.science.biology.informatics.conductor >>> >> >> > [[alternative HTML version deleted]]