Patch for the "c" issue submitted:
https://bugs.r-project.org/bugzilla3/show_bug.cgi?id=15681
On Sat, Feb 22, 2014 at 2:56 PM, Michael Lawrence <micha...@gene.com> wrote:
>
>
>
> On Sat, Feb 22, 2014 at 2:04 PM, Hervé Pagès <hpa...@fhcrc.org> wrote:
>
>>
>>
>> On 02/22/2014 01:59 PM, Michael Lawrence wrote:
>>
>>> Yea, I can see the argument for consistency. Not sure I'll get around to
>>> it though. How about we see if they take this patch first.
>>>
>>
>> Sounds good. I was hoping that fixing selectMethod would maybe
>> fix ?generic(arg) automatically but that's probably too naive.
>>
>> In the meantime I tried the patch and got:
>>
>> > library(IRanges)
>> > ?coverage(IRanges())
>> Error in args$...[[1L]] : object of type 'symbol' is not subsettable
>>
>> > library(rtracklayer)
>> > path_to_gff <- system.file("tests", "v1.gff", package = "rtracklayer")
>> > ?import(path_to_gff)
>> Error in args$...[[1L]] : object of type 'symbol' is not subsettable
>>
>> which seems to go away if I add the following line
>>
>> args <- args[fdef@signature]
>>
>> right before
>>
>> args$... <- args$...[[1L]]
>>
>> like you did in your original helpwith() proposal.
>>
>>
> Thanks for the catch. I've adapted the patch.
>
>
>> Also I found a situation that is still failing (with or without the
>> patch):
>>
>> > ?c(IRanges(),IRanges())
>> Error in .helpForCall(topicExpr, parent.frame()) :
>> no documentation for function 'c' and signature 'x = "IRanges",
>> recursive = "logical"'
>> In addition: Warning message:
>> In .helpForCall(topicExpr, parent.frame()) :
>> no method defined for function 'c' and signature 'x = "IRanges",
>> recursive = "logical"'
>>
>>
>
> Yea, we know the "c" generic is messed up. There's a bunch of variadic
> primitive generics like all, any, sum, etc, that report extra arguments in
> their signature, even though they only dispatch on the "x". I'll see if I
> can submit a patch on this.
>
>
>> Thanks again for working on this.
>>
>> H.
>>
>>
>>
>>> Michael
>>>
>>>
>>> On Sat, Feb 22, 2014 at 12:26 PM, Hervé Pagès <hpa...@fhcrc.org
>>> <mailto:hpa...@fhcrc.org>> wrote:
>>>
>>> Thanks Michael.
>>>
>>> Do you think it would be sensible to offer a similar fix for
>>> selectMethod?
>>>
>>> > setGeneric("f", function(x, y) standardGeneric("f"))
>>> > setMethod("f", c("numeric", "missing"), function(x, y) x)
>>> > selectMethod("f", "numeric")
>>> Error in selectMethod("f", "numeric") :
>>> no method found for signature numeric, ANY
>>>
>>> In a perfect world,
>>>
>>> generic(arg)
>>> selectMethod(generic, class(arg))
>>> ?generic(arg)
>>>
>>> should behave consistently.
>>>
>>> The ?import(path_to_gff) error I reported earlier was actually
>>> inspired by a sad experience I had in the past where I was trying
>>> to use selectMethod("import", ...) to find out what particular
>>> method was being called (I needed to see the code).
>>>
>>> Thanks,
>>> H.
>>>
>>>
>>> On 02/22/2014 05:59 AM, Michael Lawrence wrote:
>>>
>>> Translated it into a patch against R, submitted here:
>>> https://bugs.r-project.org/__bugzilla3/show_bug.cgi?id=__15680
>>> <https://bugs.r-project.org/bugzilla3/show_bug.cgi?id=15680>
>>>
>>>
>>> On Fri, Feb 21, 2014 at 2:53 PM, Hervé Pagès <hpa...@fhcrc.org
>>> <mailto:hpa...@fhcrc.org>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>>> wrote:
>>>
>>>
>>>
>>> On 02/21/2014 02:01 PM, Michael Lawrence wrote:
>>>
>>> This function seems to solve the problem:
>>>
>>> helpwith <- function(expr) {
>>> env <- IRanges:::top_prenv(expr)
>>> expr <- substitute(expr)
>>> fun <- eval(expr[[1L]], env)
>>> fun.name <http://fun.name> <http://fun.name> <http://fun.name>
>>> <- deparse(expr[[1L]])
>>> if (!isGeneric(fun.name <http://fun.name>
>>> <http://fun.name> <http://fun.name>,
>>> env)) {
>>>
>>> stop("'expr' must be a call to a generic")
>>> }
>>> args <- formals(fun)
>>> mc <- match.call(fun, expr, expand.dots=FALSE)
>>> args[names(mc[-1L])] <- mc[-1L]
>>> args <- args[fun@signature]
>>> args$... <- args$...[[1L]]
>>> target.sig <- vapply(args, function(arg) {
>>> .arg <- arg # nice trick
>>> if (missing(.arg)) {
>>> "missing"
>>> } else {
>>> class(eval(arg, env))[1L]
>>> }
>>> }, character(1))
>>> defined.sig <- selectMethod(fun, target.sig)@defined
>>> help(paste0(fun.name <http://fun.name>
>>> <http://fun.name> <http://fun.name>,
>>> ",", paste(defined.sig,
>>>
>>> collapse=","), "-method"))
>>> }
>>>
>>> path_to_gff <- "my.gff"
>>> helpwith(import(path_to_gff))
>>>
>>> helpwith(rbind(DataFrame(____mtcars),
>>> DataFrame(mtcars)))
>>>
>>> But where should it go?
>>>
>>>
>>> Nice fix. It's a bug in ? so it would need to go into ?
>>> itself.
>>>
>>> The man page for ? (accessible with ?`?`) says:
>>>
>>> S4 Method Documentation:
>>>
>>> [...]
>>>
>>> There are two different ways to look for
>>> documentation on a
>>> particular method. The first is to supply the
>>> 'topic' argument in
>>> the form of a function call, omitting the 'type'
>>> argument. The
>>> effect is to look for documentation on the method
>>> that would be
>>> used if this function call were actually evaluated.
>>> See the
>>> examples below. If the function is not a generic (no
>>> S4 methods
>>> are defined for it), the help reverts to
>>> documentation on the
>>> function name.
>>>
>>> Thanks,
>>> H.
>>>
>>>
>>>
>>>
>>>
>>> On Fri, Feb 21, 2014 at 11:17 AM, Hervé Pagès
>>> <hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>>>> wrote:
>>>
>>> Hi Gabriel,
>>>
>>>
>>> On 02/20/2014 05:03 PM, Gabriel Becker wrote:
>>>
>>> Herve,
>>>
>>> The help is correct (though possibly a bit
>>> pedantic),
>>> there is no
>>> method for that signature.
>>>
>>>
>>> But the dispatch mechanism is able to find one
>>> because
>>> 'import(path_to_gff)' *does* work. So the help
>>> maybe is correct
>>> but that doesn't really help the naive user.
>>>
>>> H.
>>>
>>>
>>> ?import("", "")
>>>
>>> works for me though
>>>
>>> ~G
>>>
>>>
>>> On Thu, Feb 20, 2014 at 4:51 PM, Hervé Pagès
>>> <hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>>
>>> <mailto:hpa...@fhcrc.org
>>> <mailto:hpa...@fhcrc.org> <mailto:hpa...@fhcrc.org
>>> <mailto:hpa...@fhcrc.org>>>
>>> <mailto:hpa...@fhcrc.org
>>> <mailto:hpa...@fhcrc.org> <mailto:hpa...@fhcrc.org
>>> <mailto:hpa...@fhcrc.org>>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>>>>> wrote:
>>>
>>> On 02/20/2014 04:16 PM, Michael Lawrence
>>> wrote:
>>>
>>> There's also "?coverage(ga)", so the
>>> user can
>>> see what
>>> happens
>>> specifically for their object, without
>>> worrying about
>>> typing out
>>> the class.
>>>
>>>
>>> I was never lucky with this:
>>>
>>> > library(rtracklayer)
>>> > path_to_gff <- system.file("tests",
>>> "v1.gff",
>>> package =
>>> "rtracklayer")
>>> > ?import(path_to_gff)
>>> Error in .helpForCall(topicExpr,
>>> parent.frame()) :
>>> no documentation for function
>>> 'import' and
>>> signature
>>> 'con =
>>> "character", format = "ANY", text = "ANY"'
>>> In addition: Warning message:
>>> In .helpForCall(topicExpr,
>>> parent.frame()) :
>>> no method defined for function
>>> 'import' and
>>> signature
>>> 'con =
>>> "character", format = "ANY", text = "ANY"'
>>>
>>> H.
>>>
>>>
>>> At some point it would be neat to
>>> generate S4
>>> documentation at
>>> run-time.
>>> Just popup a browser and see every
>>> method that
>>> might be
>>> dispatched with
>>> a given object. In theory, the R help
>>> server would
>>> support this.
>>>
>>>
>>> On Thu, Feb 20, 2014 at 3:13 PM,
>>> Hervé Pagès
>>> <hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>>>
>>> <mailto:hpa...@fhcrc.org
>>> <mailto:hpa...@fhcrc.org>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>>>>
>>> <mailto:hpa...@fhcrc.org
>>> <mailto:hpa...@fhcrc.org>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>>>
>>> <mailto:hpa...@fhcrc.org
>>> <mailto:hpa...@fhcrc.org> <mailto:hpa...@fhcrc.org
>>> <mailto:hpa...@fhcrc.org>>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>>>>>> wrote:
>>>
>>>
>>>
>>> On 02/20/2014 02:55 PM, Martin
>>> Morgan wrote:
>>>
>>> On 02/20/2014 02:32 PM,
>>> Hervé Pagès
>>> wrote:
>>>
>>> Hi Jesper,
>>>
>>> On 02/20/2014 02:13 PM,
>>> Jesper
>>> Gådin wrote:
>>>
>>> Very true that it is
>>> quite
>>> difficult
>>> to find the
>>> documentation when
>>> one
>>> is not aware of its
>>> existence :P
>>>
>>>
>>> Yeah, this has been a
>>> source of
>>> frustration for
>>> many people. And
>>> always a source of
>>> embarrassment
>>> (for us) when
>>> teaching our
>>> software
>>> to beginners.
>>>
>>> I've started to change
>>> this. In
>>> the upcoming
>>> version of BioC
>>> (2.14,
>>> scheduled for
>>> mid-April), when
>>> you'll do
>>> ?coverage,
>>> you'll
>>> get to
>>> choose between the 3 man
>>> pages that
>>> document coverage
>>> methods (there
>>> is one in IRanges, one in
>>> GenomicRanges,
>>> and one in
>>> GenomicAlignments).
>>>
>>> I want to generalize
>>> this to other
>>> generics that have
>>> methods spread
>>> across several packages
>>> (e.g.
>>> findOverlaps, the
>>> intra- and
>>> inter-range
>>> methods, etc...).
>>>
>>> Having to choose between
>>> several
>>> man pages
>>> every
>>> time you do
>>> e.g.
>>> ?findOverlaps is a minor
>>> annoyance
>>> compared to not
>>> being able to
>>> find the man page at
>>> all. (Note
>>> that if
>>> you already
>>> know
>>> where is
>>> your favorite man page,
>>> you'll be
>>> able to
>>> direct
>>> access it with
>>> e.g.
>>> ?GenomicRanges::findOverlaps). Nobody
>>> will
>>> ever need to use
>>> the insane
>>>
>>>
>>>
>>> ?`findOverlaps,GenomicRanges,__________GenomicRanges-method` to
>>>
>>>
>>>
>>>
>>>
>>> tab completion helps, so
>>> that you
>>> don't need
>>> to be totally
>>> insane, just
>>> insane enough to know to
>>> start with
>>>
>>> ?"cover<tab>
>>>
>>>
>>> and also insane enough to know
>>> which
>>> method you
>>> need to
>>> pick up amongst
>>> the 30+ methods listed by
>>> ?"findOverlaps<tab>
>>> Overwhelming!
>>>
>>> H.
>>>
>>>
>>>
>>> Martin
>>>
>>> open that man page
>>> again. Ever!
>>> (it will
>>> still work
>>> though...)
>>>
>>> Cheers,
>>> H.
>>>
>>>
>>> coverage() is fast
>>> and
>>> beautiful. Thanks!
>>>
>>> /Jesper
>>>
>>>
>>> On Wed, Feb 19, 2014
>>> at 9:21
>>> PM, Hervé
>>> Pagès
>>> <hpa...@fhcrc.org
>>> <mailto:hpa...@fhcrc.org>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>>
>>> <mailto:hpa...@fhcrc.org
>>> <mailto:hpa...@fhcrc.org> <mailto:hpa...@fhcrc.org
>>> <mailto:hpa...@fhcrc.org>>>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>>
>>> <mailto:hpa...@fhcrc.org
>>> <mailto:hpa...@fhcrc.org> <mailto:hpa...@fhcrc.org
>>> <mailto:hpa...@fhcrc.org>>>>
>>> <mailto:hpa...@fhcrc.org
>>> <mailto:hpa...@fhcrc.org>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>>>
>>> <mailto:hpa...@fhcrc.org
>>> <mailto:hpa...@fhcrc.org> <mailto:hpa...@fhcrc.org
>>> <mailto:hpa...@fhcrc.org>>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>>>>>
>>>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>>
>>> <mailto:hpa...@fhcrc.org
>>> <mailto:hpa...@fhcrc.org> <mailto:hpa...@fhcrc.org
>>> <mailto:hpa...@fhcrc.org>>>
>>> <mailto:hpa...@fhcrc.org
>>> <mailto:hpa...@fhcrc.org>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>>>>
>>> <mailto:hpa...@fhcrc.org
>>> <mailto:hpa...@fhcrc.org> <mailto:hpa...@fhcrc.org
>>> <mailto:hpa...@fhcrc.org>>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>>>
>>> <mailto:hpa...@fhcrc.org
>>> <mailto:hpa...@fhcrc.org>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>
>>> <mailto:hpa...@fhcrc.org <mailto:hpa...@fhcrc.org>>>>>>
>>> >
>>>
>>> wrote:
>>>
>>> Hi Jesper,
>>>
>>>
>>> On 02/19/2014
>>> 08:44 AM,
>>> Michael
>>> Lawrence
>>> wrote:
>>>
>>> On Wed, Feb
>>> 19, 2014
>>> at 8:39 AM,
>>> Jesper Gådin
>>>
>>> <jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>__>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>__>__>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>__>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>__>__>__>
>>>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>__>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>__>__>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>__>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>__>__>__>__>
>>>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>__>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>__>__>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>__>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>__>__>__>
>>>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>__>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>__>__>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>__>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>
>>> <mailto:jesper.ga...@gmail.com
>>> <mailto:jesper.ga...@gmail.com>__>__>__>__>__>>wrote:
>>>
>>>
>>>
>>>
>>> Hi
>>> Michael,
>>>
>>> Herves
>>> suggestion will
>>> probably
>>> work for my
>>> use case, but if
>>> there
>>> are any
>>>
>>> smoother ways it
>>> would be
>>> preferable.
>>>
>>> The use
>>> case is
>>> as follows:
>>>
>>> 1)
>>> calculate
>>> strand specific
>>> coverage over
>>> a region from
>>>
>>> GAlignments
>>> object (or file)
>>>
>>> At the
>>> moment I
>>> read a
>>> file using
>>>
>>> readGAlignmentsFromBam()
>>> with
>>> tag="XS",
>>> then
>>> filter it
>>> on "flag" and
>>> "mapq". Then I
>>> subset the
>>>
>>> resulting
>>> GAlignments in
>>> a minus
>>> and a plus
>>> -strand object.
>>> Then I
>>> calculate
>>> coverage
>>> by my own
>>> coverage function
>>> which
>>> uses
>>> the cigar
>>>
>>> information in the
>>> GAlignments
>>> object. This
>>> function is the
>>> one
>>> using
>>>
>>> cigarToRleList()
>>> at the
>>> moment. As I
>>> understand the
>>>
>>> coverage() function
>>> from the
>>> GenomicAlignments/IRanges
>>> package
>>> does not take
>>> into
>>> account
>>> cigars,
>>> or does it?
>>>
>>>
>>> It does
>>> take the
>>> coverage
>>> into account;
>>> specifically to
>>> exclude
>>> the introns
>>> from
>>> coverage. I think
>>> there's also an
>>> option
>>> to exclude
>>> deletions.
>>>
>>>
>>> Unfortunately
>>> the man
>>> page is not
>>> easy to
>>> access
>>> (you need to do
>>>
>>> ?`coverage,GAlignments-method`____________),
>>> but
>>>
>>>
>>> it says:
>>>
>>> The methods
>>> for
>>> GAlignments and
>>> GAlignmentPairs
>>> objects do:
>>>
>>>
>>> coverage(grglist(x), ...)
>>>
>>> And if you do
>>> grglist() on a
>>> GAlignments or
>>> GAlignmentPairs
>>> objects, the
>>> ranges you get
>>> in the
>>> returned
>>> GRangesList
>>> object
>>> are calculated
>>> based
>>> on the CIGAR.
>>>
>>> Trust but
>>> verify. Here
>>> is how you
>>> can actually
>>> verify that
>>> coverage()
>>> does take the
>>> CIGAR into
>>> account:
>>>
>>>
>>> library(RNAseqData.HNRNPC.bam.
>>> ____________chr14)
>>> gal <-
>>>
>>>
>>>
>>>
>>> readGAlignments(RNAseqData.____________HNRNPC.bam.chr14_____
>>> BAMFILES[__1])
>>>
>>>
>>>
>>> cig_op_table
>>> <-
>>> cigarOpTable(cigar(gal))
>>>
>>> First we pick
>>> up an
>>> alignment
>>> with an N in its
>>> CIGAR and do
>>> coverage()
>>> on it:
>>>
>>> > gal_with_N
>>> <-
>>> gal[which(cig_op_table[
>>> , "N"]
>>> != 0)[1]]
>>>
>>> > gal_with_N
>>> GAlignments
>>> with 1
>>> alignment and 0
>>> metadata columns:
>>>
>>> seqnames strand
>>> cigar qwidth
>>> start end
>>> width
>>>
>>> <Rle> <Rle>
>>> <character> <integer>
>>> <integer> <integer>
>>> <integer>
>>> [1]
>>> chr14 +
>>> 55M2117N17M 72
>>> 19650072 19652260
>>> 2189
>>>
>>> ngap
>>>
>>> <integer>
>>> [1]
>>> 1
>>> ---
>>> seqlengths:
>>>
>>> chr1
>>> chr10 ...
>>> chrY
>>>
>>> 249250621
>>> 135534747 ...
>>> 59373566
>>>
>>> >
>>> coverage(gal_with_N)$chr14
>>> integer-Rle
>>> of length
>>> 107349540 with 5 runs
>>> Lengths:
>>> 19650071
>>> 55
>>> 2117
>>> 17
>>> 87697280
>>> Values :
>>> 0
>>> 1
>>> 0
>>> 1
>>> 0
>>>
>>> Same thing with
>>> an
>>> alignment with
>>> an I in
>>> its CIGAR:
>>>
>>> > gal_with_I
>>> <-
>>> gal[which(cig_op_table[
>>> , "I"]
>>> != 0)[1]]
>>>
>>> > gal_with_I
>>> GAlignments
>>> with 1
>>> alignment and 0
>>> metadata columns:
>>>
>>> seqnames strand
>>> cigar qwidth
>>> start end
>>> width
>>>
>>> <Rle> <Rle>
>>> <character> <integer>
>>> <integer> <integer>
>>> <integer>
>>> [1]
>>> chr14 -
>>> 64M1I7M 72
>>> 19411677 19411747
>>> 71
>>>
>>> ngap
>>>
>>> <integer>
>>> [1]
>>> 0
>>> ---
>>> seqlengths:
>>>
>>> chr1
>>> chr10 ...
>>> chrY
>>>
>>> 249250621
>>> 135534747 ...
>>> 59373566
>>>
>>> >
>>> coverage(gal_with_I)$chr14
>>> integer-Rle
>>> of length
>>> 107349540 with 3 runs
>>> Lengths:
>>> 19411676
>>> 71 87937793 <tel:71%2087937793>
>>> <tel:71%2087937793>
>>> Values :
>>> 0
>>> 1
>>> 0
>>>
>>> Same thing with
>>> an
>>> alignment with
>>> a D in
>>> its CIGAR:
>>>
>>> > gal_with_D
>>> <-
>>> gal[which(cig_op_table[
>>> , "D"]
>>> != 0)[1]]
>>>
>>> > gal_with_D
>>> GAlignments
>>> with 1
>>> alignment and 0
>>> metadata columns:
>>>
>>> seqnames strand
>>> cigar qwidth
>>> start end
>>> width
>>>
>>> <Rle> <Rle>
>>> <character> <integer>
>>> <integer> <integer>
>>> <integer>
>>> [1]
>>> chr14 +
>>> 38M1D34M 72
>>> 19659063 19659135
>>> 73
>>>
>>> ngap
>>>
>>> <integer>
>>> [1]
>>> 0
>>> ---
>>> seqlengths:
>>>
>>> chr1
>>> chr10 ...
>>> chrY
>>>
>>> 249250621
>>> 135534747 ...
>>> 59373566
>>>
>>> >
>>> coverage(gal_with_D)$chr14
>>> integer-Rle
>>> of length
>>> 107349540 with 3 runs
>>> Lengths:
>>> 19659062
>>> 73
>>> 87690405
>>> Values :
>>> 0
>>> 1
>>> 0
>>>
>>> Seeing is
>>> believing,
>>>
>>> Cheers,
>>> H.
>>>
>>>
>>>
>>> I
>>> started to
>>> look at the
>>> applyPileups() in
>>> Rsamtools which I
>>> can get
>>> to
>>>
>>> calculate
>>> coverage using
>>> cigars,
>>> but not
>>> using the strand or
>>> flag
>>>
>>> information for
>>> filtering. That
>>> solution
>>> would start from a
>>> bam-file
>>> instead
>>> of a
>>> GAlignments
>>> object,
>>> and sure I
>>> can do the
>>>
>>> filtering outside R.
>>> But it
>>> would be very
>>> convenient to
>>> do it
>>> all from within R.
>>>
>>> If
>>> there are
>>> nice solutions
>>> starting from
>>> both a GAlignments
>>
>>
>
[[alternative HTML version deleted]]
_______________________________________________
Bioc-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/bioc-devel
