Hi Marc, Herve,
I also noticed this behaviour:
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
seqlevels(txdb, force=TRUE) <- c("chr1", "ch2")
oldLevs <- seqlevels(txdb)
seqlevels(txdb, force=TRUE) <- "chr1"
seqlevels(txdb, force=TRUE) <- oldLevs
Error in .seqinfo.TranscriptDbReplace(x, new2old = new2old, force = force,
:
The replacement value must be either a 1 to 1 replacement or a subset of
the original set when replacing the 'seqinfo' of a TranscriptDb object
But:
restoreSeqlevels(txdb)
seqlevels(txdb, force=TRUE) <- oldLevs
This is probably intentional, but I found it to be rather confusing. If
one should be able to use the seqlevels replacement method to control
active and inactive chromosomes in the TranscriptDb object, then having
this mandatory step of restoring all chromosomes to the active state
followed by another restriction to be quite cumbersome. At least a
somewhat more useful error message would help in that case. Or couldn't
the replacement method always call restoreSeqlevels internally.
It would also help to point out somewhere (I may have missed it) that
TranscriptDb seqlevels (or rather the internal GRanges) are actually pass
by reference:
txdb2 <- txdb
seqlevels(txdb)
seqlevels(txdb2)
restoreSeqlevels(txdb)
seqlevels(txdb2)
Also Herve: what is the definition of a "used" seqlevel? Does that mean
that a factor level is defined, but no range in the GRanges object is
assigned this level? Why would those exists in a TranscriptDb object at
all?
Florian
On 9/18/13 9:08 PM, "Hervé Pagès" <hpa...@fhcrc.org> wrote:
>Note that currently it's kind of inconsistent anyway because it doesn't
>look at the seqlevels that are in use. For example:
>
> library(TxDb.Hsapiens.UCSC.hg19.knownGene)
> txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
>
>Trying to drop chrUn_gl000249 (which I know is not in use):
>
> > seqlevels(txdb) <- setdiff(seqlevels(txdb), "chrUn_gl000249")
> Error in .seqinfo.TranscriptDbReplace(x, new2old = new2old, force =
>force, :
> You need to use force=TRUE if you want to drop seqlevels.
>
>'force=TRUE' should only be required when trying to drop seqlevels
>that are in use.
>
>So the choice are more between: (a) consistent, (b) convenient,
>or (c) leave it as it is (which is neither convenient or consistent).
>
>I vote for (b).
>
>H.
>
>
>On 09/18/2013 11:31 AM, Marc Carlson wrote:
>> I actually considered this, but I opted to do it this way just for the
>> sake of being consistent (which was my whole mission for implementing
>> seqlevels in here in the 1st place). Now I could make it more
>> convenient here and break consistency with how it is used elsewhere, but
>> what do people prefer?
>>
>> Consistent or convenient?
>>
>>
>> Marc
>>
>>
>>
>> On 09/18/2013 10:40 AM, Hervé Pagès wrote:
>>> Hi Marc,
>>>
>>> Wouldn't it make sense to just ignore the 'force' arg when
>>> dropping the seqlevels of a TranscriptDb?
>>>
>>> The 'force' argument is FALSE by default and this prevents
>>> seqlevels<- to shrink GRanges or other vector-like objects
>>> when the user tries to drop seqlevels that are in use.
>>> Internally seqlevels<- calls seqlevelsInUse() to get the
>>> seqlevels currently in use and see if they intersect with
>>> the seqlevels to drop.
>>>
>>> In the TranscriptDb situation, people always have to use
>>> 'force=TRUE' to drop seqlevels, regardless of whether the
>>> levels to drop are in use or not (the seqlevelsInUse()
>>> getter not being defined for TranscriptDb objects, I suspect
>>> seqlevels<- doesn't look at this).
>>>
>>> So maybe 'force' could just be ignored for TranscriptDb objects?
>>> That would make seqlevels<- a little bit more user-friendly on
>>> those objects.
>>>
>>> Thanks,
>>> H.
>>>
>>>
>>> On 09/13/2013 10:38 AM, Marc Carlson wrote:
>>>> Hi Florian,
>>>>
>>>> Yes we are trying to make things more uniform. seqlevels() lets you
>>>> rename as well as deactivate chromosomes you want to ignore, so it was
>>>> really redundant with isActiveSeq(). So we are moving away from
>>>> isActiveSeq() just so that users have less to learn about. The reason
>>>> why isActiveSeq was different from seqlevels was just because it was
>>>> born for a TranscriptDb (which is based on an annotation database)
>>>> instead of being born on a GRanges object. So seqlevels was the more
>>>> general tool.
>>>>
>>>> Marc
>>>>
>>>>
>>>>
>>>> On 09/13/2013 07:24 AM, Hahne, Florian wrote:
>>>>> Hi Marc,
>>>>> I saw these warnings in Gviz, but they stem from GenomicFeatures
>>>>>
>>>>> Warning messages:
>>>>> 1: 'isActiveSeq' is deprecated.
>>>>> Use 'seqlevels' instead.
>>>>> See help("Deprecated") and help("GenomicFeatures-deprecated").
>>>>> 2: 'isActiveSeq' is deprecated.
>>>>> Use 'seqlevels' instead.
>>>>> See help("Deprecated") and help("GenomicFeatures-deprecated").
>>>>> 3: 'isActiveSeq<-' is deprecated.
>>>>> Use 'seqlevels' instead.
>>>>> See help("Deprecated") and help("GenomicFeatures-deprecated").
>>>>> 4: 'isActiveSeq<-' is deprecated.
>>>>> Use 'seqlevels' instead.
>>>>> See help("Deprecated") and help("GenomicFeatures-deprecated").
>>>>> 5: 'isActiveSeq' is deprecated.
>>>>> Use 'seqlevels' instead.
>>>>> See help("Deprecated") and help("GenomicFeatures-deprecated").
>>>>> 6: 'isActiveSeq<-' is deprecated.
>>>>> Use 'seqlevels' instead.
>>>>> See help("Deprecated") and help("GenomicFeatures-deprecated").
>>>>>
>>>>> So has the whole idea of active chromosomes in the data base been
>>>>> dropped? I could not find anything in the change notes. Do I get it
>>>>> right that you can now do
>>>>> seqlevels(txdb, force=TRUE) <- "chr1"
>>>>> if you just want the first chromosome to be active?
>>>>>
>>>>> Florian
>>>>>
>>>>
>>>>
>>>> [[alternative HTML version deleted]]
>>>>
>>>> _______________________________________________
>>>> Bioc-devel@r-project.org mailing list
>>>> 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: hpa...@fhcrc.org
>Phone: (206) 667-5791
>Fax: (206) 667-1319
_______________________________________________
Bioc-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/bioc-devel
