r-转换/更新表中基因ID的最快方法

注意：我不是在问Bioconductor特定的问题，但我需要在示例代码中使用Bioconducter。请耐心听我说

你好，

我有许多以制表符分隔的文件，其中包含关于特定基因的各种类型的信息。一个或多个列可以是基因符号的别名，我需要将其升级到最新的基因符号注释。

我使用的是Bioconductor的org.Hs.eg.db库（尤其是org.Hs.egALIAS2EG和org.Hs.eg SYMBOL对象）。

报告的代码完成了这项工作，但速度非常慢，我想这是因为嵌套的for循环在每次迭代时都会查询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.eg SYMBOL是对象，而不是函数。

谢谢！

使用mget函数。

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

等等。这就是它应该使用的方式，而且比其他替代方案快得多。然而，也许有必要首先将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.vvector（data.frame）不起作用的原因：data.frame不是矩阵。这是一个以奇特方式显示的列表，上面定义了[]函数。

您可以使用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)

只是一个快速警告：一个别名可以映射到多个Entrez Gene 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是"正确的"。因此，你最安全的选择是在表中获得每个别名的所有可能的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 Gene 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

这是我能想到的最好的。

首先写一个函数：

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)

相关内容

最新更新

热门标签：