On Thu, Sep 17, 2009 at 4:35 PM, Patrick Aboyoun <[email protected]
<mailto:[email protected]>> wrote:
Michael,
We could add that level of flexibility, but I worry that may open
the door too widely because the aggregator will need to take the
old RangedData object and the new block RangedData object as
arguments and produce a new RangedData object for the next
iteration and users may find it difficult to construct a valid
aggregator.
Well just to clarify, I thought the aggregator would take the old and
new objects as output by the block processor. So for coverage, one
function would take a RangedData and return a RleList, and the
aggregator would take two RleLists and add them.
The call for calculating the coverage would be like:
cov <- import("reads.bed", blockSize = 500000, blockFun = coverage,
aggregateFun = `+`)
Of course, high-level functionality would be nice for routine
operations. That will help prevent the user from making mistakes. The
blockFun and aggregateFun could be bound into a single class
BlockHandler.
cov <- import("reads.bed", blockSize = 500000, blockHandler =
CoverageHandler())
I think there may also be a need for filters to allow the user to
keep only those rows that satisfy their criteria. We could supply
chromosome, strand, score, etc. filters to meet this need.
Yes, and the Filter class could inherit from BlockHandler, so that it
works incrementally. import() could be smart enough to incorporate
both a filter and a user block handler into a single "meta" handler.
So to obtain an RleList of coverage, but only over chromosome 1,
processing half a million reads at a time:
cov <- import("reads.bed", blockSize = 500000, blockHandler =
CoverageHandler(), filter = ChromosomeFilter("chr1"))
Michael
Patrick
Michael Lawrence wrote:
I don't know if we want to add explicit support for scoring,
when there are many reasons why one might want to
incrementally load a track. One solution would be a design
based on callbacks. The user specifies the block size and a
callback to be invoked for each block. The result could be a
list, with each element holding the result for a block. The
user could then aggregate that. To be even more efficient, we
could support a second callback that aggregates at each step.
On Thu, Sep 17, 2009 at 3:47 PM, Patrick Aboyoun
<[email protected] <mailto:[email protected]>
<mailto:[email protected] <mailto:[email protected]>>> wrote:
Michael,
Before people start developing their own bed file parsers, it
looks pretty straight-forward to incorporate this operation in
rtracklayer now. You or I could add a scoreRanges argument,
which
takes either a logical or perhaps a character string if we
envision there being more than one way to score a set of
ranges,
to the import* functions and when a appropriate value is
supplied,
it will internally block process the data file and output a
RangedData object with a score column.
Patrick
Michael Lawrence wrote:
On Thu, Sep 17, 2009 at 2:17 PM, Ivan Gregoretti
<[email protected] <mailto:[email protected]>
<mailto:[email protected] <mailto:[email protected]>>> wrote:
Hi Michael,
True; there are no scores in this BED file. That
brings us
to Sean's
remark.
You may have noticed that all features in the BED
file are
36 bases
long. Each feature is the position and orientation of a
Solexa read.
I want the value in each bin in the WIG file to
represent
the number
of tags contained in that bin.
So, coverage() is close to what I am looking for but I
think that
pileup() is perhaps better. Then comes the problem
of how
to place
those values in an object exportable as WIG. Lets
not talk
about
efficiency for now.
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")
I think that I can solve this problem myself but
given that:
1) there is a package that can read BED and write
WIG files
and
2) BED to WIG conversion is such a common task
I thought that there must have been an existing tool or
perhaps an
established efficient way of doing this with BioC
tools.
Why is this
important? Well, this myRegions.bed contains 4.3
million
records.
Loading them with import() used about 1GB of RAM.
My real
world BED
files are never less than 50 million reads, 150
more like.
Conclusion: I DO care about an expert's opinion of
BED-> WIG
conversion within R.
Maybe nobody does this with R. That is also acceptable.
What
tool do
you use then?
The bottleneck seems to be loading and representing
large BED
files, since
after the reduction to coverage the data size is usually
manageable.
rtracklayer was never designed to load BED files with
millions
of records.
In my experience, short reads are more commonly
represented in
more
efficient formats like those output by MAQ and bowtie,
which
are imported
with ShortRead.
There is much room for optimization in rtracklayer, but for
now one idea is
to calculate coverage on subsets of the reads, and
aggregate
the results.
The BED file can be loaded with simple low-level calls like
read.table and
scan.
Michael
Thank you,
Ivan
Ivan Gregoretti, PhD
National Institute of Diabetes and Digestive and Kidney
Diseases
National Institutes of Health
5 Memorial Dr, Building 5, Room 205.
Bethesda, MD 20892. USA.
Phone: 1-301-496-1592
Fax: 1-301-496-9878
On Thu, Sep 17, 2009 at 4:45 PM, Sean Davis
<[email protected] <mailto:[email protected]>
<mailto:[email protected] <mailto:[email protected]>>> wrote:
On Thu, Sep 17, 2009 at 4:43 PM, Michael Lawrence
<[email protected] <mailto:[email protected]>
<mailto:[email protected] <mailto:[email protected]>>>
wrote:
On Thu, Sep 17, 2009 at 8:14 AM, Ivan
Gregoretti
<[email protected]
<mailto:[email protected]> <mailto:[email protected]
<mailto:[email protected]>>>
wrote:
Hello everybody
How do you convert BED formatted files
to WIG
files?
I tried to do that with rtracklayer but it
didn't quite succeed.
This is the session's transcript:
First, you can download
http://dl.getdropbox.com/u/2051155/myRegions.bed,
which looks like
this
chr1 3002444 3002479 +
chr1 3002989 3003024 -
chr1 3017603 3017638 +
chr1 3017879 3017914 -
chr1 3018173 3018208 +
chr1 3018183 3018218 -
chr1 3018183 3018218 -
chr1 3019065 3019100 +
chr1 3019761 3019796 -
chr1 3020044 3020079 -
...
Now to R
suppressMessages(library(rtracklayer))
myRegions <- import('myRegions.bed')
So far so good. Now I try:
export(myRegions, 'myRegions.wig',
format = 'wig')
but I get:
Error in export.ucsc(object, con,
subformat,
...) :
Track not compatible with WIG format:
Overlapping features must be
on separate strands and every feature width
must be positive
I seems that the error message is a feature
rather than a bug. My
interpretation is that export() does
not like
records like lines 6 and
7:
chr1 3018183 3018218 -
chr1 3018183 3018218 -
So, how do you convert BED to WIG in your
everyday work?
Well, as rtracklayer said, this track is not
representable by WIG, which
is
meant to communicate a single data value for a
given genomic region (to
generate the bar plots in the UCSC browser, for
instance). There aren't
any
scores in your file, so why do you want to
use WIG?
Just guessing, here, but you may want to calculate
coverage() first,
Ivan?
Sean
Thank you,
Ivan
sessionInfo()
R version 2.9.2 (2009-08-24)
x86_64-redhat-linux-gnu
locale:
LC_CTYPE=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8;LC_COLLATE=en_US.UTF-8;LC_MONETARY=C;LC_MESSAGES=en_US.UTF-8;LC_PAPER=en_US.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_US.UTF-8;LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils
datasets methods base
other attached packages:
[1] rtracklayer_1.4.0 RCurl_0.94-1
loaded via a namespace (and not attached):
[1] Biobase_2.4.0 Biostrings_2.12.0
BSgenome_1.12.0
IRanges_1.2.0
[5] tools_2.9.2 XML_2.3-0
Ivan Gregoretti, PhD
National Institute of Diabetes and
Digestive
and Kidney Diseases
National Institutes of Health
5 Memorial Dr, Building 5, Room 205.
Bethesda, MD 20892. USA.
Phone: 1-301-496-1592
Fax: 1-301-496-9878
_______________________________________________
Bioc-sig-sequencing mailing list
[email protected]
<mailto:[email protected]>
<mailto:[email protected]
<mailto:[email protected]>>
https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing
[[alternative HTML version deleted]]
_______________________________________________
Bioc-sig-sequencing mailing list
[email protected]
<mailto:[email protected]>
<mailto:[email protected]
<mailto:[email protected]>>
https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing
_______________________________________________
Bioc-sig-sequencing mailing list
[email protected]
<mailto:[email protected]>
<mailto:[email protected]
<mailto:[email protected]>>
https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing
[[alternative HTML version deleted]]
_______________________________________________
Bioc-sig-sequencing mailing list
[email protected]
<mailto:[email protected]>
<mailto:[email protected]
<mailto:[email protected]>>
https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing
