哪个没有按预期工作

Question

我有一个包含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

标准化分数为:

得分本身除以得分的sqrt,其First.Protein和Second.Protein指数均为4乘以其First.Protein和Second.Protein指数均为3的得分的sqrt.

因此:

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])]))
}

Answer 1

这是对原始尝试的重写(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]之前和他们已经被更新后出现在右侧.这是一个修订版,原始列被解释为"只读"

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)来避免担心返回值的分配

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,]

这是原始数据框的第i个分数

s_ij = df$Score

查找First.Protein的df中ii并记录成绩; 同样为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()而不是多次调用which()for循环内部或测试内部的标识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()的输入列以与f2()的输入列相同的顺序得到归一化结果,因此无法与其他列表进行比较(尽管它们看起来正确'通过眼睛').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()只是查找"对角线"元素

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])
}