轉換/更新表中基因ID的最快方法？

Question

注：我不是問一個Bioconductor特定的問題，但我需要示例代碼中的Bioconductor。請跟我來。

嗨，

我有許多包含各種類型的關於特定基因信息製表符分隔的文件。一個或多個列可以是基因符號別名，我需要升級到最新的基因符號註釋。

我使用Bioconductor的org.Hs.eg.db庫來執行此操作（尤其是org.Hs.egALIAS2EG和org.Hs.egSYMBOL對象）。

代碼報告做的工作，但速度非常慢，我猜是因爲在每次迭代查詢org.Hs.eg.db數據庫的嵌套循環。是否有更快/更簡單/更智能的方法來達到相同的效果？

library(org.Hs.eg.db) 

myTable <- read.table("tab_delimited_file.txt", header=TRUE, sep="\t", as.is=TRUE) 

for (i in 1:nrow(myTable)) { 
    for (j in 1:ncol(myTable)) { 
     repl <- org.Hs.egALIAS2EG[[myTable[i,j]]][1] 
     if (!is.null(repl)) { 
      repl <- org.Hs.egSYMBOL[[repl]][1] 
      if (!is.null(repl)) { 
       myTable[i,j] <- repl 
      } 
     } 
    } 
} 

write.table(myTable, file="new_tab_delimited_file", quote=FALSE, sep="\t", row.names=FALSE, col.names=TRUE) 

我想使用的應用功能之一，但請記住，org.Hs.egALIAS2EG和org.Hs.egSYMBOL是對象，而不是功能。

謝謝！

Answer 1

這是我能想到的最好的。

首先寫一個函數：

alias2GS <- function(x) { 
    for (i in 1:length(x)) { 
     if (!is.na(x[i])) { 
      repl <- org.Hs.egALIAS2EG[[x[i]]][1] 
      if (!is.null(repl)) { 
       repl <- org.Hs.egSYMBOL[[repl]][1] 
       if (!is.null(repl)) { 
        x[i] <- repl 
       } 
      } 
     } 
    } 
    return(x) 
} 

，然後調用函數對於需要轉換的數據幀的每一列，就像這樣：

df$GeneSymbols <- alias2GS(df$GeneSymbols) 

Answer 2

可以sapply使用和名稱的幾個變量不屬於載體，如你的對象從org.Hs.eg.db庫：

library(org.Hs.eg.db) 
myTable <- read.table("tab_delimited_file.txt", header=TRUE, sep="\t", as.is=TRUE) 

myfunc <- function(idx,mytab,a2e,es){ 
      i = idx %/% nrow(mytab) + 1 
      j = idx %% ncol(mytab) + 1 
      repl <- a2e[[myTable[i,j]]][1]; 
      if (!is.null(repl)) { 
       repl <- es[[repl]][1] 
       if (!is.null(repl)) { 
       return(repl) 
       } 
      } 
      else {return("NA")} 
      } 

vec <- 0:(ncol(myTable)*nrow(myTable)-1) 
out <- sapply(vec,mytab=myTable,a2e=org.Hs.egALIAS2EG,es=org.Hs.egSYMBOL,myfunc) 
myTable <- matrix(out, nrow=nrow(myTable),ncol=ncol(myTable),byrow=T) 

Answer 3

使用mget功能。

eg[i,] <- mget(myTable[i,], org.Hs.egALIAS2EG) 
symbol[i, ] <- mget(myTable[i,], org.Hs.egSYMBOL) 

等，這是它應該被使用的方式和比ENY其他替代快得多。但是，也許這將是歡顏重塑爲myTable到基因名稱第一的向量：

v <- unique(as.vector(as.matrix(myTable))) 
v <- v[ v %in% ls(org.Hs.egALIAS2EG) ] 
eg <- unlist(mget(v, org.Hs.egALIAS2EG)) 
symbol <- unlist(mget(eg, org.Hs.egSYMBOL)) 

等上面的第二行確保您只查找實際上是在數據庫中的符號。現在您可以使用符號表再次修改表格。假設不是myTable匹配的所有元素，這是一種可以做到的方法。爲簡潔起見，我將表複製到t：

t <- as.matrix(myTable) 
names(symbol) <- v 
t[ !is.na(match(t, v)) ] <- symbol[ match(t, v) ][ ! is.na(match(t, v)) ] 

好的。假設我們正在使用矩陣（或多或少）字符。然而，坦率地說，你只有一個有兩列的數據框架，所以不需要真正自動化代碼就好像你有數百列。讓我們寫一個小函數。 （這將是簡單的，如果我們可以假設，表中的所有元素可以在org.Hs.egALIAS2EG找到）

convert2symbol <- function(x) { 
    v <- unique(as.character(x)) 
    v <- v[ v %in% ls(org.Hs.egALIAS2EG) ] 
    eg <- unlist(mget(v, org.Hs.egALIAS2EG)) 
    symbol <- unlist(mget(eg, org.Hs.egSYMBOL)) 
    m <- match(x, v) 
    v[ ! is.na(m) ] <- symbol[ m ][ ! is.na(m) ] 
    v 
} 

現在你可以

myTable$LigandGene <- convert2symbol(myTable$LigandGene) 

或

newTable <- apply(myTable, 2, convert2symbol) 

至於爲什麼as.vector（data.frame）不起作用：data.frame不是矩陣。這是一個列表，以一種奇特的方式顯示，並定義了其中的[]函數。

Answer 4

只是一個快速的警告：一個別名可以映射到多個Entrez基因ID。

因此，您當前的解決方案假定第一個列出的ID是正確的（這可能不是正確的）。

# e.g. The alias "A1B" is assumed to map to "1" and not "6641" 
mget("A1B", org.Hs.egALIAS2EG) 
# $A1B 
# [1] "1" "6641" 

如果您?org.Hs.egALIAS2EG退房的幫助下，你會看到，這是從來沒有建議使用別名或符號作爲主要基因標識。

## From the 'Details' section of the help: 
# Since gene symbols are sometimes redundantly assigned in the literature, 
# users are cautioned that this map may produce multiple matching results 
# for a single gene symbol. Users should map back from the entrez gene IDs 
# produced to determine which result is the one they want when this happens. 

# Because of this problem with redundant assigment of gene symbols, 
# is it never advisable to use gene symbols as primary identifiers. 

沒有人工調整，不可能知道哪個ID是「正確的」。因此，你的最簡單的辦法是讓所有可能的標識和符號爲表中的每個別名，同時保持相關信息哪些是受體，哪些是配體：

# your example subset with "A1B" and "trash" added for complexity 
myTable <- data.frame(
    ReceptorGene = c("A1B", "ACVR2B", "ACVR2B", "ACVR2B", "ACVR2B", "AMHR2", "BLR1", "BMPR1A", "BMPR1A", "BMPR1A", "BMPR1A", "BMPR1A"), 
    LigandGene = c("trash", "INHA", "INHBA", "INHBB", "INHBC", "AMH", "SCYB13", "BMP10", "BMP15", "BMP2", "BMP3", "BMP4"), 
    stringsAsFactors = FALSE 
) 

# unlist and rename 
my.aliases <- unlist(myTable) 
names(my.aliases) <- paste(names(my.aliases), my.aliases, sep = ".") 

# determine which aliases have a corresponding Entrez Gene ID 
has.key <- my.aliases %in% keys(org.Hs.egALIAS2EG) 

# replace Aliases with character vectors of all possible entrez gene IDs 
my.aliases[has.key] <- sapply(my.aliases[has.key], function(x) { 
    eg.ids <- unlist(mget(x, org.Hs.egALIAS2EG)) 
    symbols <- unlist(mget(eg.ids, org.Hs.egSYMBOL)) 
}) 

# my.aliases retains all pertinent information regarding the original alias 
my.aliases[1:3] 
# $ReceptorGene1.A1B 
#  1 6641 
# "A1BG" "SNTB1" 
# 
# $ReceptorGene2.ACVR2B 
#  93 
# "ACVR2B" 
# 
# $ReceptorGene3.ACVR2B 
#  93 
# "ACVR2B" 

一旦你知道哪些Entrez基因ID是適當的，您可以將它們作爲附加列存儲在桌面上。

myTable$receptor.id <- c("1", "93", "93", "93", "93", "269", "643", "657", "657", "657", "657", "657") 
myTable$ligand.id <- c(NA, "3623", "3624", "3625", "3626", "268", "10563", "27302", "9210", "650", "651", "652") 

然後，當您需要更新到最新的符號時，您可以使用Entrez基因ID（永遠不需要更新）。

has.key <- myTable$receptor.id %in% keys(org.Hs.egSYMBOL) 
myTable$ReceptorGene[has.key] <- unlist(mget(myTable$receptor.id[has.key], org.Hs.egSYMBOL)) 

has.key <- myTable$ligand.id %in% keys(org.Hs.egSYMBOL) 
myTable$LigandGene[has.key] <- unlist(mget(myTable$ligand.id[has.key], org.Hs.egSYMBOL)) 

head(myTable) 
# ReceptorGene LigandGene receptor.id ligand.id 
# 1   A1BG  trash   1  <NA> 
# 2  ACVR2B  INHA   93  3623 
# 3  ACVR2B  INHBA   93  3624 
# 4  ACVR2B  INHBB   93  3625 
# 5  ACVR2B  INHBC   93  3626 
# 6  AMHR2  AMH   269  268