It's clear we're not going to agree. I will volunteer to maintain these,
OK? I will maintain them anyway for my own use.

Michael


On Tue, May 21, 2013 at 2:34 PM, Hervé Pagès <hpages@fhcrc.org> wrote:

>
>
> On 05/21/2013 08:00 AM, Michael Lawrence wrote:
>
>>
>>
>>
>> On Mon, May 20, 2013 at 10:36 PM, Hervé Pagès <hpages@fhcrc.org
>> <mailto:hpages@fhcrc.org>> wrote:
>>
>>     Michael,
>>
>>
>>     On 05/20/2013 08:44 PM, Michael Lawrence wrote:
>>
>>
>>
>>
>>         On Mon, May 20, 2013 at 4:13 PM, Hervé Pagès <hpages@fhcrc.org
>>         <mailto:hpages@fhcrc.org>
>>         <mailto:hpages@fhcrc.org <mailto:hpages@fhcrc.org>>> wrote:
>>
>>              Hi,
>>
>>
>>              On 05/20/2013 03:15 PM, Michael Lawrence wrote:
>>
>>                  On Mon, May 20, 2013 at 3:11 PM, Martin Morgan
>>                  <mtmorgan@fhcrc.org <mailto:mtmorgan@fhcrc.org>
>>         <mailto:mtmorgan@fhcrc.org <mailto:mtmorgan@fhcrc.org>>> wrote:
>>
>>                      Hi Michael --
>>
>>
>>                      On 5/20/2013 5:56 AM, Michael Lawrence wrote:
>>
>>                          While it's cool that seqlevels<- has become so
>>         flexible,
>>                          I still claim
>>                          that
>>                          rename/keep/drop would be a lot easier to read,
>>         because
>>                          they describe the
>>                          high-level operation, and there's no reason to
>>         deprecate
>>                          them. They're
>>                          also
>>                          easier to remember. For example, for dropping
>> with
>>                          seqlevels<-, the user
>>                          needs
>>                          to remember that "force" is necessary to drop.
>> Much
>>                          easier to just say
>>                          "dropSeqlevels(), please". And reimplementing
>>                          keepSeqlevels is still too
>>                          complicated. Is it such a maintenance burden to
>>         have
>>                          these simple
>>                          wrappers sit
>>                          on top of the low-level manipulators?
>>
>>                          Another suggestion: renameSeqlevels should not
>>         require a
>>                          named vector (in
>>                          cases
>>                          where it is unnamed, require that it is
>>         parallel to the
>>                          existing
>>                          seqlevels, as
>>                          with seqlevels<-).
>>
>>
>>                      I didn't really indicate what drove my desire to see
>>                      keepSeqlevels
>>                      deprecated. keepSeqlevels, seqlevels<-, and
>>         isActiveSeq were
>>                      developed more
>>                      or less independently, and have different
>>         contracts. The
>>                      different
>>                      contracts are confusing to the user, as is creating
>>         a usable
>>                      help system
>>                      when there are multiple end points for what is a
>> common
>>                      operation. The help
>>                      pages of each were inadequate in different ways.
>>         And because
>>                      the code was
>>                      developed independently, support for different
>>         objects was
>>                      inconsistent. So
>>                      actually, yes, the maintenance (and use) burden for
>> the
>>                      previous state of
>>                      affairs was substantial.
>>
>>                      On the other hand, I agree that keepSeqlevels is
>>         convenient
>>                      as a simple
>>                      wrapper around seqlevels<-, in the same way that
>>         setNames
>>                      and names<- are
>>                      both useful.
>>
>>                      So we could iterate to
>>
>>                          keepSeqlevels <-
>>                              function(x, value, ...)
>>                          {
>>                              seqlevels(x, force=TRUE) <- value
>>                              x
>>                          }
>>
>>                      but I'd be less enthusiastic about maintaining the
>>         original
>>                      contract of
>>                      keepSeqlevels, where keepSeqlevels(gr1, gr2), would
>>         have
>>                      worked for two
>>                      GRanges objects.
>>
>>
>>              Why would this be called keepSeqlevels() if, depending on
>>         how it's used,
>>              it will either add, drop, rename, or permute the seqlevels?
>>              Couldn't this be called setSeqlevels?
>>
>>
>>         I thought that new2old was necessary for permuting.
>>
>>
>>     The seqlevels setter has no 'new2old' arg.
>>
>>
>>         You're right that
>>         adding should be disallowed for keepSeqlevels(). Adding
>>         seqlevels is not
>>         a common operation. The two common operations are:
>>         * Permuting, usually because the data were imported in
>> non-canonical
>>         order (seqnameStyle was designed to address this, no?).
>>
>>
>>     Permuting is straightforward with seqlevels<-:
>>
>>        > seqlevels(gr)
>>        [1] "chr1" "chr2" "chr3"
>>
>>        > seqlevels(gr) <- rev(seqlevels(gr))
>>
>>        > seqlevels(gr)
>>        [1] "chr3" "chr2" "chr1"
>>
>>
>>         * Subsetting, either by keeping or dropping.
>>
>>
>>     Also straightforward:
>>
>>        > seqlevels(gr, force=TRUE) <- "chr2"
>>
>>        > seqlevels(gr)
>>        [1] "chr2"
>>
>>
>>         Two main reasons: taking a
>>         small slice of the data for prototyping, or removing problematic
>>         chromosomes (sex, circular). This goes back to bsapply and the
>>         exclude
>>         argument.
>>
>>
>>     There are other important use cases:
>>
>>        - Dropping seqlevels that are *not* in use. This happens for
>> example
>>          when subsetting a BAM file with Rsamtools::filterBam to keep only
>>          the alignments located on a given chromosome. The entire sequence
>>          dictionary (in the header) is propagated to the resulting BAM
>> file
>>          so when you read the file with readGAlignments() you end up with
>> a
>>          bunch of useless seqlevels that you generally start by dropping.
>>          You don't want to drop any alignment so force=FALSE is your
>> friend.
>>
>>
>> This is sort of tangential to the discussion, but do you really want to
>> do this? I would preserve the universe as given by the BAM.
>>
>>        - An even more common operation is renaming: 90% of the times
>> that's
>>          what I use seqlevels<- for, and that's what I tell people to use
>>          when they need to rename their chromosomes. Also straightforward:
>>
>>            > seqlevels(gr)
>>            [1] "chr1" "chr2" "chr3" "chrM"
>>
>>            > seqlevels(gr) <- sub("^chr", "", seqlevels(gr))
>>            > seqlevels(gr)
>>            [1] "1" "2" "3" "M"
>>
>>            > seqlevels(gr)[seqlevels(gr) == "M"] <- "MT"
>>            > seqlevels(gr)
>>            [1] "1"  "2"  "3"  "MT"
>>
>>          Note that this is simpler to use than renameSeqlevels() which
>>          always required you to build the entire right value as a named
>>          vector.
>>
>>
>> Sure, renaming is a common use case. Not sure how I forgot that.
>>
>> As I wrote earlier in the thread, renameSeqlevels should be changed so
>> as not to require naming of the vector.
>>
>>     So the added-value of keepSeqlevels() seems to boil down to:
>>
>>        (a) it always uses force=TRUE
>>
>>        (b) it preserves the original object and returns the modified
>> object
>>
>>     If you want to restrict the use of keepSeqlevels() to permuting and
>>     dropping, you'll also have to disallow renaming (in addition to
>> adding).
>>     Then its name will more or less reflect what it does (if "keep" means
>>     "permute" or "drop"). The final result will be something that does
>> less
>>     than setSeqlevels() and that is not necessarily easier to read: both
>>     will set on the object whatever vector of seqlevels they are supplied,
>>     except that one will fail when doing so would actually result in
>> adding
>>     or renaming seqlevels.
>>
>>
>> So it looks like seqlevels<- is pretty powerful. Now I see that if I
>> scroll down to the examples in the man page, I find the actual
>> documentation. It's cool to learn that we can replace seqlevels on TxDb
>> objects.
>>
>
> Has been working on TranscriptDb objects for more than
>
>
>  My argument has always been that since these low-level
>> operations are so powerful, it's nice to have high-level operations to
>> clarify the code. We have to think hard to know what the RHS will do in
>> that replacement.
>>
>
> When you do something like seqlevels(gr) <- new_seqlevels, you're saying
> "can you please replace the current seqlevels by those new ones".
> That's all. The request/contract is simple and all the complexity is
> taken care of internally by seqlevels<-. It might be that, depending on
> what's currently on 'gr', this will actually result in dropping, adding,
> permuting. Or *a combination of the 3*. Like here:
>
>
>   > seqlevels(gr)
>   [1] "chr1" "chr2" "chr3"
>
>   > seqlevels(gr) <- c("chr4", "chr3", "chr1")
>   > seqlevels(gr)
>   [1] "chr4" "chr3" "chr1"
>
>   # one seqlevel was drop, one was added, and the 2 seqlevels that were
>   # preserved were reordered.
>
> When you want to subset a vector to keep all elements with an even
> index you do x[c(FALSE, TRUE)]. Is that keeping or dropping elements?
> Do we need to have a subsetting operator for keeping elements, another
> one for dropping elements, and another one for permuting them (like
> in x[sample(length(x))])?
>
> H.
>
>  It's probably one of the most complex replacement
>> operations; far more complex than the typical one, including levels<- on
>> factor. There's nothing wrong with that; it's just complexity that we
>> should be able to hide.
>>
>> Michael
>>
>>     H.
>>
>>
>>         Michael
>>
>>
>>              H.
>>
>>
>>
>>                  Well, I wasn't even aware of that feature, so you've
>>         made your
>>                  point about
>>                  the documentation ;) Sounds like a good compromise to me.
>>
>>                  Thanks for understanding,
>>                  Michael
>>
>>
>>
>>                      Martin
>>
>>
>>                          Michael
>>
>>
>>
>>
>>
>>
>>                          On Sat, May 18, 2013 at 6:00 PM, Martin Morgan
>>                          <mtmorgan@fhcrc.org <mailto:mtmorgan@fhcrc.org>
>>         <mailto:mtmorgan@fhcrc.org <mailto:mtmorgan@fhcrc.org>>
>>                          <mailto:mtmorgan@fhcrc.org
>>         <mailto:mtmorgan@fhcrc.org> <mailto:mtmorgan@fhcrc.org
>>         <mailto:mtmorgan@fhcrc.org>>>>
>>
>>                          wrote:
>>
>>                                On 05/18/2013 05:42 PM, Martin Morgan
>> wrote:
>>
>>                                    Some of the most common operations are
>>                          straight-forward, e.g.,
>>
>>                                          > gr = GRanges(paste0("chr",
>>         c(1:22,
>>                          "X", "Y")), IRanges(1,
>>                          100))
>>                                          > seqlevels(gr) = sub("chr",
>> "ch",
>>                          seqlevels(gr))
>>
>>
>>                                To get a more comprehensive example I
>>         should have
>>                          followed my own
>>                          advice and
>>                                grabbed from the help page!
>>
>>                                      ## Rename:
>>                                      seqlevels(txdb) <- sub("chr", "",
>>                          seqlevels(txdb))
>>                                      seqlevels(txdb)
>>
>>                                      seqlevels(txdb) <- paste0("CH",
>>         seqlevels(txdb))
>>                                      seqlevels(txdb)
>>
>>
>>         seqlevels(txdb)[seqlevels(__****____txdb) ==
>>
>>
>>                          "CHM"] <- "M"
>>
>>
>>                                      seqlevels(txdb)
>>
>>
>>                                --
>>                                Computational Biology / Fred Hutchinson
>>         Cancer
>>                          Research Center
>>                                1100 Fairview Ave. N.
>>                                PO Box 19024 Seattle, WA 98109
>>
>>                                Location: Arnold Building M1 B861
>>                                Phone: (206) 667-2793
>>         <tel:%28206%29%20667-2793> <tel:%28206%29%20667-2793>
>>                          <tel:%28206%29%20667-2793>
>>
>>
>>
>>
>>                      --
>>                      Dr. Martin Morgan, PhD
>>
>>                      Fred Hutchinson Cancer Research Center
>>                      1100 Fairview Ave. N.
>>                      PO Box 19024 Seattle, WA 98109
>>
>>
>>                           [[alternative HTML version deleted]]
>>
>>                  ______________________________**_____________________
>>         Bioc-devel@r-project.org <mailto:Bioc-devel@r-project.**org<Bioc-devel@r-project.org>
>> >
>>         <mailto:Bioc-devel@r-project._**_org
>>         <mailto:Bioc-devel@r-project.**org <Bioc-devel@r-project.org>>>
>>                  mailing list
>>         https://stat.ethz.ch/mailman/_**___listinfo/bioc-devel<https://stat.ethz.ch/mailman/____listinfo/bioc-devel>
>>         <https://stat.ethz.ch/mailman/**__listinfo/bioc-devel<https://stat.ethz.ch/mailman/__listinfo/bioc-devel>
>> >
>>
>>
>>                  <https://stat.ethz.ch/mailman/**__listinfo/bioc-devel<https://stat.ethz.ch/mailman/__listinfo/bioc-devel>
>>         <https://stat.ethz.ch/mailman/**listinfo/bioc-devel<https://stat.ethz.ch/mailman/listinfo/bioc-devel>
>> >>
>>
>>
>>              --
>>              Hervé Pagès
>>
>>              Program in Computational Biology
>>              Division of Public Health Sciences
>>
>>              Fred Hutchinson Cancer Research Center
>>              1100 Fairview Ave. N, M1-B514
>>              P.O. Box 19024
>>              Seattle, WA 98109-1024
>>
>>              E-mail: hpages@fhcrc.org <mailto:hpages@fhcrc.org>
>>         <mailto:hpages@fhcrc.org <mailto:hpages@fhcrc.org>>
>>
>>              Phone: (206) 667-5791 <tel:%28206%29%20667-5791>
>>         <tel:%28206%29%20667-5791>
>>              Fax: (206) 667-1319 <tel:%28206%29%20667-1319>
>>         <tel:%28206%29%20667-1319>
>>
>>
>>
>>     --
>>     Hervé Pagès
>>
>>     Program in Computational Biology
>>     Division of Public Health Sciences
>>     Fred Hutchinson Cancer Research Center
>>     1100 Fairview Ave. N, M1-B514
>>     P.O. Box 19024
>>     Seattle, WA 98109-1024
>>
>>     E-mail: hpages@fhcrc.org <mailto:hpages@fhcrc.org>
>>     Phone: (206) 667-5791 <tel:%28206%29%20667-5791>
>>     Fax: (206) 667-1319 <tel:%28206%29%20667-1319>
>>
>>
>>
> --
> Hervé Pagès
>
> Program in Computational Biology
> Division of Public Health Sciences
> Fred Hutchinson Cancer Research Center
> 1100 Fairview Ave. N, M1-B514
> P.O. Box 19024
> Seattle, WA 98109-1024
>
> E-mail: hpages@fhcrc.org
> Phone:  (206) 667-5791
> Fax:    (206) 667-1319
>

	[[alternative HTML version deleted]]

