On Wed, Sep 30, 2009 at 3:15 PM, Michael Lawrence <mflawren@fhcrc.org>wrote:

>
>
> On Wed, Sep 30, 2009 at 11:20 AM, Chris Seidel <seidel@phaget4.org> wrote:
>
>> Hi Michael,
>>
>> Thanks for the feedback. It's helpful.
>>
>> I think for ChIP Seq analysis a common activity will be going from
>> mapped reads to a trackfile, with some exploration of normalization,
>> peaks, smoothing, etc, in between. So I'm trying to figure out the
>> simplest way to do this, and traverse all the objects correctly.
>>
>> -----Original Message-----
>> From: lawremi@gmail.com [mailto:lawremi@gmail.com] On Behalf Of Michael
>> Lawrence
>> Sent: Wednesday, September 30, 2009 1:22 AM
>> To: Seidel, Christopher
>> Cc: bioc-sig-sequencing@r-project.org
>> Subject: Re: [Bioc-sig-seq] WIG from coverage
>>
>> On Tue, Sep 29, 2009 at 7:23 PM, Michael Lawrence <mflawren@fhcrc.org>
>> wrote:
>>
>> >> ip1.cov <- coverage(ip1, width=230208)
>> >> export(ip1.cov, "coverage.wig")
>> >> Am I doing something wrong?
>>
>> > Note that in general your code will not work, because you do not tell
>> rtracklayer the name of the chromosome. This makes me worry that
>> generating a default chromosome name will often hide user errors.
>>
>> I would have specified the chromosome if I knew it was necessary and
>> could have figured out how to do it. The current worflow goes from:
>>
>> export.txt file -> readAligned Object (one lane) -> GenomeData Object
>> (one lane) -> GenomeDataList (2 lanes) -> IRanges (from extendReads) ->
>> Rle (from coverage).
>>
>> To go from coverage to a track it would seem I have to make a RangedData
>> object from the coverage result, rather than use the result directly, is
>> that right? When I use the (modified) workflow below, now I get the
>> chromosome names correct:
>>
>> # read in some IP data
>> aln <- readAligned(spath, pattern="s_1_export.txt", type="SolexaExport",
>> filter=filt)
>> # convert to GenomeData object
>> s1 <- as(aln,"GenomeData")
>>
>> # read in some WCE data
>> aln <- readAligned(spath, pattern="s_3_export.txt", type="SolexaExport",
>> filter=filt)
>> # convert to GenomeData object
>> s3 <- as(aln,"GenomeData")
>>
>> # create GenomeDataList from ip and wce
>> ypd <- GenomeDataList(list(ip=s1, wce=s3))
>>
>> # extend reads
>> ip1 <- extendReads(ypd$ip$chrI, seqLen=150)
>> wce1 <- extendReads(ypd$wce$chrI, seqLen=150)
>>
>>
> The most convenient path would be to perform the analysis over all
> chromosomes and lanes at once. The extendReads() function should do this for
> you, e.g:
>
> extendedReads <- extendReads(ypd, seqLen=150)
>
> The 'extendedReads' above is a GenomeDataList. You could then do something
> like:
>
> cov <- gdapply(extendedReads, coverage)
>
> Now 'cov' is another GenomeDataList. Just subset out one lane and export
> it:
>
> export(cov$wce, "wce.wig")
>
>
>
>
>> # find coverage for chromosome 1
>> ip1.cov <- coverage(ip1, width=230208)
>> wce1.cov <- coverage(wce1, width=230208)
>>
>> # make a track for just chrI of IP
>> range1 <- IRanges(start=start(ip1.cov), end=end(ip1.cov))
>> score <- runValue(ip1.cov)
>> chrom <- rep("chrI", length(range1))
>> # put all these things together in a GenomicData object for rtracklayer
>> export
>> x <- GenomicData(range1, score, chrom=chrom, genome="sacCer2")
>>
>
> The coerce method does most of this for you, as in:
> x <-  as(ip1.cov, "RangedData")
>
> Now you just need to set the chromosome:
> names(x) <- "chr1"
>
> And then export.
>
>
>> export(x, "ip1_coverage.wig")
>>
>> So now this works to give chromosome names:
>>
>> track name="R Track" type=wiggle_0
>> chrI    0       2       169
>> chrI    2       4       176
>> chrI    4       5       178
>> chrI    5       7       179
>>
>> But it feels awkward to extract the ranges and values from the coverage
>> result, just to put them into a new kind of object. I know I'm doing
>> something wrong...any help?
>>
>> Also, according to UCSC this is a bedGraph file, not a wig file. So why
>> is the result of the wig export a bedGraph file? I don't get it.
>>
>>
> Well yes it looks like UCSC has changed the specs. WIG used to encompass
> the BED lines, and it's still supported, though undocumented. Now they call
> it a separate format "bedGraph". I will add that to rtracklayer.
>
>
As of rtracklayer 1.5.18, WIG is never exported BED-style. To do that, you
need to use the newly supported bedGraph format.

Btw, your tracks should now be transparently output as variableStep, which
> should be more compact.
>
> Michael
>
>
>> Thanks for the help.
>>
>> -Chris
>>
>>
>

	[[alternative HTML version deleted]]

