1

我正在尝试使用Bioconductor 的 sva 包应用代理变量分析。小插图中的示例工作正常，但是当我尝试使用真实数据时，出现“下标越界”错误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 与我的训练集相同。）这导致 aV为 100 x 0；“下标越界”错误是因为irwsva.build尝试索引到V. 我的数据一定有某些东西导致了这种行为——但是什么？

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

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

这有效，但我需要随后调用fsva. 那失败了，因为调用svawithmethod="two-step"导致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得足以说出为什么会发生这种情况？我能做些什么来让它在我的数据上工作吗？

4

2 回答 2

3

失败的最接近的原因是irwa.build使用快速奇异值分解，它只返回矩阵的正奇异值，如?fast.svd. 在您的数据中，唯一的值是零，这不是正数，因此您必须使用 plainsvd而不是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]

您可以在此处获取代码：

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

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

分解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 值，这就是它不返回错误的原因。

于 2012-07-04T22:06:34.193 回答
0

来自Bioconductor 邮件列表的 John Leek（作者sva）的回复：

这个问题可能是因为您正在考虑的基因/特征数量较少（453）和响应变量的高维度（12）。由于响应变量有这么多不同的水平，许多特征可能与响应显着相关。sva 算法中的部分迭代是降低与响应强相关的特征的权重，因此整个数据集的权重被降低到 0。

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

于 2012-07-05T19:00:53.830 回答