在 R 中映射基因型矩阵的有效代码答案

【问题标题】：Efficient code to map genotype matrix in R在 R 中映射基因型矩阵的有效代码
【发布时间】：2018-02-10 20:51:51
【问题描述】：

您好，我想将编码为三元组的基因型矩阵转换为编码为 0、1、2 的矩阵，即

c(0,0,1) <-> 0; c(0,1,0) <-> 1; c(0,0,1) <-> 2

首先这里是一些代码来生成需要减少的矩阵。

# generate genotypes
expand.G = function(n,p){
  probs = runif(n = p)
  G012.rows = matrix(rbinom(2,prob = probs,n=n*p),nrow = p)
  colnames(G012.rows) = paste('s',1:n,sep = '')
  rownames(G012.rows) = paste('g',1:p, sep = '')
  G012.cols = t(G012.rows)

  expand.geno = function(g){
    if(g == 0){return(c(1,0,0))}
    if(g == 1){return(c(0,1,0))}
    if(g == 2){return(c(0,0,1))}
  }

  gtype = c()
  for(i in 1:length(c(G012.cols))){
    gtype = c(
      gtype,
      expand.geno(c(G012.cols)[i])
    )
  }

  length(gtype)

  G = matrix(gtype,byrow = T, nrow = p)
  colnames(G) = paste('s',rep(1:n,each = 3),c('1','2','3'),sep = '')
  rownames(G) = paste('g',1:p, sep = '')
  print(G[1:10,1:15])
  print(G012.rows[1:10,1:5])

  return(G)
}

输出有 3n 列和 p 行，其中 n 是样本大小，p 是基因型数。现在我们可以使用以下函数将矩阵还原为 0,1,2 编码

reduce012 = function(x){
  if(identical(x, c(1,0,0))){
    return(0)
  } else if(identical(x, c(0,1,0))){
    return(1)
  } else if(identical(x,  c(0,0,1))){
    return(2)
  } else { 
    return(NA)
  }
}

reduce.G = function(G.gen){
  G.vec = 
    mapply(function(i,j) reduce012(as.numeric(G.gen[i,(3*j-2):(3*j)])), 
           i=expand.grid(1:(ncol(G.gen)/3),1:nrow(G.gen))[,2], 
           j=expand.grid(1:(ncol(G.gen)/3),1:nrow(G.gen))[,1]
    )

  G = matrix(G.vec, nrow = ncol(G.gen)/3, ncol = nrow(G.gen))
  colnames(G) = rownames(G.gen)
  return(G)
}

reduce.G.loop = function(G.gen){
  G = matrix(NA,nrow = ncol(G.gen)/3, ncol = nrow(G.gen))
  for(i in 1:nrow(G.gen)){
    for(j in 1:(ncol(G.gen)/3)){
      G[j,i] = reduce012(as.numeric(G.gen[i,(3*j-2):(3*j)]))
    }
  }
  colnames(G) = rownames(G.gen)
  return(G)
}

输出是 n 行 x p 列。顺便说一下，编码为 0,1,2 的矩阵是编码为三元组的矩阵的转置。

代码不是特别快。困扰我的是时间与n ^ 2有关。你能解释或提供更有效的代码吗？

G = expand.G(1000,20)
system.time(reduce.G(G))
system.time(reduce.G.loop(G))

G = expand.G(2000,20)
system.time(reduce.G(G))
system.time(reduce.G.loop(G))

G = expand.G(4000,20)
system.time(reduce.G(G))
system.time(reduce.G.loop(G))

【问题讨论】：

我无法运行您的代码，因为 all.equal(x, c(1,0,0)) 并不总是返回逻辑值。您可以使用isTRUE() 包装它，但请提供一个可重现的示例。
是的，这是以前的版本，我已更新为相同（）。我将编辑以反映以上内容。
如果你喜欢冒险，这里有一些 awk 代码我为此写的：github.com/vforget/gen2gemma/blob/master/gen2mgf.awk

标签： r performance matrix genome

【解决方案1】：

您可以简单地制作一个访问器查找表：

decode <- array(dim = c(3, 3, 3))
decode[cbind(1, 0, 0) + 1] <- 0
decode[cbind(0, 1, 0) + 1] <- 1
decode[cbind(0, 0, 1) + 1] <- 2

然后，就这样做：

matrix(decode[matrix(t(G + 1), ncol = 3, byrow = TRUE)], ncol = nrow(G))

这个完整的矢量化 R 版本将为您提供相同的矩阵，没有暗名和超快。

然而，如果你有更大的矩阵，你真的应该使用 Rcpp 来解决内存和时序问题。

【讨论】：

【解决方案2】：

这似乎比你的版本快大约三倍（重命名为reduce.G.orig）：

reduce.G <- function(G) {
  varmap = c("100"=0, "010"=1, "001"=2)
  result <- do.call(rbind, lapply(1:(ncol(G)/3)-1, function(val) 
    varmap[paste(G[,3*val+1], G[,3*val+2], G[,3*val+3], sep="")]))
  colnames(result) <- rownames(G)
  result
}

system.time(reduce.G(G))
#   user  system elapsed 
#  0.156   0.000   0.155 

system.time(reduce.G.orig(G))
#   user  system elapsed 
#  0.444   0.000   0.441 

identical(reduce.G(G), reduce.G.orig(G))
# [1] TRUE

【讨论】：