[Bioc-sig-seq] WIG from coverage
Seidel, Christopher
CWS at stowers.org
Wed Sep 30 01:28:24 CEST 2009
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?). Shouldn't rtracklayer be putting out a variableStep format for this? Am I doing something wrong?
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
More information about the Bioc-sig-sequencing
mailing list