[BioC] Quickest way to convert IDs in a data frame?
James W. MacDonald
jmacdon at uw.edu
Thu Jul 25 23:57:35 CEST 2013
Hi Enrico and Herve,
This has to do with duplicate entries, but only when the duplicate entry
maps to many ENTREZID:
> select(org.Hs.eg.db, rep("ADORA2A", 4), "ENTREZID", "ALIAS")
ALIAS ENTREZID
1 ADORA2A 135
2 ADORA2A 135
3 ADORA2A 135
4 ADORA2A 135
> select(org.Hs.eg.db, rep("AGT", 4), "ENTREZID", "ALIAS")
ALIAS ENTREZID
1 AGT 183
2 AGT 189
Warning message:
In .generateExtraRows(tab, keys, jointype) :
'select' and duplicate query keys resulted in 1:many mapping between
keys and return rows
> select(org.Hs.eg.db, "AGT", "ENTREZID", "ALIAS")
ALIAS ENTREZID
1 AGT 183
2 AGT 189
Warning message:
In .generateExtraRows(tab, keys, jointype) :
'select' resulted in 1:many mapping between keys and return rows
So in the instances where a gene symbol maps to more than one ENTREZID,
the output gets truncated, whereas if it is a one-to-one mapping, it
does not.
Best,
Jim
On 7/25/2013 5:06 PM, Enrico Ferrero wrote:
> Hi,
>
> Hervé, that's exactly what I'm trying to say.
>
> Attached to this email is a tab delimited file with two columns of
> GeneSymbols (or Aliases), and here is some simple code to reproduce
> the unexpected behaviour:
>
> library(org.Hs.eg.db)
> mydf<- read.table("testdata.txt", sep="\t", header=TRUE, as.is=TRUE)
> mytest<- select(org.Hs.eg.db, key=mydf$GeneSymbol1, keytype="ALIAS",
> cols=c("SYMBOL","ENTREZID","ENSEMBL"))
> # check that mytest has less rows than mydf
> nrow(mydf)
> nrow(mytest)
> # pick a random row: they don't match
> mydf[250,]
> mytest[250,]
>
> Ideally, mytest should have the same number and position of rows of
> mydf so that I can then cbind them.
> If mytest has more rows because of multiple mappings that's also fine:
> I can always use merge(mydf, mytest), right?
>
> Thanks a lot to both for your help, it's very appreciated.
> Best,
>
>
> On 25 July 2013 21:32, Hervé Pagès<hpages at fhcrc.org> wrote:
>> Hi Enrico,
>>
>>
>> On 07/25/2013 01:20 PM, James W. MacDonald wrote:
>>> Hi Enrico,
>>>
>>> Please don't take things off-list (e.g., use reply-all).
>>>
>>>
>>> On 7/25/2013 2:17 PM, Enrico Ferrero wrote:
>>>> Hi James,
>>>>
>>>> Thanks very much for your help.
>>>> There is an issue that needs to be solved before thinking about what's
>>>> the best approach in my opinion.
>>>>
>>>> I don't understand why, but the object created with the call to select
>>>> (test in my example, first.two in yours) has a different number of
>>>> rows from the original object (df in my example). Specifically it has
>>>> *less* rows.
>>
>> I'm surprised it has less rows. It can definitely have more, when some
>> of the keys passed to select() are mapped to more than 1 row, but my
>> understanding was that select() would propagate unmapped keys to the
>> output by placing them in rows stuffed with NAs. So maybe I
>> misunderstood how select() works, or its behavior was changed, or
>> there is a bug somewhere. Could you please send the code that allows
>> us to reproduce this? Thanks.
>>
>> H.
>>
>>
>>> If all symbols were converted to all possible Entrez IDs,
>>>> I would expect it to have more rows, not less. To me, it looks like
>>>> not all rows are looked up and returned.
>>>>
>>>> Do you see what I mean?
>>>
>>> Sure. You could be using outdated gene symbols. Or perhaps you are using
>>> a mixture of symbols and aliases. Which is even cooler than just all
>>> symbols:
>>>
>>> > symb<- c(Rkeys(org.Hs.egSYMBOL)[1:10],
>>> Rkeys(org.Hs.egALIAS2EG)[31:45])
>>> > symb
>>> [1] "A1BG" "A2M" "A2MP1" "NAT1" "NAT2" "AACP"
>>> [7] "SERPINA3" "AADAC" "AAMP" "AANAT" "AAMP" "AANAT"
>>> [13] "DSPS" "SNAT" "AARS" "CMT2N" "AAV" "AAVS1"
>>> [19] "ABAT" "GABA-AT" "GABAT" "NPD009" "ABC-1" "ABC1"
>>> [25] "ABCA1"
>>> > select(org.Hs.eg.db, symb, "ENTREZID","SYMBOL")
>>> SYMBOL ENTREZID
>>> 1 A1BG 1
>>> 2 A2M 2
>>> 3 A2MP1 3
>>> 4 NAT1 9
>>> 5 NAT2 10
>>> 6 AACP 11
>>> 7 SERPINA3 12
>>> 8 AADAC 13
>>> 9 AAMP 14
>>> 10 AANAT 15
>>> 11 AAMP 14
>>> 12 AANAT 15
>>> 13 DSPS<NA>
>>> 14 SNAT<NA>
>>> 15 AARS 16
>>> 16 CMT2N<NA>
>>> 17 AAV<NA>
>>> 18 AAVS1 17
>>> 19 ABAT 18
>>> 20 GABA-AT<NA>
>>> 21 GABAT<NA>
>>> 22 NPD009<NA>
>>> 23 ABC-1<NA>
>>> 24 ABC1<NA>
>>> 25 ABCA1 19
>>> > select(org.Hs.eg.db, symb, "ENTREZID","ALIAS")
>>> ALIAS ENTREZID
>>> 1 A1BG 1
>>> 2 A2M 2
>>> 3 A2MP1 3
>>> 4 NAT1 9
>>> 5 NAT1 1982
>>> 6 NAT1 6530
>>> 7 NAT1 10991
>>> 8 NAT2 10
>>> 9 NAT2 81539
>>> 10 AACP 11
>>> 11 SERPINA3 12
>>> 12 AADAC 13
>>> 13 AAMP 14
>>> 14 AANAT 15
>>> 15 DSPS 15
>>> 16 SNAT 15
>>> 17 AARS 16
>>> 18 CMT2N 16
>>> 19 AAV 17
>>> 20 AAVS1 17
>>> 21 ABAT 18
>>> 22 GABA-AT 18
>>> 23 GABAT 18
>>> 24 NPD009 18
>>> 25 ABC-1 19
>>> 26 ABC1 19
>>> 27 ABC1 63897
>>> 28 ABCA1 19
>>> Warning message:
>>> In .generateExtraRows(tab, keys, jointype) :
>>> 'select' and duplicate query keys resulted in 1:many mapping between
>>> keys and return rows
>>> > mget(c("1982","6530","10991"), org.Hs.egGENENAME)
>>> $`1982`
>>> [1] "eukaryotic translation initiation factor 4 gamma, 2"
>>>
>>> $`6530`
>>> [1] "solute carrier family 6 (neurotransmitter transporter,
>>> noradrenalin), member 2"
>>>
>>> $`10991`
>>> [1] "solute carrier family 38, member 3"
>>>
>>> Best,
>>>
>>> Jim
>>>
>>>> On 25 July 2013 18:17, James W. MacDonald<jmacdon at uw.edu> wrote:
>>>>> Hi Enrico,
>>>>>
>>>>>
>>>>> On 7/25/2013 12:56 PM, Enrico Ferrero wrote:
>>>>>> Dear James,
>>>>>>
>>>>>> Thanks very much for your prompt reply.
>>>>>> I knew the problem was the for loop and the select function is indeed
>>>>>> a lot faster than that and works perfectly with toy data.
>>>>>>
>>>>>> However, this is what happens when I try to use it with real data:
>>>>>>
>>>>>>> test<- select(org.Hs.eg.db, keys=df$GeneSymbol, keytype="ALIAS",
>>>>>>> cols=c("SYMBOL","ENTREZID","ENSEMBL"))
>>>>>> Warning message:
>>>>>> In .generateExtraRows(tab, keys, jointype) :
>>>>>> 'select' and duplicate query keys resulted in 1:many mapping
>>>>>> between
>>>>>> keys and return rows
>>>>>>
>>>>>> which is probably the warning you mentioned.
>>>>>
>>>>> That's not the warning I mentioned, but it does point out the same
>>>>> issue,
>>>>> which is that there is a one to many mapping between symbol and
>>>>> entrez gene
>>>>> ID.
>>>>>
>>>>> So now you have to decide if you want to be naive (or stupid,
>>>>> depending on
>>>>> your perspective) or not. You could just cover your eyes and do this:
>>>>>
>>>>> first.two<- first.two[!duplicated(first.two$SYMBOL),]
>>>>>
>>>>> which will choose for you the first symbol -> gene ID mapping and
>>>>> nuke the
>>>>> rest. That's nice and quick, but you are making huge assumptions.
>>>>>
>>>>> Or you could decide to be a bit more sophisticated and do something like
>>>>>
>>>>> thelst<- tapply(1:nrow(first.two), first.two$SYMBOL, function(x)
>>>>> first.two[x,])
>>>>>
>>>>> At this point you can take a look at e.g., thelst[1:10] to see what
>>>>> we just
>>>>> did
>>>>>
>>>>> thelst<- do.call("rbind", lapply(thelst, function(x) c(x[1,1],
>>>>> paste(x[,2],
>>>>> collapse = "|")))
>>>>>
>>>>> and here you can look at head(thelst).
>>>>>
>>>>> Then you can check to ensure that the first column of thelst is
>>>>> identical to
>>>>> the first column of df, and proceed as before.
>>>>>
>>>>> But there is still the problem of the multiple mappings. As an example:
>>>>>
>>>>>> thelst[1:5]
>>>>> $HBD
>>>>> SYMBOL ENTREZID
>>>>> 2535 HBD 3045
>>>>> 2536 HBD 100187828
>>>>>
>>>>> $KIR3DL3
>>>>> SYMBOL ENTREZID
>>>>> 17513 KIR3DL3 115653
>>>>> 17514 KIR3DL3 100133046
>>>>>
>>>>>> mget(as.character(thelst[[1]][,2]), org.Hs.egGENENAME)
>>>>> $`3045`
>>>>> [1] "hemoglobin, delta"
>>>>>
>>>>> $`100187828`
>>>>> [1] "hypophosphatemic bone disease"
>>>>>
>>>>>> mget(as.character(thelst[[2]][,2]), org.Hs.egGENENAME)
>>>>> $`115653`
>>>>> [1] "killer cell immunoglobulin-like receptor, three domains, long
>>>>> cytoplasmic tail, 3"
>>>>>
>>>>> $`100133046`
>>>>> [1] "killer cell immunoglobulin-like receptor three domains long
>>>>> cytoplasmic
>>>>> tail 3"
>>>>>
>>>>>
>>>>> So HBD is the gene symbol for two different genes! If this gene
>>>>> symbol is in
>>>>> your data, you will now have attributed your data to two genes that
>>>>> apparently are not remotely similar. if KIR3DL3 is in your data, then it
>>>>> worked out OK for that gene.
>>>>>
>>>>> Best,
>>>>>
>>>>> Jim
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>> The real problem is that the number of rows is now different for the 2
>>>>>> objects:
>>>>>>> nrow(df); nrow(test)
>>>>>> [1] 573
>>>>>> [1] 201
>>>>>>
>>>>>> So I obviously can't put the new data into the original df. My
>>>>>> impression is that when the 1 to many mapping arises, the select
>>>>>> functions exits, with that warning message. As a result, my test
>>>>>> object is incomplete.
>>>>>>
>>>>>> On top of that, and I can't really explain this, the row positions are
>>>>>> messed up, e.g.
>>>>>>
>>>>>>> all.equal(df[100,],test[100,])
>>>>>> returns FALSE.
>>>>>>
>>>>>> How can I work around this?
>>>>>>
>>>>>> Thanks a lot!
>>>>>>
>>>>>> Best,
>>>>>>
>>>>>> On 25 July 2013 16:58, James W. MacDonald<jmacdon at uw.edu> wrote:
>>>>>>> Hi Enrico,
>>>>>>>
>>>>>>>
>>>>>>> On 7/25/2013 11:35 AM, Enrico Ferrero wrote:
>>>>>>>> Hello,
>>>>>>>>
>>>>>>>> I often have data frames where I need to perform ID conversions on
>>>>>>>> one
>>>>>>>> or
>>>>>>>> more of the columns while preserving the order of the rows, e.g.:
>>>>>>>>
>>>>>>>> GeneSymbol Value1 Value2
>>>>>>>> GS1 2.5 0.1
>>>>>>>> GS2 3 0.2
>>>>>>>> ..
>>>>>>>>
>>>>>>>> And I want to obtain:
>>>>>>>>
>>>>>>>> GeneSymbol EntrezGeneID Value1 Value2
>>>>>>>> GS1 EG1 2.5 0.1
>>>>>>>> GS2 EG2 3 0.2
>>>>>>>> ..
>>>>>>>>
>>>>>>>> What I've done so far was to create a function that uses
>>>>>>>> org.Hs.eg.db to
>>>>>>>> loop over the rows of the column and does the conversion:
>>>>>>>>
>>>>>>>> library(org.Hs.eg.db)
>>>>>>>> alias2EG<- function(x) {
>>>>>>>> for (i in 1:length(x)) {
>>>>>>>> if (!is.na(x[i])) {
>>>>>>>> repl<- org.Hs.egALIAS2EG[[x[i]]][1]
>>>>>>>> if (!is.null(repl)) {
>>>>>>>> x[i]<- repl
>>>>>>>> }
>>>>>>>> else {
>>>>>>>> x[i]<- NA
>>>>>>>> }
>>>>>>>> }
>>>>>>>> }
>>>>>>>> return(x)
>>>>>>>> }
>>>>>>>
>>>>>>> I should first note that gene symbols are not unique, so you are
>>>>>>> taking a
>>>>>>> chance on your mappings. Is there no other annotation for your data?
>>>>>>>
>>>>>>> In addition, you should note that it is almost always better to
>>>>>>> think of
>>>>>>> objects as vectors and matrices in R, rather than as things that
>>>>>>> need to
>>>>>>> be
>>>>>>> looped over (e.g., R isn't Perl or C).
>>>>>>>
>>>>>>> first.two<- select(org.Hs.eg.db, as.character(df$GeneSymbol),
>>>>>>> "ENTREZID",
>>>>>>> "SYMBOL")
>>>>>>>
>>>>>>> Note that there used to be a warning or an error (don't remember
>>>>>>> which)
>>>>>>> when
>>>>>>> you did something like this, stating that gene symbols are not unique,
>>>>>>> and
>>>>>>> that you shouldn't do this sort of thing. Apparently this warning has
>>>>>>> been
>>>>>>> removed, but the issue remains valid.
>>>>>>>
>>>>>>> ## check yourself
>>>>>>>
>>>>>>> all.equal(df$GeneSymbol, first.two$SYMBOL)
>>>>>>>
>>>>>>> ## if true, proceed
>>>>>>>
>>>>>>> df<- data.frame(first.two, df[,-1])
>>>>>>>
>>>>>>> Best,
>>>>>>>
>>>>>>> Jim
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>> and then call the function like this:
>>>>>>>>
>>>>>>>> df$EntrezGeneID<- alias2GS(df$GeneSymbol)
>>>>>>>>
>>>>>>>> This works well, but gets very slow when I need to do multiple
>>>>>>>> conversions
>>>>>>>> on large datasets.
>>>>>>>>
>>>>>>>> Is there any way I can achieve the same result but in a quicker, more
>>>>>>>> efficient way?
>>>>>>>>
>>>>>>>> Thank you.
>>>>>>>>
>>>>>>> --
>>>>>>> James W. MacDonald, M.S.
>>>>>>> Biostatistician
>>>>>>> University of Washington
>>>>>>> Environmental and Occupational Health Sciences
>>>>>>> 4225 Roosevelt Way NE, # 100
>>>>>>> Seattle WA 98105-6099
>>>>>>>
>>>>> --
>>>>> James W. MacDonald, M.S.
>>>>> Biostatistician
>>>>> University of Washington
>>>>> Environmental and Occupational Health Sciences
>>>>> 4225 Roosevelt Way NE, # 100
>>>>> Seattle WA 98105-6099
>>>>>
>>>>
>> --
>> 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
>
>
--
James W. MacDonald, M.S.
Biostatistician
University of Washington
Environmental and Occupational Health Sciences
4225 Roosevelt Way NE, # 100
Seattle WA 98105-6099
More information about the Bioconductor
mailing list