1

我有一个大字符矩阵（15000 x 150），格式如下：

       A     B       C     D
[1,] "0/0" "0/1"   "0/0" "1/1"
[2,] "1/1" "1/1"   "0/1" "0/1"
[3,] "1/2" "0/3"   "1/1" "2/2"
[4,] "0/0" "0/0"   "2/2" "0/0"
[5,] "0/0" "0/0"   "0/0" "0/0"

我需要在列之间进行成对比较并获得行的比例

  • 被分隔的字符串都不'/'相等（编码为0）；
  • 只有一个被分隔的字符串'/'相等（编码为1）；
  • 分隔的两个字符串'/'相等（编码为 2）。

上述样本 5 x 4 矩阵的预期输出为

          0    1    2
A    B    0.2  0.2  0.6
A    C    0.2  0.4  0.4  
A    D    0.2  0.4  0.4
B    C    0.4  0.4  0.2
B    D    0.2  0.4  0.4  
C    D    0.6  0.0  0.4

我尝试过使用pmatch，但无法进行成对比较以获得上述输出。任何帮助表示赞赏。

修改后的问题

是否可以排除两对之间的值“0/0”以获得比例？即当比较 A 和 B 时排除 A=B= 0/0 并获得其余部分的比例？

4

3 回答 3

2

这是我到目前为止可以提供的：

fun1 <- function (S) {
  n <- ncol(S)
  ref2 <- combn(colnames(S), 2)
  ref1 <- paste(ref2[1, ], ref2[2, ], sep = "&")
  z <- matrix(0, choose(n, 2), 3L, dimnames = list(ref1, 0:2))
  k <- 1L
  for (j in 1:(n - 1)) {
    x <- scan(text = S[, j], what = integer(), sep = "/", quiet = TRUE)
    for (i in (j + 1):n) {
      y <- scan(text = S[, i], what = integer(), sep = "/", quiet = TRUE)
      count <- tabulate(.colSums(x == y, 2L, length(x) / 2L) + 1L)
      z[k, ] <- count / sum(count)
      k <- k + 1L
      }
    }
  z
  }

scan它看起来很糟糕，因为它有一个用 R 编写的双循环嵌套，但最里面的内核通过使用,.colSums和非常有效tabulate。对于您的 150 列矩阵，迭代的总数是choose(ncol(S), 2)，并不算多。如果你愿意，我可以用fun1Rcpp 版本替换。

## your data
S <- structure(c("0/0", "1/1", "1/2", "0/0", "0/0", "0/1", "1/1", 
"0/3", "0/0", "0/0", "0/0", "0/1", "1/1", "2/2", "0/0", "1/1", 
"0/1", "2/2", "0/0", "0/0"), .Dim = c(5L, 4L), .Dimnames = list(
NULL, c("A", "B", "C", "D")))

fun1(S)
#      0   1   2
#A&B 0.2 0.2 0.6
#A&C 0.2 0.4 0.4
#A&D 0.2 0.4 0.4
#B&C 0.4 0.4 0.2
#B&D 0.2 0.4 0.4
#C&D 0.6 0.0 0.4

表现

哈，当我在 15000 x 150 矩阵上实际测试我的函数时，我发现：

  1. 我可以scan跳出循环嵌套以加快速度，也就是说，我可以一次将字符矩阵扫描成整数矩阵；
  2. scan(text = blabla)需要很长时间，而scan(file = blabla)速度很快，因此值得从文本文件中读取数据；
  3. 使用文本文件对文件的格式很敏感，因此编写健壮的代码很棘手。

我制作了一个fun2具有文件访问权限的版本，以及一个fun3使用 Rcpp 进行循环嵌套的版本。事实证明：

  1. 从文件读取确实更好（但我们必须以“/”分隔格式提供文件）；
  2. 循环的 Rcpp 实现是有益的。

我回来把它们贴在这里（见修订 2），我看到user20650 的开头是strsplit. 当我开始时，我排除strsplit了我的选择，因为我认为使用字符串的操作可能很慢。是的，它很慢，但仍然比scan. 所以我写了一个fun4usingstrsplit和一个对应fun5于 Rcpp （见修订 3）。Profiling 说 60% 的执行时间都花在了，strsplit所以它确实是一个性能杀手。然后我用一个更简单的 C++ 实现替换了strsplit、unlist和。它产生10倍的提升！好吧，如果您仔细考虑，这是有道理的。通过使用C 库中的（或）as.integermatrixatoistrtol<stdlib.h>，我们可以直接把字符串翻译成整数，这样就省去了所有的字符串操作！

长话短说，我只提供最终的、最快的版本。

library(Rcpp)

cppFunction("IntegerMatrix getInt (CharacterMatrix Char) {
  int m = Char.nrow(), n = Char.ncol();
  IntegerMatrix Int(2 * m, n);
  char *s1, *s2;
  int i, *iptr = &Int(0, 0);
  for (i = 0; i < m * n; i++) {
    s1 = (char *)Char[i]; s2 = s1;
    while(*s2 != '/') s2++; *iptr++ = atoi(s1);
    s2++; *iptr++ = atoi(s2);
    }
  return Int;
  }")

cppFunction('NumericMatrix pairwise(NumericMatrix z, IntegerMatrix Int) {
  int m = Int.nrow() / 2, n = Int.ncol();
  int i, j, k, *x, *y, count[3], *end; bool b1 = 0, b2 = 0;
  double M = 1 / (double)m;
  for (k = 0, j = 0; j < (n - 1); j++) {
    end = &Int(2 * m, j);
    for (i = j + 1; i < n; i++, k++) {
      x = &Int(0, j); y = &Int(0, i);
      count[0] = 0; count[1] = 0; count[2] = 0;
      for (; x < end; x += 2, y += 2) {
        b1 = (x[0] == y[0]);
        b2 = (x[1] == y[1]);
        count[(int)b1 + (int)b2]++;
        }
      z(k, 0) = (double)count[0] * M;
      z(k, 1) = (double)count[1] * M;
      z(k, 2) = (double)count[2] * M;
      }
    }
  return z;
  }')

fun7 <- function (S) {
  ## separate rows using Rcpp; `Int` is an integer matrix
  n <- ncol(S)
  Int <- getInt(S)
  m <- nrow(Int) / 2
  ## initialize the resulting matrix `z`
  ref2 <- combn(colnames(S), 2)
  ref1 <- paste(ref2[1, ], ref2[2, ], sep = "&")
  z <- matrix(0, choose(n, 2), 3L, dimnames = list(ref1, 0:2))
  ## use Rcpp for pairwise summary
  pairwise(z, Int)
  }

让我们生成一个随机的 15000 x 150 矩阵并尝试一下。

sim <- function (m, n) {
  matrix(sample(c("0/0", "0/1", "1/0", "1/1"), m * n, TRUE), m, n,
         dimnames = list(NULL, 1:n))
  }

S <- sim(15000, 150)
system.time(oo <- fun7(S))
#   user  system elapsed 
#  1.324   0.000   1.325

哦，这是闪电快！

是否可以排除两对之间的值“0/0”以获得比例？即当比较 A 和 B 时排除 A=B= 0/0 并获得其余部分的比例？

这种适应在 C/C++ 级别上很简单。只是一个附加if测试。

## a new C++ function `pairwise_exclude00`
cppFunction('NumericMatrix pairwise_exclude00(NumericMatrix z, IntegerMatrix Int) {
  int m = Int.nrow() / 2, n = Int.ncol();
  int i, j, k, *x, *y, count[3], size, *end;
  bool b1 = 0, b2 = 0, exclude = 0;
  double M; 
  for (k = 0, j = 0; j < (n - 1); j++) {
    end = &Int(2 * m, j);
    for (i = j + 1; i < n; i++, k++) {
      x = &Int(0, j); y = &Int(0, i);
      count[0] = 0; count[1] = 0; count[2] = 0; size = 0;
      for (; x < end; x += 2, y += 2) {
        b1 = (x[0] == y[0]);
        b2 = (x[1] == y[1]);
        exclude = (x[0] == 0) & (x[1] == 0) & b1 & b2;
        if (!exclude) {
          count[(int)b1 + (int)b2]++;
          size++;
          }
        }
      M = 1 / (double)size;
      z(k, 0) = (double)count[0] * M;
      z(k, 1) = (double)count[1] * M;
      z(k, 2) = (double)count[2] * M;
      }
    }
  return z;
  }')

## re-define `fun7` with a new logical argument `exclude00`
fun7 <- function (S, exclude00) {
  ## separate rows using Rcpp; `Int` is an integer matrix
  n <- ncol(S)
  Int <- getInt(S)
  m <- nrow(Int) / 2
  ## initialize the resulting matrix `z`
  ref2 <- combn(colnames(S), 2)
  ref1 <- paste(ref2[1, ], ref2[2, ], sep = "&")
  z <- matrix(0, choose(n, 2), 3L, dimnames = list(ref1, 0:2))
  ## use Rcpp for pairwise summary
  if (exclude00) pairwise_exclude00(z, Int)
  else pairwise(z, Int)
  }

使用S您问题中的示例：

fun7(S, TRUE)
#            0         1         2
#A&B 0.3333333 0.3333333 0.3333333
#A&C 0.3333333 0.6666667 0.0000000
#A&D 0.3333333 0.6666667 0.0000000
#B&C 0.5000000 0.5000000 0.0000000
#B&D 0.3333333 0.6666667 0.0000000
#C&D 0.7500000 0.0000000 0.2500000
于 2018-08-18T00:11:35.883 回答
1

这使用了来自李哲源的回答的想法，特别是tabulate-- 加快了速度。在我的旧笔记本电脑上，15000x160 的数据需要大约 14 秒

# split strings and form matrix for each column
ap =  matrix(unlist(strsplit(m, "/")), nc=2, byrow=TRUE)
ap = split.data.frame(ap, rep(colnames(m), each=nrow(m))) # maybe a way to use array?

# get 2-way combination of column names
co = combn(colnames(m), 2)

# test equality of each matrix
ap = apply(co, 2, function(x) tabulate(rowSums(ap[[x[1]]]==ap[[x[2]]])+1, 3))

# output
data.frame(t(co), t(ap)/nrow(m))

数据

m = as.matrix(read.table(header=T, text='       A     B       C     D
 "0/0" "0/1"   "0/0" "1/1"
 "1/1" "1/1"   "0/1" "0/1"
 "1/2" "0/3"   "1/1" "2/2"
 "0/0" "0/0"   "2/2" "0/0"
 "0/0" "0/0"   "0/0" "0/0"'))

m = do.call(cbind, replicate(40 , m, simplify = FALSE))
m = do.call(rbind, replicate(3000, m, simplify = FALSE))
colnames(m) =  paste0("A", 1:160)
于 2018-08-18T01:56:09.473 回答
0

您可以创建 3 个函数来指示 0、1、2 个条件，然后遍历列名以获得不同的对并应用函数来创建结果 data.frame：

library(tidyr)
matrix <- read.csv("matrix.csv", stringsAsFactors = F)
n <-nrow(matrix)
c <- ncol(matrix)
zero <- function(A, B){ res <- sum(!grepl("0", A) & !grepl("0", B))/n }
one <- function(A, B) {
  A <- unlist(str_split(A, "/"))
  B <- unlist(str_split(B, "/"))
  comp <-data.frame(cbind(A==B, c(1,2), id= sort(rep(1:n,2))))%>%spread(V2, V1)
  res <- sum(sum(comp[,2]+comp[,3])>0)/n} 
two <- function(A, B){res <- sum(A==B)/n}  

res <-data.frame()
k <-1
for (i in 1:(c-1)){
  for (j in (i+1):c){
    A<-matrix[,i]
    B<-matrix[,j]
    res[k,1] <- colnames(matrix)[i]
    res[k,2] <- colnames(matrix)[j]
    res[k,3] <- zero(A,B)
    res[k,4] <- one(A,B)
    res[k,5] <- two(A,B)
    k <-k+1
  }
}
colnames(res) <-c("G1", "G2", "0", "1", "2")
于 2018-08-17T23:57:54.897 回答