受试者之间和内部的设计与批量校正
1
0
进入编辑模式
demisa•0
@demisa - 22970
最后一次见到是8天前
美国

你好,

我已经搜索了很多关于它的东西,如果我错过了一个解决方案,很抱歉。

我试图在DESeq2中做一个成对的分析,意思是有一个成对的样本设计。

我有一个包含两组的数据集(这是我的批处理),每组都有不均匀的主题数量,每个主题只能在其中一个组中,有一个或两个样本(所以一个或两个条件)。注意,我的计数是用kallisto估算的。只是一个玩具的例子,如何我的colData看起来像:

1 A S1 2样品2 1 B S1 3样品3 1 A S2 4样品4 1 B S2 5样品5 1 A S3 6样品6 1 B S3 7样品7 1 A S4 8样品8 1 B S4 9样品9 1 A S5 10样品10 1 B S5 11样品11 2 A S6 12样品12 2 B S6 13样品13 2 A S7 14样品14 2 B S7 15样品15 2 A S8

在这个例子中,batch == 1有5个受试者,每个受试者都有两个条件,而batch == 2有3个受试者,其中一个只有一个条件。我通过保持平衡的配对样本对条件进行简化,所以我过滤了样本15。

我的目标是在控制主题效果的同时测试条件效果。

所以一开始我认为我的模型应该是~批处理+主体+条件。我要查看的resultName是'condition_B_vs_A',以查看条件效果(同时已经控制了批处理和主题效果)。这种模型设计导致了“模型矩阵不满秩”的误差。

dds = DESeqDataSetFromMatrix(countData = counts.)mat, colData = sample.summary。balanced, design = ~ batch + subject + condition) dds <- DESeq(dds) checkFullRank(modelMatrix)中的Error in checkFullRank(modelMatrix):模型矩阵不是满秩的,因此模型不能按指定的方法拟合。设计公式中的一个或多个变量或交互项是其他变量或交互项的线性组合,必须删除。

问题在于批次和主题的线性关系。小插图中的线性例子并不完全符合我的理解。

虽然我已经尝试了一些基于这里的对话的想法(比如我最后的方法是首先应用批处理校正,转换为正整数,然后使用~ subject + condition运行DESeq2,结果是零deg),但我尝试过的任何方法都不起作用。

顺便说一下,如果我不做一个成对的设计,只是有~batch+条件,模型工作得很好。但我想利用这一点每个科目都有两个条件。

任何见解都将非常感激!

RNASeqDataBatchEffect配对DESeq2•258年的观点
添加评论
0
进入编辑模式
@mikelove
最后一次见到是1天前
美国

这部分:

每个受试者只能在其中一个组中有一个或两个样本(一种或两种情况)

在之前的文章中已经讨论过。在DESeq2中,我们使用了固定效果,所以对pair和条件的固定效果,但limma的duplicateCorrelation方法更灵活,所以你可以研究一下。它将使用所有配对或未配对的样本,并建立主题内部相关性的模型。

在DESeq2中,我建议使用成对的样本,您可以在小插图中使用嵌套在组中的成对样本的策略。

0
进入编辑模式

谢谢Michael的回复。

虽然我已经看到了使用duplicateCorrelation推荐limma的讨论,但我想知道是否可以使用DESeq2来实现。

我的困惑正是在“群体特定条件的影响,个体嵌套在群体中”这部分,因为我想到了主题。n是把两个不同的科目分组,但我现在知道它到底是做什么的了,所以没关系。但我并不是在寻找特定于群体(或特定于批次)的条件效应,而是在调整群体/批次和受试者差异时的条件效应。所以在添加了subject之后,我尝试了下面的模型。n变量:

~ batch + batch:主题。n +条件

但让我再往前走一步,我在简单模型中遇到了一个收敛问题,为此我尝试了两种方法,但都没有解决问题。首先将DESeq函数分解成步骤并增加maxit,其次过滤更多低表达基因。

#对象构造的简化模型:dds = DESeqDataSetFromMatrix(countData = counts.)mat, colData = sample.summary。平衡,设计= ~主体+条件)dds <- DESeq(dds) *估计尺寸因素估计色散基因方面色散估计平均-色散关系最终色散估计拟合模型和测试21行不收敛beta,标记在mcols(对象)$betaConv。所以我重复分成步骤和过滤更多低表达基因:dds = DESeqDataSetFromMatrix(countData = counts.)mat, colData = sample.summary。平衡,设计= ~主题+条件)dds <- estimateSizeFactors(dds) nc <- counts(dds, normalized=TRUE) filter <- rowsum (nc >= 10) >= 4 dds <- dds[filter,] dds <- estimatedispersion (dds) *基因-wise色散估计均值-色散关系最终色散估计* dds <- nbinomWaldTest(dds, maxit=2000) *23行没有收敛在beta,标记在mcols(object)$betaConv。使用更大的maxit参数nbinomWaldTest* #没有解决收敛问题,我尝试向前。#添加组合变量dds ColData: dds$subject。N = ifelse(dds$batch == 1, factor(rep(rep(1:batch1.subjects,each=2),2)), factor(rep(rep(1:batch2.subjects,each=2),2))) dds$subject。N = as.factor(dds$subject. N)样品批次条件subject subject。n 1 1 sample1 S1 2 sample2 1 B S1 1 3 sample3 1 S2 2 4 sample4 1 B S2 2 5 1 sample5 S3 3 1 6 sample6 1 B S3 3 7 sample7 S4 4 8 sample8 1 B S4 4 9 sample9 1 S5 5 10 sample10 1 B S5 5 11 sample11 2 S6 1 12 sample12 2 B S6 1 13 sample13 2 S7 2 14 sample14 2 B S7 2 #饲料新colData DESeqDataSet和编辑设计:dds = DESeqDataSetFromMatrix (countData =。mat, colData = colData(dds),设计= ~ batch + batch*subject。n +条件) *converting counts to integer mode Error in checkFullRank(modelMatrix) : the model matrix is not full rank, so the model cannot be fit as specified. One or more variables or interaction terms in the design formula are linear combinations of the others and must be removed.* # This is because I have unbalanced data in respect to batch: uneven subjects in the two batches. # So I removed the levels without samples, like for example phase2:subject.n3: m = model.matrix(~batch + batch:subject.n + condition, colData(dds)) all.zero <- apply(m, 2, function(x) all(x==0)) length(all.zero[all.zero == TRUE]) idx <- which(all.zero) m <- m[,-idx] # Feed clean model matrix to DESeq: dds = DESeq(dds, full = m) *using supplied model matrix using pre-existing size factors estimating dispersions gene-wise dispersion estimates mean-dispersion relationship final dispersion estimates fitting model and testing 121 rows did not converge in beta, labelled in mcols(object)$betaConv. Use larger maxit argument with nbinomWaldTest* resultsNames(dds) [1] "Intercept" "batch2" "conditionB" [4] "batch1.subject.n2" "batch2.subject.n2" "batch1.subject.n3" [7] "batch1.subject.n4" "batch1.subject.n5"

如何解决收敛问题?我甚至重复了单独的步骤,并添加了useOptim选项:

dds <- estimateSizeFactors(dds) dds <- estimatedistrisions (dds) dds <- nbinomWaldTest(dds, maxit=2000, useOptim = T)

另外,为了得到对批次和受试者差异进行校正的条件效果,是res1还是res2?

res1 = results(dds,反差= c("condition", "B", "A")) res2 = results(dds,反差=list(c("condition_B", "batch1.subject "))n2”、“batch2.subject。n2”、“batch1.subject。n3”、“batch1.subject。陶瓷”、“batch1.subject。它们被”、“batch2”)))
添加回复
0
进入编辑模式

跳到最上面,如果你不想要群体特定的条件效果,那么只需使用~sample +条件。你不能(也不能)控制在比样本更高的分辨率。

0
进入编辑模式

如果我们说的特定群体指的是特定批次(在这种情况下),那么我肯定不这么认为。我只需要修改这批货。

当你说“你不能(也不能)控制在比样本更高的分辨率”时,你是指DESeq2还是一般情况?如果我的模型是~sample + condition,那么我如何考虑sample1和sample2来自同一个主题?或者样品1和样品14是不同批次的?

添加回复
0
进入编辑模式

对不起,上面你有主题.使用~主题+条件.您不需要为嵌套主题的对象添加任何控件。

0
进入编辑模式

嗯,我也有这个想法,但它不是在学科之间和学科内部吗?所以对于between部分,我不需要修改batch吗?当你说“在哪个主题内嵌套”是什么意思?现在是嵌套了吗批处理:subject.n不在这个模型里吗?

还有,还有~主题+条件(这就是我上面提到的模型)有些基因是不会收敛的。

添加回复
1
进入编辑模式

主题是在批内—也许与统计学家讨论这是如何发生的。

关于收敛,我建议如下:

keep <- rowsum (counts(dds) >= 10) >= x dds <- dds[keep,] #然后DESeq()下面

对于选择x来避免有很多0的基因,比如7。

0
进入编辑模式

谢谢你,迈克尔,现在一切都清楚了。

事实上,收敛是通过提高来解决的X = 2/3的样本.谢谢!

添加回复
0
进入编辑模式

如果我可以再问一个问题,为什么即使我达到了收敛,我也没有发现deg ?

dds = DESeqDataSetFromTximport(计数,sample.summary.)平衡,~主题+条件)min.cutoff = round(nrow(sample.summary.balanced)/1.2) keep <- rowsum (counts(dds) >= min.cutoff) >= min.cutoff dds <- dds[keep,] dds <- DESeq(dds) *估计大小因素估计分散基因方面的分散估计平均-分散关系最终分散估计拟合模型和测试* res =结果(dds,对比= c("条件","B", "A"), pAdjustMethod = "BH",alpha = 0.05)摘要(res) *在10735中,非零总读计数调整p值< 0.05 LFC > 0(上):1,0.0093% LFC < 0(下):0,0%异常值[1]:0,0%低计数[2]:0,0%(平均计数< 145)*

得到deg的唯一方法是只通过p值进行过滤。即使是FC。这里是log2FC, padj和p值的范围,以及p值和log2FC分布的火山。

总结(res$log2FoldChange)最小第1曲。中值平均第3曲。最大值。-0.654707 -0.110581 -0.004492 -0.007384 0.097655 0.668882 summary(res$padj)第1个曲数中位数平均第3个曲数最大值0.02358 0.69443 0.77172 0.79575 0.89926 0.99999 summary(res$pvalue)第1个曲数中位数平均第3个曲数最大值0.0000022 0.1736602 0.3860016 0.4302314 0.6745380 0.9999893

用p值绘制火山图

添加回复
0
进入编辑模式

数据集中没有大的变化。你不能拒绝这些基因的null。

您可能需要与统计学家讨论如何增强未来的数据集。

0
进入编辑模式

我完全同意。但当我不做成对分析,而我使用一个简单的~批处理+条件模型时,这种能力非常好。我通过FDR和4FC进行过滤,得到数百个deg。我不知道为什么配对这么差。

添加回复
0
进入编辑模式

抱歉,我在你最初的帖子中错过了这部分。

我猜你可能在100的意义边缘,通过添加所有额外的协变量的主题术语,你失去了自由度。

我想指出的是,考虑到患者的基线,这些LFC非常小。所以当它告诉你变化很小的时候,我不会太强调数据。

0
进入编辑模式

我理解你的观点并同意你的观点。

但当我使用1/10的样本进行配对分析时,我得到了1000个具有p0.01和4FC阈值的deg。不过fdr还是不够好。但我确实认为有更多的样本会增加功率,而不是相反。我觉得很奇怪。但也许你是对的,也许添加所有这些主题词是在最小化任何影响。

你觉得limma会更合适吗?

我当然不想把数据推得太紧,只是想确保我做的方法是正确的。

添加回复
0
进入编辑模式

在软件方面,我没有更多的建议。

0
进入编辑模式

非常感谢您的评论,谢谢

添加回复

登录在添加答案之前。

流量:上个小时访问了253个用户

使用本网站即构成接受我们的用户协议和隐私政策

由的2.3.6版本