代理变量分析因“下标越界”而失败

Question

我正在尝试使用Bioconductor's sva package 应用代理变量分析。 the vignette 中的示例运行良好，但是当我尝试使用真实数据时，irwsva.build 中出现“下标越界”错误：

$ R

R version 2.15.0 (2012-03-30)
…
> trainData <- read.table('http://www.broadinstitute.org/~ljosa/svaproblem/trainData.txt')
> trainpheno <- read.table('http://www.broadinstitute.org/~ljosa/svaproblem/trainpheno.txt')
> testData <- read.table('http://www.broadinstitute.org/~ljosa/svaproblem/testData.txt')
> trainData <- as.matrix(trainData)
> testData <- as.matrix(testData)
> library(sva)
> trainMod <- model.matrix(~as.factor(label), trainpheno)
> num.sv(trainData, trainMod)
[1] 8
> trainMod0 <- model.matrix(~1, trainpheno)
> trainSv <- sva(trainData, trainMod, trainMod0)
Number of significant surrogate variables is:  8 
Iteration (out of 5 ):1  2  3  4  5  Error in irwsva.build(dat = dat, mod = mod, mod0 = mod0, n.sv = n.sv,  : 
  subscript out of bounds

使用debug() 缩小范围的尝试显示fast.svd 正在被一个453 x 100 的全零矩阵调用。 （尺寸 453 x 100 与我的训练集相同。）这导致 V 为 100 x 0； “下标越界”错误是因为 irwsva.build 试图索引到 V。我的数据一定有某些东西导致了这种行为——但是什么？

作为一种可能的解决方法，我尝试使用method="two-step" 调用sva：

> trainSv <- sva(trainData, trainMod, trainMod0, method='two-step')
Number of significant surrogate variables is:  8 

这行得通，但我需要随后致电fsva。失败是因为用method="two-step" 调用sva 导致trainSv$pprob.b 为NULL。

那么我的数据与小插图中的数据有何不同？在这两种情况下，训练和测试数据都是矩阵。在小插图中，训练矩阵为 22283 x 30；在我的例子中，它是 453 x 100。在小插图中，感兴趣的变量 (cancer) 是二进制的；在我的例子中，因变量可以取 12 个不同的值。

最后一个区别似乎很重要，因为如果我将范围缩小到 [0, 7]，它会起作用：

> trainMod <- model.matrix(~as.factor(label), trainpheno %% 8)
> trainSv <- sva(trainData, trainMod, trainMod0)
Number of significant surrogate variables is:  9 
Iteration (out of 5 ):1  2  3  4  5  > 

考虑到 100 个样本（列）对于 12 个类来说可能是不够的，我尝试了一个包含 293 个列的类似数据集。 （数据来自同一个实验，但分析的是 293 个单独的样本而不是 100 个处理。）它没有帮助：

> trainData <- read.table('http://www.broadinstitute.org/~ljosa/svaproblem/trainData3.txt')
> trainpheno <- read.table('http://www.broadinstitute.org/~ljosa/svaproblem/trainpheno.txt')
> trainData <- as.matrix(trainData)
> trainMod <- model.matrix(~as.factor(label), trainpheno)
> trainMod0 <- model.matrix(~1, trainpheno)
> trainSv <- sva(trainData, trainMod, trainMod0)
Number of significant surrogate variables is:  11 
Iteration (out of 5 ):1  2  3  4  5  Error in irwsva.build(dat = dat, mod = mod, mod0 = mod0, n.sv = n.sv,  : 
  subscript out of bounds

如果我将 sva 限制为一次迭代，它可以运行完成，但我不知道我是否可以相信结果：

> trainSv <- sva(trainData, trainMod, trainMod0, B=1)
Number of significant surrogate variables is:  11 
Iteration (out of 1 ):1  > 

有人理解irwsva 足以说明为什么会这样吗？我可以做些什么来让它在我的数据上运行？

Answer 1

失败的近端原因是irwa.build 使用快速奇异值分解，它只返回矩阵的正 奇异值，如?fast.svd 中所述。在您的数据中，唯一的值是零，这不是正数，因此您必须使用普通的 svd 而不是 fast.svd。

我创建了一个修补函数sva.patched，它稍微修补了irwa.build 和sva 函数来处理这种外部情况。我基本上在irwa.build中更改了一行：

# Before
sv = fast.svd(dats, tol = 0)$v[, 1:n.sv]
# After
if(any(dats!=0)) sv = fast.svd(dats, tol = 0)$v[, 1:n.sv]
else sv=svd(dats)$v[, 1:n.sv]

您可以领取代码here：

但真正的问题是，为什么这些数据最终会产生一个零值矩阵？我对这种方法了解不多，但我可以给你一些线索。

据我所知，您正确使用了这些功能。但是，如果您检查循环 irwsva.build 函数，您会发现如果 edge.ldfr 函数返回 0，它将返回一个零矩阵。该函数仅在 f.pvalue 没有返回 p 值时返回零高于 0.8。

分解irwa.build，这是从您的数据开始的：

dat=trainData
mod=trainMod
mod0=trainMod0
Id <- diag(ncol(dat))
resid <- dat %*% (Id - mod %*% solve(t(mod) %*% mod) %*% t(mod))
uu <- eigen(t(resid) %*% resid)
# Iterations begin.
mod.b <- cbind(mod, uu$vectors[, 1:n.sv])
mod0.b <- cbind(mod0, uu$vectors[, 1:n.sv])
ptmp <- f.pvalue(dat, mod.b, mod0.b)
which(ptmp>0.8)
# Only one value

现在，第一次循环时，只有一个 p 值高于 0.8。到第二次迭代时，没有，这是所有零的原因。

如果你在小插图数据上运行相同的代码，你会发现它有很多高于 0.8 的 p 值，这就是它不返回错误的原因。

Answer 2

来自 John Leek（sva 的作者）on the Bioconductor mailing list 的回复：

这个问题很可能是因为基因/特征数量少 您正在考虑 (453) 和响应的高维度 变量 (12)。有这么多不同级别的响应变量， 许多特征可能与响应显着相关。 sva 算法中的部分迭代是降低特征的权重 与响应密切相关，因此整个数据集正在 权重为 0。

我建议只运行一次 sva 迭代。通常需要一个 收敛的迭代次数非常少，并且由于您的数据是 特征数量的维数相对较低，这可能是 如果您正在进行工件发现，您可以做的最好的事情。