[BioC] export Views of Rle in rtracklayer
Jason Lu
jasonlu68 at gmail.com
Thu Nov 24 01:19:28 CET 2011
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=toshow))
>> >> 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,ranges(exon.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, as(exon.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
>> >> >
>> >> >
>> >
>> >
>
>
More information about the Bioconductor
mailing list