求助：字符型矩阵的运算

有一个基因数据矩阵mat，约20000行*300列，格式如下，取值为A、T、G、C四个值，其中T为高频值，其余均为低频值。
现在需要做类似矩阵相乘的运算：mat %*% t(mat)，得到一个20000*20000的矩阵matrix，运算时定义若两个值相等，则乘积为0，如A*A=0；若两个值不等，则乘积为1，如A*T=1。
例如：

请教高手在R语言里面有没有办法实现，万分感谢。


数据示例：

MrkID	y1	y2	y3	y4
x1	C	T	T	T
x2	T	T	T	T
x3	T	T	T	T
x4	T	T	T	T
x5	T	C	T	T
x6	T	T	T	T
x7	T	T	T	T
x8	T	T	T	T
x9	T	T	T	T
x10	T	T	T	T
x11	G	T	T	T
x12	T	T	T	T
x13	G	T	T	T
x14	T	T	T	T
x15	T	T	T	A
x16	T	T	T	T
x17	T	T	T	T
x18	T	G	T	T
x19	T	T	T	T

1. dat <- matrix(sample(c("A", "T", "G", "C"), 10*3, replace = TRUE), ncol = 3)
   
2. 
   
3. pairwise.comparison <- matrix(nrow = nrow(dat), ncol = nrow(dat))
   
4. for (j in 1:nrow(dat)) {
   
5.   equal <- (dat[j, ] == t(dat))
   
6.   equal.cnt <- colSums(equal)
   
7.   pairwise.comparison[j, ] <- equal.cnt
   
8. }

复制代码

qoiqpwqr 发表于 2018-4-13 22:30

厉害，谢谢高手。还有一个问题，由于我的数据有20000行，这个for循环耗时比较长，而我这个计算矩阵的过程又会被循环使用很多次，这样就会耗时很长，有没有办法提升这个for循环的效率，或者可否使用apply族函数来代替这个for循环？

kanzongxuan 发表于 2018-4-15 09:51
厉害，谢谢高手。还有一个问题，由于我的数据有20000行，这个for循环耗时比较长，而我这个计算矩阵的过程 ...

程序跑成功没有？

2楼的程序跑成功了，只是耗时比较长，完成矩阵的输出近2个小时，还在琢磨怎样提高效率，因为要反复循环算很多个类似的矩阵。

cheetahfly 发表于 2018-4-17 15:16
程序跑成功没有？

2楼的程序跑成功了，只是耗时比较长，完成矩阵的输出近2个小时，还在琢磨怎样提高效率，因为要反复循环算很多个类似的矩阵。

kanzongxuan 发表于 2018-4-17 15:41
2楼的程序跑成功了，只是耗时比较长，完成矩阵的输出近2个小时，还在琢磨怎样提高效率，因为要反复循环算 ...

如果你的电脑内存足够大，到底多大我没测算过，反正我的16GB内存没有跑成功20000*300，你可以试试下面的代码：

1. . <- unique(as.vector(dat))
   
2. . <- lapply(., function(x) tcrossprod(dat == x))
   
3. . <- Reduce(`+`, .)
   
4. result1 <- ncol(dat) - .

复制代码

我用2000*300的数据测算过，如果2楼的代码耗时为1的话，上面的方法大概耗时0.013。

如果你电脑的内存不大的话，就可能要忍受较长的执行时间了，可以试试用C++代码的函数：

1. library(Rcpp)
   
2. cppFunction('NumericMatrix charMatProd(CharacterMatrix x, CharacterMatrix y) {
   
3.     int nrow = x.nrow(), ncol = x.ncol();
   
4.     NumericMatrix out(nrow, nrow);
   
5.     for (int i = 0; i < nrow; i++) {
   
6.         for (int j = 0; j < nrow; j++) {
   
7.             int total = 0;
   
8.             for (int k = 0; k < ncol; k ++) {
   
9.                 if (x(i, k) != y(k, j)) {
   
10.                     total += 1;
    
11.                 }
    
12.             }
    
13.             out(i, j) = total;
    
14.         }
    
15.     }
    
16.     return out;
    
17. }')
    
18. result2 <- charMatProd(dat, t(dat))

复制代码

用2000*300的数据测算过，如果2楼的代码耗时为1的话，上面的方法大概耗时0.85左右。但是好处是，对内存的利用比较有效率，这需要你tradeoff了。

用v_outer函数，需要安装java的。但是速度比较快，可以试一下哈！
v_outer函数在qdap包里。
本来我用outer做了一个，速度比for循环要好一些。但这个应该是比较快的。
测了一下10000行，300列的数据。

library(qdap)
dat <- matrix(sample(c('A','T','G','C'),10000*300,replace=TRUE),ncol=300)
newadd <- function(x,y){
    sum(!(x==y))
}
res <- v_outer(t(dat),newadd)

内存16g,处理器i7,6分钟多吧！

jiangbeilu 发表于 2018-4-17 19:53
用v_outer函数，需要安装java的。但是速度比较快，可以试一下哈！
v_outer函数在qdap包里。
本来我用out ...

我找到进一步改进的方法了——
为什么我第一种方法耗内存，其实主要是中间变量太多，每一个20000*20000的中间变量都占用3.2G内存，会使后面的运算越来越不够资源，那么，我可以这样：

1. out <- tcrossprod(dat == "A")
   
2. out <- tcrossprod(dat == "C") + out
   
3. out <- tcrossprod(dat == "G") + out
   
4. out <- tcrossprod(dat == "T") + out
   
5. out <- ncol(dat) - out

复制代码

在我的16G内存上测试20000*300的矩阵成功了，耗时20多秒（这个和我之前的理论推算倒是一致的）；如果内存比较小的，可以试试在上面代码的每一句之间加上：gc()试试。

多谢楼上各位大神的指点。我是在一台八核服务器上跑的，2楼的方法耗时约2个小时左右。
后来我采用了并行的方式，耗时约6分钟可以跑完，以下是我的代码。
有时间再试试楼上各位大神的方法。
初学R，各位大神的指点受益匪浅。

1. gdmatrix<-function(dat){
   
2.   func<-function(x,y){
   
3.     mat<-colSums(x != t(y))
   
4.     return(mat)
   
5.   }
   
6.   cl.cores <- detectCores()
   
7.   cl <- makeCluster(cl.cores-1)
   
8.   aa<-parApply(cl,dat,1,func,y=dat)
   
9.   stopCluster(cl)
   
10.   return(aa)
    
11. }

复制代码