[Bioc-devel] deprecation of keepSeqlevels and renameSeqlevels

Sounds great. Could we have a more intuitive name? Like naturalOrder(). And
then the method to fix the GRanges could be restoreNaturalOrder(). Just
some suggestions...




On Tue, Sep 17, 2013 at 10:58 AM, Herv? Pag?s <hpages at fhcrc.org> wrote:

For sorting the chromosome names in "natural" order, you can use
the makeSeqnameIds() helper function in GenomicRanges:

 seqlevels(gr) <- seqlevels(gr)[makeSeqnameIds(**seqlevels(gr))]

It recognizes several naming conventions (chr-prefixed or not,
roman, etc...) e.g.:

makeSeqnameIds(c("Y", "1", "9", "M", "2"))

 [1] 4 1 3 5 2

makeSeqnameIds(c("chrY", "chrI", "chrIX", "chrM", "chrII"))

 [1] 4 1 3 5 2

I still need to improve the documentation of this function, put
more examples, and add unit tests. It's used internally by the
makeTranscriptDb* functions in GenomicFeatures to assign internal
ids to the chromosomes so that all the TranscriptDb extractors
can then return GRanges and GRangesList objects with the seqlevels
in the "natural" order.

Cheers,
H.



On 09/17/2013 10:23 AM, Kasper Daniel Hansen wrote:

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

        <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 at fhcrc.org <mailto:hpages at 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 at fhcrc.org
Phone:  (206) 667-5791
Fax:    (206) 667-1319

        [[alternative HTML version deleted]]

        [[alternative HTML version deleted]]

       [[alternative HTML version deleted]]

--
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 at fhcrc.org
Phone:  (206) 667-5791
Fax:    (206) 667-1319

   [[alternative HTML version deleted]]