[BioC] BAM files to Genomic Ranges object

[jluis.lavin at unavarra.es]

Thank you very much Hervé,

Your code really got the job done.

And thanks also to Steve for his nice comments and Tim for his kind
suggestions.

With best wishes

JL


El Mie, 23 de Enero de 2013, 7:40, Hervé Pagès escribió:
> Hi,
>
> On 01/22/2013 08:29 AM, Steve Lianoglou wrote:
>> Hi,
>>
>> On Tue, Jan 22, 2013 at 10:54 AM,  <jluis.lavin at unavarra.es> wrote:
>>> Hello Steve,
>>>
>>> (1) What is the output of `is(cpdens)` ? If it's not a data.frame, this
>>> won't work
>>>
>>> is(cpdens)
>>> [1] "list"           "vector"         "AssayData"      "vectorORfactor"
>>>
>>> head(cpdens)
>>> .
>>> .
>>> .
>>>
>>> [99841]  1  3  3  3  4  4  2  2  2  1  0  0  0  0  0  1  1  1  1  1  1
>>> 1
>>> 0  0
>>> [99865]  2  2  0  0  0  0  0  0  0  3  1  1  0  0  0  0  0  0  0  2  1
>>> 1
>>> 1  0
>>> [99889]  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1  1
>>> 0
>>> 0  0
>>> [99913]  0  3  3  3  3  3  3  3  0  0  0  3  1 13  0  3  2  2  1  0  0
>>> 0
>>> 1  1
>>> [99937]  0  0  0  2  2  2  2  3  3  1  1  1  1  1  0  0  0  0  0  0  0
>>> 2
>>> 0  2
>>> [99961]  0  1  0  0  2  2  2  0  0  0  0  0  0  0  0  0  0  0  0  0  1
>>> 1
>>> 1  1
>>> [99985]  1  0  0  0  0  0  0  1  0  0  0  0  0  0  0
>>>   [ reached getOption("max.print") -- omitted 3301343 entries ]
>>>
>>> It doesn't seem to be data.frame as you pointed out...
>>> Is it possible write this to an output file of some kind?
>>
>> Yes, of course.
>>
>> You just have to ask yourself how you want to serialize it.
>>
>> ?writeLines can be your friend, here, if you want to save it in some
>> ASCII format.
>
> Or you can use some higher-level facilities like export() from the
> rtracklayer package to export to a BED file. That means that you
> first need to stick the numeric vectors of CpG densities back into
> the GRanges object they correspond to ('cpdens' should have the same
> shape as the original GRangesList you passed to cpgDensityCalc()
> i.e. it should be list of numeric vectors where the i-th list
> element has the length of the i-th GRangesList element). Then
> you export each of the GRanges object separately. Could be done
> with something like:
>
>      library(rtracklayer)
>      for (i in seq_along(gr_list)) {
>          gr <- gr_list[[i]]
>          ## Need to set the "score" metadata col (see ??export.bed).
>          mcols(gr)$score <- cpdens[[i]]
>          filepath <- paste("CpG", i, ".bed", sep="")
>          export.bed(gr, filepath)
>      }
>
> Granted that 'gr_list' was obtained by loading BAM files, you should
> end up with 1 BED file per original BAM file.
>
> HTH,
> H.
>
>>
>> -steve
>>
>
> --
> Hervé Pagès
>
> Program in Computational Biology
> Division of Public Health Sciences
> Fred Hutchinson Cancer Research Center
> 1100 Fairview Ave. N, M1-B514
> P.O. Box 19024
> Seattle, WA 98109-1024
>
> E-mail: hpages at fhcrc.org
> Phone:  (206) 667-5791
> Fax:    (206) 667-1319
>


-- 
Dr. José Luis Lavín Trueba

Dpto. de Producción Agraria
Grupo de Genética y Microbiología
Universidad Pública de Navarra
31006 Pamplona
Navarra
SPAIN