[BioC] export Views of Rle in rtracklayer
Jason Lu
jasonlu68 at gmail.com
Wed Nov 23 18:35:56 CET 2011
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