考虑到我有两个向量，一个叫residues，第二个叫scores，它们有31个分数，每个残基一个，都是正数。为了说明，得到两个向量如下图：
residues <- 1:31
scores <- runif(n = 31, min = 0.35, max = 3.54)
我正在考虑一个随机序列来举例说明。
如果我绘制 residues x scores 我将得到以下图形：
我想要做的是：我将考虑 15 个残基（以下称为 15mer）的特定组合，跳过一个残基（即 1:15、2:16、3:17 一直到 17 个： 31) 我想计算所有这 17 种组合的曲线下面积 (AUC)。我的最终目标是选择具有最高 AUC 的 15mer。
可以使用 zoo 包中的 rollmean 函数计算 AUC，如this question 所示。但是，正如我在此示例中所拥有的 17 种可能的组合一样，我正在尝试找到一个脚本来自动化该过程。 提前致谢。
【问题讨论】：
这是我使用过的以下功能
scores <- runif(n = 31, min = 0.35, max = 3.54)
fun <- function(dat, n) {
require(zoo)
N <- which(max(rollmean(dat, n)) == rollmean(dat, n))
output <- matrix(0, length(N), n)
for (i in 1:length(N)) {
output[i, ] <- dat[N[i]:(N[i] + n - 1)]
}
output
}
fun(scores, 15)
让我们从里到外跑吧
rollmean(dat, n)
您提到的 zoo 包为我们提供了滚动平均值，我们
max(rollmean(dat, n))
找到滚动平均值的最大值
max(rollmean(dat, n)) == rollmean(dat, n)
返回一个 TRUE/FALSE 向量，其中滚动均值等于最大值
N <- which(max(rollmean(dat, n)) == rollmean(dat, n))
返回最大值的索引。根据您的数据，您可能有多个序列获得最大值，我们决定使用以下循环返回所有序列
for (i in 1:length(N)) {
output[i, ] <- dat[N[i]:(N[i] + n -1)]
}
结果：
set.seed(12345)
scores <- runif(n = 31, min = 0.35, max = 3.54)
fun(scores, 15)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.588179 1.633928 0.9208938 3.385791 1.797393 1.39234
[,7] [,8] [,9] [,10] [,11] [,12]
[1,] 3.429675 2.606867 2.406091 1.593553 2.578354 2.085545
[,13] [,14] [,15]
[1,] 1.07243 1.895739 2.879693
fun(rpois(1000, 1), 10)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 1 4 2 1 1 3 3 2 2
[2,] 1 4 2 1 1 3 3 2 2 1
[3,] 4 2 1 1 3 3 2 2 1 1
【讨论】：
library(zoo)
set.seed(555)
residues <- 1:31
scores <- runif(n = 31, min = 0.35, max = 3.54)
which.max(sapply(1:17, function(x){sum(diff(residues[x:(x+14)])*rollmean(scores[x:(x+14)],2))}))
# result 7 i.e. 7:21
或
sapply(1:17, function(x){sum(diff(residues[x:(x+14)])*rollmean(scores[x:(x+14)],2))}) # gives you the AUCs
# result [1] 28.52530 29.10203 28.52847 27.65325 27.19925 28.77782 29.29373 28.13133 28.23705 27.68724 25.75294 25.27226 25.44963 25.81201 25.49907 23.48632
#[17] 22.45763
或使用自定义函数
f_AUC <- function(x, y, lngth){
sapply(1:(length(x)-lngth+1), function(z) sum(diff(x[z:(z+lngth-1)])*rollmean(y[z:(z+lngth-1)],2)))
}
f_AUC(x=residues, y=scores, lngth=15)
【讨论】：