On Wed, Sep 30, 2009 at 11:20 AM, Chris Seidel <[email protected]> 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: [email protected] [mailto:[email protected]] On Behalf Of Michael > Lawrence > Sent: Wednesday, September 30, 2009 1:22 AM > To: Seidel, Christopher > Cc: [email protected] > Subject: Re: [Bioc-sig-seq] WIG from coverage > > On Tue, Sep 29, 2009 at 7:23 PM, Michael Lawrence <[email protected]> > 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. 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]] _______________________________________________ Bioc-sig-sequencing mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing
