[BioC] BAM files to Genomic Ranges object

[Hervé Pagès]

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