On Tue, Sep 29, 2009 at 7:23 PM, Michael Lawrence <[email protected]>wrote:
> > > On Tue, Sep 29, 2009 at 4:28 PM, Seidel, Christopher <[email protected]>wrote: > >> Hello, >> >> In an earlier thread: BED to WIG format conversion (Thu 9/17/2009 4:50 PM) >> Michael Lawrence mentioned: >> >> Ignoring efficiency, it's pretty straight-forward. In the development >> version of IRanges, one can convert an Rle (or RleList for multiple >> chromosomes) as output to coverage() to a RangedData. In the development >> version of rtracklayer, this is done implicitly: >> >> > export(cov, "coverage.wig") >> >> as a way to create a wig file from a coverage() result. When I try this, >> my wig file looks like the following: >> >> track name="R Track" type=wiggle_0 >> . 0 2 168 >> . 2 4 175 >> . 4 5 177 >> . 5 7 178 >> . 7 8 182 >> . 8 12 185 >> . 12 14 188 >> . 14 15 191 >> etc. >> >> This isn't compatible with any of the UCSC wig formats that I'm aware of >> (is it?). > > > No, it should be outputting the chromosome names instead of ".", following > the BED-line WIG variant. > > The coercion from Rle to RangedData yields a single-chromosome RangedData > with NULL names. This becomes <NA> upon coercion to a UCSCData (which tries > to generate conforming names but fails). <NA> shows up as "." in BED. > > New policy: RangedData cannot have NULL names. Made a validity check for > this. Fixed Rle->RangedData coercion to yield default of "1". This mimics > the behavior of RangedData(). > > Shouldn't rtracklayer be putting out a variableStep format for this? > > > Well, come to think of it, that should be possible. Right now, rtracklayer > only uses variableStep when the span is the same across an entire > chromosome. Should be possible to split the chromosome into multiple blocks > though. > I just made this work. It generates much smaller WIG files for coverage-like tracks. Thanks for the suggestion. Michael > > >> 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. > > > >> By the way, here is my work flow for going from solexa export to coverage: >> >> # 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) >> >> # find coverage for chromosome 1 >> ip1.cov <- coverage(ip1, width=230208) >> wce1.cov <- coverage(wce1, width=230208) >> >> # make a track >> export(ip1.cov, "ip1_coverage.wig") >> >> > sessionInfo() >> R version 2.10.0 Under development (unstable) (2009-09-03 r49555) >> x86_64-unknown-linux-gnu >> >> locale: >> [1] C >> >> attached base packages: >> [1] stats graphics grDevices utils datasets methods base >> >> other attached packages: >> [1] rtracklayer_1.5.13 RCurl_1.2-0 bitops_1.0-4.1 >> chipseq_0.1.25 >> [5] ShortRead_1.3.34 lattice_0.17-25 BSgenome_1.13.11 >> Biostrings_2.13.39 >> [9] IRanges_1.3.75 >> >> loaded via a namespace (and not attached): >> [1] Biobase_2.5.6 XML_2.6-0 grid_2.10.0 hwriter_1.1 >> >> -Chris >> >> Chris Seidel, Ph.D. >> http://research.stowers-institute.org/cws >> >> _______________________________________________ >> Bioc-sig-sequencing mailing list >> [email protected] >> https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing >> > > [[alternative HTML version deleted]] _______________________________________________ Bioc-sig-sequencing mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing
