2

我正在寻找一种方法来识别在各种条件下显着稳定的基因。换句话说，与标准 DE 分析相反。

标准 DE 将基因分为两类：一侧显着变化，其他一切，“其余”，另一侧。然而，“其余的”既包含实际上不改变的基因，也包含对改变的置信度不足以称它们为差异基因的基因。
我想要的是找到那些不会改变的，或者换句话说，那些我可以自信地说我的条件没有改变的。

我知道这在 DEseq 中是可能的，方法是提供一个替代的零假设，但我必须将它作为一个额外的步骤集成到已经使用 limma 的其他人的管道中，我想坚持下去。理想情况下，我想以类似的方式测试 DE 和非改变基因，这在概念上类似于改变 DEseq 中的 H0。

目前，测试 DE 的代码如下：

# shaping data
comparison <- eBayes(lmFit(my_data, weights = my.weights^2))
results <- limma::topTable(my_data, sort.by = "t",  
                     coef = 1, number = Inf)

作为一个例子，我喜欢下面的东西，但任何概念上相似的东西都可以。

comparison <- eBayes(lmFit(my_data, weights = my.weights^2), ALTERNATIVE_H0 = my_H0)

我知道treat（）允许通过提供折叠更改来指定区间零假设，引用手册：“它使用区间零假设，其中区间为[-lfc，lfc]”。
但是，这仍然会测试从 0 左右的中心间隔开始的变化，而我想测试的间隔是 [-inf,-lfc] + [lfc,inf]。

我有什么选择吗？

谢谢！

4

1 回答 1

2

您可以尝试使用 logFC 的置信区间来选择您的基因，但我必须说这在很大程度上取决于您拥有的样本数量，以及生物学方差的强度。下面我展示了一个例子：

首先我们使用 DESeq2 生成一个示例数据集，我们设置 betaSD 以便我们有一小部分基因应该显示条件之间的差异

library(DESeq2)
library(limma)
set.seed(100)
dds = makeExampleDESeqDataSet(n=2000,betaSD=1)
#pull out the data
DF = colData(dds)
# get out the true fold change
FC = mcols(dds)

现在我们可以在这个数据集上运行 limma-voom，

V = voom(counts(dds),model.matrix(~condition,data=DF))
fit = lmFit(V,model.matrix(~condition,data=DF))
fit = eBayes(fit)
# get the results, in this case, we are interested in the 2nd coef
res = topTable(fit,coef=2,n=nrow(V),confint=TRUE)

所以有一个选项可以在函数 topTable 中收集倍数变化的 95% 置信区间。我们这样做并与真正的 FC 进行比较：

# fill in the true fold change
res$true_FC = FC[rownames(res),"trueBeta"]

我们可以看看估计值和真实值有何不同：

plot(res$logFC,res$true_FC)

在此处输入图像描述

因此，假设我们想要找到基因，我们确信其中的倍数变化 < 1，我们可以这样做：

tabResults = function(tab,fc_cutoff){
true_unchange = abs(tab$true_FC)<fc_cutoff
pred_unchange = tab$CI.L>(-fc_cutoff) & res$CI.R <fc_cutoff
list(
X = table(pred_unchange,true_unchange),
expression_distr = aggregate(
tab$AveExpr ~ pred_unchange+true_unchange,data=tab,mean
))
}

tabResults(res,1)$X

             true_unchange
pred_unchange FALSE TRUE
        FALSE   617 1249
        TRUE      7  127

上述结果告诉我们，如果我们将其限制为 95% 置信度在 +/- 1 FC 内的基因，我们会得到 134 个命中，其中 7 个是错误的（实际倍数变化 > 1）。

我们错过了一些真正不变的基因的原因是因为它们的表达略低，而我们正确预测为不变的大部分基因却具有高表达：

tabResults(res,1)$expression_distr

  pred_unchange true_unchange tab$AveExpr
1         FALSE         FALSE    7.102364
2          TRUE         FALSE    8.737670
3         FALSE          TRUE    6.867615
4          TRUE          TRUE   10.042866

我们可以降低 FC，但我们最终也会减少基因：

tabResults(res,0.7)

             true_unchange
pred_unchange FALSE TRUE
        FALSE   964 1016
        TRUE      1   19

置信区间很大程度上取决于您拥有的样本数量。因此，一个数据集的截止值 1 对另一个数据集意味着不同的东西。

所以我想说，如果你手头有一个数据集，你可以先在数据集上运行 DESeq2，获得平均方差关系并像我一样模拟数据，或多或少猜测，什么折叠变化截止就可以了，多少你能得到，并从那里做出决定。

于 2019-12-10T01:28:44.567 回答