时间戳

首页 > 解决方案 > edgeR：设计矩阵未满秩

r - edgeR：设计矩阵未满秩

问题描述

我正在分析具有四个时间点和两个处理的 RNA seq 数据。在多个时间点测量样品并始终得到相同的处理。我想做一个配对模型。以下代码总是给我一个错误说明，即设计矩阵不是满秩的。为什么这个设计矩阵不是满秩的？

|--------|-----------|------------|
| sample | treatment | time_point |
|--------|-----------|------------|
|   1    |    a      |      1     |
|   1    |    a      |      2     |
|   1    |    a      |      3     |
|   1    |    a      |      4     |
|   2    |    b      |      1     |
|   2    |    b      |      2     |
|   2    |    b      |      3     |
|   2    |    b      |      4     |
|  ...   |   ...     |     ...    |
|--------|-----------|------------|

merged_vars <- factor(paste(treatment, time_point, sep="."))
design <- stats::model.matrix(~sample+merged_vars)
y <- edgeR::estimateDisp(y, design = design)

标签： r

解决方案

当您的一个预测变量可以表示为另一个的线性组合时，就会出现错误。

在您的示例中，如果您查看 model.matrix，则 b.1 + b.2 + b.3 + b.4 会给您示例 2 的效果：

df = data.frame(sample=rep(1:2,each=4),
treatment=rep(letters[1:2],each=4),time_point=rep(1:4,2))
df$merged_vars= with(df,factor(paste(treatment, time_point, sep=".")))

model.matrix(~sample +merged_vars,data=df)

  (Intercept) sample merged_varsa.2 merged_varsa.3 merged_varsa.4
1           1      1              0              0              0
2           1      1              1              0              0
3           1      1              0              1              0
4           1      1              0              0              1
5           1      2              0              0              0
6           1      2              0              0              0
7           1      2              0              0              0
8           1      2              0              0              0
  merged_varsb.1 merged_varsb.2 merged_varsb.3 merged_varsb.4
1              0              0              0              0
2              0              0              0              0
3              0              0              0              0
4              0              0              0              0
5              1              0              0              0
6              0              1              0              0
7              0              0              1              0
8              0              0              0              1

此外，如果您这样做merged_vars <- factor(paste(treatment, time_point, sep="."))，则每个级别有 1 个观察值，这意味着没有重复，您不能对它进行 egdeR。

我建议这样做，因为它看起来像样本和治疗混淆：

design <- stats::model.matrix(~treatment*time_point)

推荐阅读