RNAseq 数据的广义线性模型或名义 Logistic同事,
想象一下,我有在一次测序实验中在 2 个环境中测量的读数计数形式的 RNAseq 数据,每个环境有 3 个重复。 像这样:
环境读取
一个 9
一个 10
一个 11
乙4
乙5
乙6
当然,我可以将每个读取计数视为一个自变量并在测试时执行,我的功效将非常低。 相反,我希望能够将每次读取视为一次观察,然后注意到在一次给定的测序工作中,我在来自环境 B 的样本中获得了 15 个读取,但在来自环境 A 的样本中获得了 30 个读取。 这不太可能是偶然观察到的,与 t 检验不同,我检测环境 A 和 B 之间差异的能力要高得多,并且会随着测序工作的进行而增加。 (在典型的 RAseq 实验中,读数数量要高得多,对于某些基因来说,每个样本通常有数百个读数;2 倍计数差异将非常显着)。
我希望能够使用更复杂的模型进行应急测试,也许使用 3 个不同的自变量(并作为块效应进行复制)。
简单的解决方案是这样的:我创建一个额外的分类变量,我们将其称为“模型”,其值为“观察到”和“预期”,“观察到”对应于上面的行,“预期”对应于相同数量的行,但具有以下列:如果没有显着影响,则包含预期计数的读取,即,在这种情况下,平均读取计数(假设每个重复中的工作量相等)。
然后运行这个:
拟合模型( Freq( :reads ), Y( :Model ), 效果( :Env, :rep ), 个性( 名义Logistic ), 运行( 似然比检验( 1 ), Wald 检验( 0 ) ) )
...或者这个(这给了我相同的结果,我不太明白为什么,但我对此感到满意):
拟合模型( Freq( :reads ), Y( :Model ), 效果( :Env, :rep ), 个性( 广义线性模型 ), GLM 分布( 二项式 ), 链接函数( Logit ), 过度离散检验和区间( 0 ) , 名称(“Firth 偏差调整估计值”)(0), 运行 )
我是否做了一些可笑的错误,或者这作为一个有效的测试是否有意义?
如果有意义,是否有一种方法可以实现相同的目的,而无需创建具有预期计数的额外行?
提前致谢!!
扬波尔斯基
数据探索和可视化质量和工艺工程回复:RNAseq 数据的广义线性模型或名义 Logistic谢谢,这很有帮助,现在更有意义了。我会看一下(我希望还有其他人)对该方法进行健全性检查。
回复:RNAseq 数据的广义线性模型或名义 Logistic@l_yampolsky谢谢你的提问。当将响应从二项式转换为二项式时,广义线性混合模型中也会出现同样的问题。在这种情况下,必须声明原始样本是随机效应,否则结果会完全误导。
不幸的是,您的第二次分析也是如此,因为您实际上假设每个观察都是独立的,而实际上它们是根据原始样本进行聚类的。要纠正第二个分析,请将代表效应声明为随机,如果您有 JMP Pro 并使响应连续(我建议对其进行 log2 转换),您还需要将个性更改为标准最小二乘法或混合模型。或者,如果您想保留泊松响应,Meichen Dong 有一个 GLMM 插件,可以从社区这里下载,她甚至在文档中提供了一个 RNASeq 示例:广义线性混合模型插件
底线是,你实际上只有 16 个犯错的自由度,而且你无法欺骗它。
您对 6K 转录本的多次测试的担忧是正确的。如果某些转录本是已知的并预先指定,您可以将单个转录本分析显示为第一个或修改后的第二个。对于全套,制作一个火山图,看看你有什么。另一个想法是构建一个或多个基因或通路评分(对相关转录本进行平均)并对其进行分析。
回复:RNAseq 数据的广义线性模型或名义 Logistic我在上面的线程中发布了一个示例。
预先感谢您的宝贵意见!
扬波尔斯基
回复:RNAseq 数据的广义线性模型或名义 Logistic这是一个例子:
具有标称响应示例的广义线性模型
数据文件:Example1Fp52.43_cH8.jmp
研究了人类线粒体载体蛋白 Aralar1 的甲壳类同源物在 4 个克隆中的表达水蚤暴露于急性缺氧(3次缺氧和3次对照重复)。 4 个克隆 x 2 次处理 x 3 个重复 = 24 行。
首先,使用 RPKM 作为响应变量,缺氧和克隆作为因子进行 2 路方差分析:
适合模型(
Y(:RPKM),
效果(:克隆,:急性缺氧,:克隆*:急性缺氧),
个性(标准最小二乘法),
强调(最少报告),
跑步(
:RPKM << {缺乏拟合( 0 ),根据预测绘制实际值( 0 ),
绘制回归 ( 0 ),绘制预测残差 ( 0 ),
情节效果杠杆( 0 )}
)
)
结果是:
响应RPKM
效果测试
源 Nparm DF 平方和 F 比率 Prob > F
克隆 3 3 10222.900 3.6649 0.0349
急性缺氧 1 1 42430.209 45.6339 <.0001
克隆*急性缺氧 3 3 12083.132 4.3318 0.0205
好吧,事情很重要,但克隆和相互作用的影响将永远无法在任何类型的多重测试校正中幸存(这只是 RNAseq 数据中约 6K 转录本的一个例子)。
然后,使用包含预期读取数的读取列创建 24 个虚拟行,以防没有克隆或缺氧的影响,并且在每种克隆/缺氧组合中的 3 个重复之间也没有差异 - 即,实际24行。 虚拟行在 :ObsExp 列中标记为“Exp”,原始实际观测值标记为“Obs”。 现在我们想知道 Obs 是否与 Exp 不同,如果是,则由于哪些自变量而不同。 为了实现这一点,我运行(观察比现在原始:读取,而不是 RPKM 用作频率:
适合模型(
频率(:读取),
Y(:观测值),
效果(
:克隆,
:急性缺氧,
:急性缺氧*:克隆,
:代表[:克隆,:急性缺氧]
),
个性(名义逻辑),
Run( 似然比检验 ( 1 ), Wald 检验 ( 0 ) ),
SendToReport( Dispatch( {}, "参数估计", OutlineBox, {Close( 1 )} ) )
)
正如预期的那样,其结果是更强大:
效应似然比检验
来源 Nparm DF LR 卡方 Prob>ChiSq
克隆 3 3 208.563245 <.0001
急性缺氧 1 1 114.216936 <.0001
急性缺氧*克隆 3 3 54.2042677 <.0001
代表[克隆,急性缺氧] 16 16 92.4764669 <.0001
注意:使用以下方法可以获得完全相同的结果:
人格(广义线性模型),
GLM 分布(二项式),
链接函数(Logit)
回复:RNAseq 数据的广义线性模型或名义 Logistic感谢您的回答。
我认为你无法回避进行类似 t 检验的分析。 重复是变异性估计和执行分析的自由度的来源,不需要嵌套。 JMP 术语中的嵌套如下:
如果一种效应 (B) 的水平仅出现在另一种效应 (A) 的单个水平内,则称 B 为嵌套的在 A (B[A]) 内。
这意味着您必须有一个因素 (B),使得 B 级别的某个级别只能应用于 Env 的某个级别内。 所以。因子 B 的级别 1 和 2 仅出现在 Env A 中,因子 B 的级别 3 和 4 仅出现在 Env 中。 B. 如果 Env 是模型效应,则意味着重复位于 Env 的关联水平内。 因此无需将它们显式分配到每个 Env 级别。
JMP 会让你这样做,但我不认为它在统计上有效或有意义。
此外,也没有必要将重复放入随机样本,因为这也是暗示的(三个重复应该是随机样本,并且在广义线性模型或广义回归中是不允许的)。
在这种情况下,Reads 是 Y,我认为 Freq 选项不会执行您期望的操作(通常避免该选项)。 泊松分布将是使用广义线性模型或负二项式或泊松进行广义回归 (JMP Pro) 的分布。
我想你需要为数千个基因做这个?
感谢您发送示例表。 我创建了一个表,其中包含保存的分析脚本,以便我们继续讨论。 我创建了一个附加列(复制和复制 ID)供本次讨论使用。
回复:RNAseq 数据的广义线性模型或名义 Logistic复制品来自不同的个体;生物复制。 你是绝对正确的,将它们集中起来并不是一个好主意。 然而,我可以将复制作为随机自变量纳入我的模型中,嵌套在环境处理中,本质上将此设计视为块设计 - 相同的实验重复 3 次。
我更担心一般方法 - 在零假设下创建具有预期“计数”的虚拟观察(即,各处的读取计数相等),并使用列联表/名义物流模型来查看实际计数是否不同。
我将按照另一位用户的建议,发布一个包含实际数据的示例。
回复:RNAseq 数据的广义线性模型或名义 Logistic你好@l_yampolsky ,
需要明确的是,对于每个环境,您都有三个样本(重复)。 但是,您不想组合每个环境的重复项并进行类似 t 检验的分析,而是希望将每个样本视为计数代表该样本转录本群体中特定转录本(转录本 A)的观察数量。 来自同一环境的下一个样本就像从同一转录本群体中提取一组新的观察结果,并观察有多少属于相同的转录本 A。 我理解正确吗?
问题:每个样本是来自同一个体(即该个体内细胞群体的样本,因此技术重复)还是来自同一环境中的不同个体(即每个个体的样本,因此代表样本的不同群体)经历过相同环境的细胞)?
如果样本来自不同的个体,我认为这不是有效的方法。 计数是对环境某些变化的响应。 这就像在测试环境中对来自光源的特定波长的光子数进行计数,然后将其与来自另一个具有相似光源且也在相同环境中的不同设备的相同波长的光子数进行汇总。 它们不是同一光源(或设备),而是同一光源的代表,必须相互独立对待。 这将代表设备之间的差异,并且必须单独跟踪,因为每个设备都不是另一个设备的精确副本。 如果是的话,那么测试就不需要多个设备,因为所有光子都来自单个源,并且每个样本都将是该源的代表性样本。
RNAseq 的相似之处在于它是对样本中的实体或转录本的丰度进行定量的努力。 下一个样本通常不是来自同一源总体的样本。 它通常是来自不同个体的样本,因此必须以这种方式处理。 个体与个体之间存在差异。 即使来自同一个体的相同组织的样本也具有变异性(细胞处于不同的分裂状态、功能分层等),并且不应将它们汇集在一起,就好像它们来自同一细胞群一样。
最终,这一切实际上取决于所提出的科学问题以及用于发现该问题答案的实验设计。 您也许可以将它们集中起来,做您想做的事,也可能不做。 这取决于....
回复:RNAseq 数据的广义线性模型或名义 Logistic如果您附上包含不同分析的示例数据表,可能会有所帮助。 我发现这比尝试从文本描述中想象数据表更容易理解。