[Bioc-devel] deprecation of keepSeqlevels and renameSeqlevels

Actually, what I (and I think several others) just want is
   fixSeqnames()
which sorts and fixes them to some convention.  I don't really care which
one, but I want to do some easy harmonization, preferably in line with
other bioc tools.  I know this is hard to do for all organisms, but having
something that deals with the major ones (human, mouse, fruit fly, yeast,
some monkeys) would be extremely valuable.  If this function existed, I
could just throw all my GRanges at it.

Kasper


On Tue, Sep 17, 2013 at 1:14 PM, Michael Lawrence <lawrence.michael at gene.com

wrote:

When I hit this I usually just do:

seqlevels(gr) <- paste0("chr", c(1:22, "X", "Y"))

It would be nice if there were a function for sorting chromosome names
according to a specified naming convention. Like sortSeqlevels(gr,
UCSCNaming()) or the alternative replacement syntax: seqinfo(gr) <-
sort(seqinfo(gr), UCSCNaming()). Although sort is only single-dispatch.

Michael






On Tue, Sep 17, 2013 at 8:59 AM, Vincent Carey
<stvjc at channing.harvard.edu>wrote:

I am replying to this as Michael mentions that this is a powerful
operation.  Here's my
problem that I believe is related, but I do not see a straightforward
solution

bhmm19uni

GRanges with 559456 ranges and 4 metadata columns:
          seqnames                 ranges strand   |              name
score
             <Rle>              <IRanges>  <Rle>   |       <character>
<numeric>
        1     chr1 [ 67229025,  67341224]      *   | 13_Heterochrom/lo
     0
        2     chr1 [203057955, 203060154]      *   |      12_Repressed
     0
        3     chr1 [  8458427,   8486226]      *   | 13_Heterochrom/lo
     0
        4     chr1 [ 16904427,  16904826]      *   | 5_Strong_Enhancer
     0
        5     chr1 [ 25289427,  25293426]      *   |       11_Weak_Txn
     0
      ...      ...                    ...    ... ...               ...
   ...
   570766    chr22   [51100931, 51101330]      *   |   2_Weak_Promoter
     0
   570767    chr22   [51101331, 51101530]      *   |   6_Weak_Enhancer
     0
   570768    chr22   [51101531, 51101730]      *   |   2_Weak_Promoter
     0
   570769    chr22   [51101731, 51178130]      *   | 13_Heterochrom/lo
     0
   570770    chr22   [51178131, 51178330]      *   | 15_Repetitive/CNV
     0
                           thick     numcode
                       <IRanges> <character>
        1 [ 67001613,  67113812]          13
        2 [201324578, 201326777]          12
        3 [  8381014,   8408813]          13
        4 [ 16777014,  16777413]           5
        5 [ 25162014,  25166013]          11
      ...                    ...         ...
   570766   [49447797, 49448196]           2
   570767   [49448197, 49448396]           6
   570768   [49448397, 49448596]           2
   570769   [49448597, 49524996]          13
   570770   [49524997, 49525196]          15
   ---
   seqlengths:
         chr1     chr10     chr11      chr5 ...      chr2     chr21
  chr4
    249250621 135534747 135006516 180915260 ... 243199373  48129895
191154276

If I try to make a karyogram with ggbio, the chromosomes are in an
unnatural order.  What is the right approach to getting the seqinfo in a
natural order with respect to chromosome names and lengths?  Reassigning
seqlevels seems dangerous but I haven't experimented fully.  It must be a
common use case but I do not see it in doc.



On Tue, May 21, 2013 at 11:00 AM, Michael Lawrence <
lawrence.michael at gene.com> wrote:

On Mon, May 20, 2013 at 10:36 PM, Herv? Pag?s <hpages at 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 at fhcrc.org
<mailto:hpages at 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 at fhcrc.org <mailto:mtmorgan at 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. 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. 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 at fhcrc.org <mailto:mtmorgan at fhcrc.org>
                 <mailto:mtmorgan at fhcrc.org <mailto:

mtmorgan at 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>




             --
             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 at r-project.org>

         mailing list

         [[alternative HTML version deleted]]

	[[alternative HTML version deleted]]