Thanks Michael. Hope you have a Happy Thanksgiving! Jason On Wed, Nov 23, 2011 at 4:58 PM, Michael Lawrence <lawrence.michael at="" gene.com=""> wrote: > At this point, > > cvg.select <- subsetByOverlaps(as(cvg, "GRanges"), toshow) > > Should work, but it is not very efficient. We just need a RleViewsList -> > RangedData/GRanges coercion. On the TODO list. > > On Wed, Nov 23, 2011 at 9:35 AM, Jason Lu <jasonlu68 at="" gmail.com=""> wrote: >> >> Hi Michael, >> >> Thanks for the quick response. >> I tried the keepSeqlevels and your solution using the seqselect. It >> works! The only caveat is that the genomic coordinates are lost ?in >> the bedGraph. Here is what I have tried >> >> # >> toshow = keepSeqlevels(toshow,seqlevels(ss)) >> cvg = coverage(ss) >> cvg.select=seqselect(cvg,as(toshow,"RangesList")) >> export(cvg.select,"test.bedGraph") >> >> # output of test.bedGraph [genomic coordinates were lost] >> track name="R Track" type=bedGraph >> chr5 ? ?0 ? ? ? 103 ? ? 0 >> chr5 ? ?103 ? ? 104 ? ? 71 >> chr5 ? ?104 ? ? 105 ? ? 72 >> chr5 ? ?105 ? ? 106 ? ? 74 >> chr5 ? ?106 ? ? 109 ? ? 76 >> chr5 ? ?109 ? ? 110 ? ? 77 >> >> I may need to do shift with each exon start positions, or you may have >> more elegant solutions. >> Thanks. >> >> Jason >> >> >> On Wed, Nov 23, 2011 at 11:38 AM, Michael Lawrence >> <lawrence.michael at="" gene.com=""> wrote: >> > >> > >> > On Wed, Nov 23, 2011 at 7:18 AM, Jason Lu <jasonlu68 at="" gmail.com=""> wrote: >> >> >> >> I want to thank Michael and Martin for being helpful. >> >> >> >> My apology for not showing the full script last email. The purpose >> >> here is to get a figure showing read coverage in the exon regions of >> >> my gene of interest (one gene). I tried some additional steps here and >> >> get an error below. >> >> >> >> # >> >> nm = "NM_000997" >> >> txs ?= exonsBy(txdb,by="tx",use.name=TRUE) >> >> toshow = txs[[nm]] >> >> ss <- >> >> readBamGappedAlignments("myalign.bam",param=ScanBamParam(which=t oshow)) >> >> strand(ss) = "*" >> >> cvg ?= coverage(ss) ?# cvg: SimpleRleList >> >> >> >> cvg.view <- Views(cvg, as(toshow, "RangesList")) >> >> Error in .Method(..., FUN = FUN, MoreArgs = MoreArgs, SIMPLIFY = >> >> SIMPLIFY, >> >> ?: >> >> ?all object lengths must be multiple of longest object length >> >> >> > >> > The problem here is likely that "toshow" is on a wider universe (set of >> > sequences) than your coverage vector from the BAM file. You can use >> > keepSeqlevels(toshow, seqlevels(ss)). You want to make sure the >> > seqlevels >> > are in the same order; not sure if keepSeqlevels() does this -- perhaps >> > it >> > should? >> > >> >> >> >> Also, I have been unsuccessful in exporting the RleViewsList to a file >> >> (wig or bedGraph format). Could you give a hint on this as well? >> >> Thanks again. >> >> >> > >> > Now that I think about it, Views are a bit of an overkill for this >> > specific >> > case. You could just do a >> > >> > cvg.select <- seqselect(cvg, as(toshow, "RangesList")) >> > >> > And then export that: >> > >> > export(cvg.select, "coverage.bedGraph") >> > >> >> Jason >> >> >> >> >> >> On Wed, Nov 23, 2011 at 8:20 AM, Michael Lawrence >> >> <lawrence.michael at="" gene.com=""> wrote: >> >> > >> >> > >> >> > On Tue, Nov 22, 2011 at 9:39 AM, Jason Lu <jasonlu68 at="" gmail.com=""> >> >> > wrote: >> >> >> >> >> >> Hi list, >> >> >> >> >> >> I have a question and would like to get your help. >> >> >> What I would like to get is to show the read coverage in the exon >> >> >> regions (only) in the ucsc genome browser. My question is how to >> >> >> export the Views so I can load into ucsc browser. >> >> >> Here is what I have: >> >> >> >> >> >> # >> >> >> >cvg = coverage(ranges(ss)) >> >> > >> >> > Careful here: you are not distinguishing chromosomes when calculating >> >> > this >> >> > coverage. It's all going into one vector. Instead, you want to call >> >> > coverage >> >> > directly on your GRanges 'ss'. >> >> > >> >> > cvg <- coverage(ss) >> >> > >> >> >> >> >> >> >cvg.view = Views(cvg,rangesexon.gr)) ?# exon.gr is a GRanges >> >> > >> >> > Since 'cvg' is now an RleList (with one element per chromosome), you >> >> > want a >> >> > RangesList to parallel it in your RleViewsList. >> >> > >> >> > cvg.view <- Views(cvg, asexon.gr, "RangesList")) >> >> > >> >> > I think we should try to make that easier and have Views() take a >> >> > GRanges. >> >> > Will make a note. >> >> > >> >> >> >> >> >> >cvg.view >> >> >> Views on a 22109556-length Rle subject >> >> >> >> >> >> views: >> >> >> ? ? ? ?start ? ? ?end width >> >> >> ?[1] 22109317 22109688 ? 372 [5 5 7 5 6 6 6 5 5 5 5 5 4 4 5 8 8 8 9 >> >> >> 9 9 >> >> >> 9 >> >> >> ...] >> >> >> ?[2] 22084156 22084276 ? 121 [4 4 4 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 >> >> >> 4 4 >> >> >> 4 >> >> >> ...] >> >> >> ?[3] 22083039 22083195 ? 157 [2 1 1 0 0 0 1 1 1 1 2 2 2 2 2 2 2 2 2 >> >> >> 3 3 >> >> >> 4 >> >> >> ...] >> >> >> ?[4] 22079485 22079612 ? 128 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 2 2 2 >> >> >> 2 2 >> >> >> 2 >> >> >> ...] >> >> >> ?[5] 22079020 22079144 ? 125 [6 6 6 6 7 7 7 7 7 7 7 7 7 7 6 6 6 6 6 >> >> >> 6 5 >> >> >> 4 >> >> >> ...] >> >> >> >> >> >> I don't know how to modify this view such as adding chromosome info >> >> >> and export this 'cvg.view' in rtracklayer. >> >> > >> >> > >> >> > You can now export this RleViewsList directly to a file for upload. >> >> > For >> >> > upload to the UCSC browser, how about: >> >> > >> >> > session <- browserSession() >> >> > session$coverage <- cvg.view >> >> > browserView(session, exon.gr[1]) >> >> > >> >> > That should show you the first exon. >> >> > >> >> >> Thanks in advance. >> >> >> Jason >> >> >> >> >> >> > sessionInfo() >> >> >> R version 2.14.0 (2011-10-31) >> >> >> Platform: x86_64-unknown-linux-gnu (64-bit) >> >> >> >> >> >> locale: >> >> >> [1] C >> >> >> >> >> >> attached base packages: >> >> >> [1] stats ? ? graphics ?grDevices utils ? ? datasets ?methods ? base >> >> >> >> >> >> other attached packages: >> >> >> ?[1] rtracklayer_1.14.3 ? ?RCurl_1.7-0 ? ? ? ? ? bitops_1.0-4.1 >> >> >> ?[4] Rsamtools_1.6.2 ? ? ? Biostrings_2.22.0 >> >> >> GenomicFeatures_1.6.4 >> >> >> ?[7] AnnotationDbi_1.16.4 ?Biobase_2.14.0 ? ? ? ?GenomicRanges_1.6.3 >> >> >> [10] IRanges_1.12.2 ? ? ? ?BiocInstaller_1.2.1 >> >> >> >> >> >> loaded via a namespace (and not attached): >> >> >> [1] BSgenome_1.22.0 DBI_0.2-5 ? ? ? RSQLite_0.10.0 ?XML_3.4-3 >> >> >> [5] biomaRt_2.10.0 ?tools_2.14.0 ? ?zlibbioc_1.0.0 >> >> >> > >> >> >> >> >> >> _______________________________________________ >> >> >> Bioconductor mailing list >> >> >> Bioconductor at r-project.org >> >> >> https://stat.ethz.ch/mailman/listinfo/bioconductor >> >> >> Search the archives: >> >> >> http://news.gmane.org/gmane.science.biology.informatics.conductor >> >> > >> >> > >> > >> > > >