error with pileup function in ShortRead
3
0
Entering edit mode
@davide-cittaro-3332
Last seen 9.6 years ago
Hi, I'm trying to build a pileup vector from eland results. p<-pileup(position(foo), 185, chromLen[chrom], svector, width(foo)) where foo is my AlignData class (filtered on chrom). It happens that I have a read on the + strand near the end of the chromosome (156 bp afar) that is greater than the specified fragment size (185). pileup complains that the chromosome length is too small... I can try to run with a shorter fragment size, but is there a way to tell pileup that the rightmost read can be smaller than the specified average size? thanks d Davide Cittaro davide.cittaro at ifom-ieo-campus.it
• 1.0k views
ADD COMMENT
0
Entering edit mode
@herve-pages-1542
Last seen 1 day ago
Seattle, WA, United States
Hi Davide, You can use coverage() from the IRanges package to achieve this. Simple example where all reads are supposed to be on the + strand and extend to the right starting at 'start': start <- c(8:20, 76) fraglength <- 15 chrlength <- 99 pileup(start, fraglength, chrlength) as.integer(coverage(IRanges(start=start, width=fraglength), 1, chrlength)) Note that coverage() returns an Rle object which is a compact representation of the returned vector of integers. as.integer() coerce this back to a standard integer vector so you get exactly the same as with pileup(). However, unlike coverage(), pileup() knows how to pileup reads that are described by a start, a fraglength and a direction (+ or -) to which the fragment extends (in case of -, the end of the read is supposed to be at start + readlength - 1): set.seed(23) dir <- strand(sample(c("+", "-"), length(start), replace=TRUE)) readlength <- 10 pileup(start, fraglength, chrlength, dir=dir, readlength=readlength) Here is the IRanges + coverage() solution: a) Make the IRanges object describing how the reads map the reference sequence: x_start <- start ii <- dir == "-" # which ranges need to be shifted ii_readlength <- if (length(readlength) > 1) readlength[ii] else readlength ii_fraglength <- if (length(fraglength) > 1) fraglength[ii] else fraglength x_start[ii] <- x_start[ii] + ii_readlength - ii_fraglength + 1L x <- IRanges(start=x_start, width=fraglength) b) Call coverage() on it: as.integer(coverage(x, 1, chrlength)) This will not complain if chrlength looks too small: as.integer(coverage(x, 1, 20)) Hope this helps, H. Davide Cittaro wrote: > Hi, I'm trying to build a pileup vector from eland results. > > p<-pileup(position(foo), 185, chromLen[chrom], svector, width(foo)) > > where foo is my AlignData class (filtered on chrom). It happens that I > have a read on the + strand near the end of the chromosome (156 bp afar) > that is greater than the specified fragment size (185). pileup complains > that the chromosome length is too small... I can try to run with a > shorter fragment size, but is there a way to tell pileup that the > rightmost read can be smaller than the specified average size? > > thanks > > d > > Davide Cittaro > davide.cittaro at ifom-ieo-campus.it > > _______________________________________________ > Bioconductor mailing list > Bioconductor at stat.math.ethz.ch > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: > http://news.gmane.org/gmane.science.biology.informatics.conductor -- Hervé Pagès Program in Computational Biology Division of Public Health Sciences Fred Hutchinson Cancer Research Center 1100 Fairview Ave. N, M2-B876 P.O. Box 19024 Seattle, WA 98109-1024 E-mail: hpages at fhcrc.org Phone: (206) 667-5791 Fax: (206) 667-1319
ADD COMMENT
0
Entering edit mode
Simon Anders ▴ 150
@simon-anders-2626
Last seen 9.6 years ago
Dear Davide Davide Cittaro wrote: > p<-pileup(position(foo), 185, chromLen[chrom], svector, width(foo)) > > where foo is my AlignData class (filtered on chrom). It happens that I > have a read on the + strand near the end of the chromosome (156 bp afar) > that is greater than the specified fragment size (185). pileup complains > that the chromosome length is too small... I can try to run with a > shorter fragment size, but is there a way to tell pileup that the > rightmost read can be smaller than the specified average size? Instead of the scalar '185', you can pass a vector that gives an explicit fragment length estimate for each read, e.g.: p <- pileup( position(foo), ifelse( position(foo) < chromLen[chrom] - 185, 185, chromLen[chrom] - position(foo) - 1 ), chromLen[chrom], svector, width(foo) ) I don't want to claim that this is a very elegant solution. Best regards Simon +--- | Dr. Simon Anders, Dipl. Phys. | European Bioinformatics Institute, Hinxton, Cambridgeshire, UK | office phone +44-1223-492680, mobile phone +44-7505-841692 | preferred (permanent) e-mail: sanders at fs.tum.de
ADD COMMENT
0
Entering edit mode
@martin-morgan-1513
Last seen 6 weeks ago
United States
Hi Davide -- a third, different, response! Since you mention ELAND and AlignedRead, it's likely that you can accomplish your task with cvg <- coverage(foo, 1, chromLen[chrom], extend=fraglength- readlength) this also loops over chromosomes. (There was an off-by-one error, fixed in version 1.1.45, that added an extra nucleotide to reads; this has not quite reached the biocLite repository but is in svn). The relevant help page is ?"AlignedRead-class"; I'll add a section to the vignette on use of coverage in the not too distant future. Martin Hervé Pagès <hpages at="" fhcrc.org=""> writes: > Hi Davide, > > You can use coverage() from the IRanges package to achieve this. > > Simple example where all reads are supposed to be on the + strand and > extend to the right starting at 'start': > > start <- c(8:20, 76) > fraglength <- 15 > chrlength <- 99 > > pileup(start, fraglength, chrlength) > > as.integer(coverage(IRanges(start=start, width=fraglength), 1, chrlength)) > > Note that coverage() returns an Rle object which is a compact representation > of the returned vector of integers. as.integer() coerce this back to a > standard integer vector so you get exactly the same as with pileup(). > > However, unlike coverage(), pileup() knows how to pileup reads that are > described by a start, a fraglength and a direction (+ or -) to which > the fragment extends (in case of -, the end of the read is supposed to > be at start + readlength - 1): > > set.seed(23) > dir <- strand(sample(c("+", "-"), length(start), replace=TRUE)) > readlength <- 10 > > pileup(start, fraglength, chrlength, dir=dir, readlength=readlength) > > Here is the IRanges + coverage() solution: > > a) Make the IRanges object describing how the reads map the > reference sequence: > > x_start <- start > ii <- dir == "-" # which ranges need to be shifted > ii_readlength <- if (length(readlength) > 1) readlength[ii] else readlength > ii_fraglength <- if (length(fraglength) > 1) fraglength[ii] else fraglength > x_start[ii] <- x_start[ii] + ii_readlength - ii_fraglength + 1L > x <- IRanges(start=x_start, width=fraglength) > > b) Call coverage() on it: > > as.integer(coverage(x, 1, chrlength)) > > This will not complain if chrlength looks too small: > > as.integer(coverage(x, 1, 20)) > > Hope this helps, > > H. > > > Davide Cittaro wrote: >> Hi, I'm trying to build a pileup vector from eland results. >> p<-pileup(position(foo), 185, chromLen[chrom], svector, width(foo)) >> where foo is my AlignData class (filtered on chrom). It happens that >> I have a read on the + strand near the end of the chromosome (156 bp >> afar) that is greater than the specified fragment size (185). pileup >> complains that the chromosome length is too small... I can try to >> run with a shorter fragment size, but is there a way to tell pileup >> that the rightmost read can be smaller than the specified average >> size? >> thanks >> d >> Davide Cittaro >> davide.cittaro at ifom-ieo-campus.it >> _______________________________________________ >> Bioconductor mailing list >> Bioconductor at stat.math.ethz.ch >> https://stat.ethz.ch/mailman/listinfo/bioconductor >> Search the archives: >> http://news.gmane.org/gmane.science.biology.informatics.conductor -- Martin Morgan Computational Biology / Fred Hutchinson Cancer Research Center 1100 Fairview Ave. N. PO Box 19024 Seattle, WA 98109 Location: Arnold Building M2 B169 Phone: (206) 667-2793
ADD COMMENT

Login before adding your answer.

Traffic: 891 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6