[BioC] Any package to calculate NGS coverage and plot it?

[Martin Morgan]

On 05/28/2011 06:24 AM, Steve Lianoglou wrote:
> Hi,
>
> On Fri, May 27, 2011 at 9:24 PM, Xiaohui Wu<wux3 at muohio.edu>  wrote:
>> Dear list,
>>
>> I have some NGS sequences that enriched in different regions along
>> the genome. I want to calculate the coverage of the short sequences
>> and plot the sequences  in a specific region. To plot, I just want
>> a function to simply plot the tags and chromosome (it is OK if also
>> plot the gene there), but not other things like in genome browser.
>> Is there any package can do this?

A simple approach is to calculate coverage

   library(GenomicRanges)
   aln = readGappedAlignments(<bamfile>)
   cvg = coverage(aln)

then slice to find views onto islands

   isles = slice(cvg, <min-depth>)

or view selected coordinates

   v = Views(cvg, <interesting-GRanges>)

then visualize as a simple 'skyline', e.g.,

   plot(as.integer(isles[[1]]))

or as a 'movie'

   for (i in seq_len(length(isles)))
       plot(as.integer(isles[[i]]), type="l")

None of this tested, so hopefully not too misleading.

Martin

>
> You can try GenomeGraphs as a start:
> http://www.bioconductor.org/packages/release/bioc/html/GenomeGraphs.html
>
>  With a bit of work I reckon you can get what you want. Read through
> the vignette -- last I remember there was an NGS example there.
>
> I don't imagine you'll get very good results when you try to plot an
> entire genome, so you'll have to identify your regions of interest
> first and plot them individually (using slice() over your coverage
> vector would be an easy first-attempt to find them).
>
> -steve
>


-- 
Computational Biology
Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N. PO Box 19024 Seattle, WA 98109

Location: M1-B861
Telephone: 206 667-2793