R 中的哪个函数没有给出所需的输出

The which function in R is not giving the desired output

我有一个包含 3 列且总共 10,000 个元素的矩阵。第一列和第二列是索引,第三列是分数。我想根据这个公式规范化分数列:

Normalized_score_i_j = score_i_j / ((sqrt(score_i_i) * (sqrt(score_j_j))

score_i_j = 当前分数本身

score_i_i = 查看第一列中当前分数的索引,并在数据集中查找第一列和第二列中均具有该索引的分数

score_j_j = 在第二列中查看当前分数的索引,并在数据集中查找在第一列和第二列中均具有该索引的分数

举个例子,如果df如下:

df <- read.table(text = "
First.Protein,Second.Protein,Score
1,1,25
1,2,90
1,3,82
1,4,19
2,1,90
2,2,99
2,3,76
2,4,79
3,1,82
3,2,76
3,3,91
3,4,33
4,1,28
4,2,11
4,3,99
4,4,50
", header = TRUE, sep = ",")

如果我们规范化这一行:

First.Protein Second.Protein Score
4             3              99

标准化分数将为:

分数本身除以其 First.Protein 和 Second.Protein 索引均为 4 的分数的平方根乘以其 First.Protein 和 [=39] 的分数的平方根=] 索引都是 3.

因此:

Normalized =  99 / (sqrt(50) * sqrt(91)) = 1.467674

我有下面的代码,但它的行为很奇怪,给我的值根本没有规范化,实际上很奇怪:

for(i in 1:nrow(Smith_Waterman_Scores))
{
  Smith_Waterman_Scores$Score[i] <- 
    Smith_Waterman_Scores$Score[i] / 
    (sqrt(Smith_Waterman_Scores$Score[which(Smith_Waterman_Scores$First.Protein==Smith_Waterman_Scores$First.Protein[i] & Smith_Waterman_Scores$Second.Protein==Smith_Waterman_Scores$First.Protein[i])])) *
    (sqrt(Smith_Waterman_Scores$Score[which(Smith_Waterman_Scores$First.Protein==Smith_Waterman_Scores$Second.Protein[i] & Smith_Waterman_Scores$Second.Protein==Smith_Waterman_Scores$Second.Protein[i])]))
}

你做这件事的方式可能很迂回。你能看看这对你有用吗:

R> xx
    First Second Score
1      1      1    25
2      1      2    90
3      1      3    82
4      1      4    19
5      2      1    90
6      2      2    99
7      2      3    76
8      2      4    79
9      3      1    82
10     3      2    76
11     3      3    91
12     3      4    33
13     4      1    28
14     4      2    11
15     4      3    99
16     4      4    50
R> contingency = xtabs(Score ~ ., data=xx)
R> contingency
    Second
First  1  2  3  4
    1 25 90 82 19
    2 90 99 76 79
    3 82 76 91 33
    4 28 11 99 50
R> diagonals <- unname(diag(contingency))
R> diagonals
[1] 25 99 91 50

R> normalize <- function (i, j, contingencies, diagonals) {
+      contingencies[i, j] / (sqrt(diagonals[i]) * sqrt(diagonals[j]))
+  }

R> normalize(4, 3, contingency, diagonals)
[1] 1.467674

使用应用遍历行:

#compute
df$ScoreNorm <- 
  apply(df, 1, function(i){
    i[3] /
      (
        sqrt(df[ df$First.Protein == i[1] &
                   df$Second.Protein == i[1], "Score"]) *
          sqrt(df[ df$First.Protein == i[2] &
                     df$Second.Protein == i[2], "Score"])
      )
  })

#test output
df[15, ]
#    First.Protein Second.Protein Score ScoreNorm
# 15             4              3    99  1.467674

您可以通过连接来实现,这里是一个使用 data.table:

的例子
library(data.table)
dt <- data.table(df)

dt.lookup <- dt[First.Protein == Second.Protein]
setkey(dt,"First.Protein" )
setkey(dt.lookup,"First.Protein" )
colnames(dt.lookup) <- c("First.Protein","Second.Protein","Score1")
dt <- dt[dt.lookup]
setkey(dt,"Second.Protein" )
setkey(dt.lookup,"Second.Protein")
colnames(dt.lookup) <- c("First.Protein","Second.Protein","Score2")
dt <- dt[dt.lookup][
   , Normalized :=  Score / (sqrt(Score1) * sqrt(Score2))][
  , .(First.Protein, Second.Protein, Normalized)]

请确保您不使用 for 循环。

这里是你最初尝试的重写(which() 不是必需的;只需使用逻辑向量进行子设置;with() 允许你在数据框中引用变量而无需必须重新输入 data.frame 的名称 -- 更容易阅读但也更容易出错)

orig0 <- function(df) {
    for(i in 1:nrow(df)) {
        df$Score[i] <- with(df, {
            ii <- First.Protein == First.Protein[i] &
                Second.Protein == First.Protein[i]
            jj <- First.Protein == Second.Protein[i] &
                Second.Protein == Second.Protein[i]
            Score[i] / (sqrt(Score[ii]) * sqrt(Score[jj]))
        })
    }
    df$Score
}

问题是 Score[ii]Score[jj] 在更新前后都出现在右侧。这是一个修订版,其中原始列被解释为 'read-only'

orig1 <- function(df) {
    normalized <- numeric(nrow(df))     # pre-allocate
    for(i in 1:nrow(df)) {
        normalized[i] <- with(df, {
            ii <- First.Protein == First.Protein[i] &
                Second.Protein == First.Protein[i]
            jj <- First.Protein == Second.Protein[i] &
                Second.Protein == Second.Protein[i]
            Score[i] / (sqrt(Score[ii]) * sqrt(Score[jj]))
        })
    }
    normalized
}

我认为现在的结果是正确的(见下文)。更好的实现是使用 sapply(或 vapply)来避免担心 return 值

的分配
orig2 <- function(df) {
    sapply(seq_len(nrow(df)), function(i) {
        with(df, {
            ii <- First.Protein == First.Protein[i] &
                Second.Protein == First.Protein[i]
            jj <- First.Protein == Second.Protein[i] &
                Second.Protein == Second.Protein[i]
            Score[i] / (sqrt(Score[ii]) * sqrt(Score[jj]))
        })
    })
}

现在结果是正确的,我们可以询问性能。您的解决方案需要每次通过循环扫描,例如 First.Protein。 First.Protein 有 N=nrow(df) 个元素,你要经过 N 次循环,所以你将进行 N * N = N^2 次比较——如果你增加数据框的大小从 10 行到 100 行,所用时间将从 10 * 10 = 100 个单位变为 100 * 100 = 10000 个时间单位。

一些答案试图避免多项式缩放。我的答案是在值向量上使用 match() 来做到这一点;这可能缩放为 N(每次查找都在恒定时间内发生,并且有 N 次查找),这比多项式要好得多。

创建具有相同第一和第二蛋白质的数据子集

ii = df[df$First.Protein == df$Second.Protein,]

这是原始数据框中的第 ij 个分数

s_ij = df$Score

ii中查找df的First.Protein并记录得分;同样 Second.Protein

s_ii = ii[match(df$First.Protein, ii$First.Protein), "Score"]
s_jj = ii[match(df$Second.Protein, ii$Second.Protein), "Score"]

标准化后的分数是

> s_ij / (sqrt(s_ii) * sqrt(s_jj))
 [1] 1.0000000 1.8090681 1.7191871 0.5374012 1.8090681 1.0000000 0.8007101
 [8] 1.1228571 1.7191871 0.8007101 1.0000000 0.4892245 0.7919596 0.1563472
[15] 1.4676736 1.0000000

这会很快,使用对 match() 的单个调用,而不是在 for 循环内多次调用 which() 或在 apply() 内测试身份——两者后者进行 N^2 次比较,因此缩放比例很差。

我将一些建议的解决方案总结为

f0 <- function(df) {
    contingency = xtabs(Score ~ ., df)
    diagonals <- unname(diag(contingency))
    i <- df$First.Protein
    j <- df$Second.Protein
    idx <- matrix(c(i, j), ncol=2)
    contingency[idx] / (sqrt(diagonals[i]) * sqrt(diagonals[j]))
}

f1 <- function(df) {
    ii = df[df$First.Protein == df$Second.Protein,]
    s_ij = df$Score
    s_ii = ii[match(df$First.Protein, ii$First.Protein), "Score"]
    s_jj = ii[match(df$Second.Protein, ii$Second.Protein), "Score"]
    s_ij / (sqrt(s_ii) * sqrt(s_jj))
}

f2 <- function(dt) {
    dt.lookup <- dt[First.Protein == Second.Protein]
    setkey(dt,"First.Protein" )
    setkey(dt.lookup,"First.Protein" )
    colnames(dt.lookup) <- c("First.Protein","Second.Protein","Score1")
    dt <- dt[dt.lookup]
    setkey(dt,"Second.Protein" )
    setkey(dt.lookup,"Second.Protein")
    colnames(dt.lookup) <- c("First.Protein","Second.Protein","Score2")
    dt[dt.lookup][
      , Normalized :=  Score / (sqrt(Score1) * sqrt(Score2))][
      , .(First.Protein, Second.Protein, Normalized)]
}

f3 <- function(dt) {
    eq = dt[First.Protein == Second.Protein]
    dt[eq, Score_ii := i.Score, on = "First.Protein"]
    dt[eq, Score_jj := i.Score, on = "Second.Protein"]
    dt[, Normalised := Score/sqrt(Score_ii * Score_jj)]
    dt[, c("Score_ii", "Score_jj") := NULL]
}

我知道如何以编程方式检查前两个生成的结果是否一致;我不太了解 data.table 以与 f2() 的输入列相同的顺序获得规范化结果,因此无法与其他列进行比较(尽管它们看起来正确 'by eye')。 f3() 产生数值相似但不相同的结果

> identical(orig1(df), f0(df))
[1] TRUE
> identical(f0(df), f1(df))
[1] TRUE
> identical(f0(df), { f3(dt3); dt3[["Normalized"]] })  # pass by reference!
[1] FALSE
> all.equal(f0(df), { f3(dt3); dt3[["Normalized"]] })
[1] TRUE

存在性能差异

library(data.table)    
dt2 <- as.data.table(df)
dt3 <- as.data.table(df)

library(microbenchmark)
microbenchmark(f0(df), f1(df), f2(dt2), f3(dt3))

> microbenchmark(f0(df), f1(df), f2(df), f3(df))
Unit: microseconds
   expr      min        lq      mean    median       uq      max neval
 f0(df)  967.117  992.8365 1059.7076 1030.9710 1094.247 2384.360   100
 f1(df)  176.238  192.8610  210.4059  207.8865  219.687  333.260   100
 f2(df) 4884.922 4947.6650 5156.0985 5017.1785 5142.498 6785.975   100
 f3(df) 3281.185 3329.4440 3463.8073 3366.3825 3443.400 5144.430   100

解决方案 f0 - f3 很可能适用于真实数据(尤其是 data.table);时间以微秒为单位的事实可能意味着速度并不重要(现在我们没有实现 N^2 算法)。

经过深思熟虑,f1() 的更直接实现只是查找 'diagonal' 元素

f1a <- function(df) {
    ii = df[df$First.Protein == df$Second.Protein, ]
    d = sqrt(ii$Score[order(ii$First.Protein)])
    df$Score / (d[df$First.Protein] * d[df$Second.Protein])
}    

以下是我使用 data.table 的方法。希望@MartinMorgan 发现这更容易理解:-)。

require(data.table) # v1.9.6+
dt = as.data.table(df) # or use setDT(df) to convert by reference
eq = dt[First.Protein == Second.Protein]

到目前为止,我刚刚创建了一个新的 data.table eq,其中包含两列相等的所有行。

dt[eq, Score_ii := i.Score, on = "First.Protein"]
dt[eq, Score_jj := i.Score, on = "Second.Protein"]

我们在这里 添加Score_iiScore_jj 同时加入 First.ProteinSecond.Protein。由于 on= 参数,应该清楚这是一个连接操作。 i. 指的是 i- 参数中提供的 data.table 中的 Score 列(这里是 eqScore)。

Note that we can use match() here as well. But that wouldn't work if you've to lookup directly (and as efficiently) based on more than one column. Using on=, we can extend this quite easily, and is also much easier to read/understand.

一旦我们获得了所有必需的列,任务就是获取最终的 Normalised 列(如果不需要,则删除中间列)。

dt[, Normalised := Score/sqrt(Score_ii * Score_jj)]
dt[, c("Score_ii", "Score_jj") := NULL] # delete if you don't want them

我将省略微秒和毫秒基准,因为我对它们不感兴趣。


PS:Score_iiScore_jj 列是在假设您可能需要它们的情况下有意添加的。如果你根本不需要它们,你也可以这样做:

Score_ii = eq[dt, Score, on = "First.Protein"] ## -- (1)
Score_jj = eq[dt, Score, on = "Second.Protein"]

(1) 读取:对于 dt 中的每一行,获取 eq 中的匹配行,同时匹配列 First.Protein 并提取与该匹配行对应的 eq$Score

然后,我们可以直接添加Normalised列为:

dt[, Normalised := Score / sqrt(Score_ii * Score_jj)]