[BioC] export Views of Rle in rtracklayer

[Jason Lu]

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

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.

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
>
>