用于基因表达的Limma包

Question

我有以下格式的rma归一化矩阵：

ID_REF GSM362180    GSM362181  GSM362188    GSM362189  GSM362192
244901 5.094871713 4.626623079 4.554272515 4.748604391 4.759221647
244902 5.194528083 4.985930299 4.817426064 5.151654407 4.838741605
244903 5.412329253 5.352970877 5.06250609  5.305709079 8.365082403
244904 5.529220594 5.28134657  5.467445095 5.62968933  5.458388909
244905 5.024052699 4.714631878 4.792865831 4.843975286 4.657188246
244906 5.786557533 5.242403911 5.060605782 5.458148567 5.890061836

其中不同的列对应于四种不同类型的启动子，并且四种启动子中的每一种都具有生物复制，所以总共有8列。

我尝试使用Limma软件包来寻找几个启动子(带有重复)的差异表达基因，但我总是得到一个错误，因为我对r来说是新的，无法完全理解它。

这是我使用的代码：

Group <- factor(c("p1", "p1", "p2", "p2", "p3","p3","p3","p4","p4"), levels = c("GSM362180","GSM362181","GSM362188","GSM362189","GSM362192","GSM362193","GSM362194","GSM362197","GSM362198"))

design <- model.matrix(~0 + Group)

colnames(design) <- c("GSM362180","GSM362181","GSM362188","GSM362189","GSM362192","GSM362193","GSM362194","GSM362197","GSM362198")
fit <- lmFit(modified, design)

其中modified是以上述格式输入的rma归一化数据矩阵。我得到以下错误：

Coefficients not estimable: GSM362180 GSM362181 GSM362188 GSM362189 GSM362192 GSM362193 GSM362194 GSM362197 GSM362198 

lmfit(设计，t(M))错误:0(非NA)情况

Answer 1

应该是：

Group <- factor(c("p1", "p1", "p2", "p2", "p3", "p3", "p3", "p4", "p4"))

分组类别顺序应与样本在colnames(modified)中的顺序对应

在您的原始代码中，您为Group中的每个因子分配了不同的唯一levels，因此基本上您没有分组类别，没有重复，也无法估计系数。

创建您的设计矩阵：

design <- model.matrix(~0 + Group)

然后，您希望将samples的名称分配给设计矩阵的rownames，而不是colnames，因为列代表分组。

rownames(design) <- c("GSM362180","GSM362181","GSM362188","GSM362189","GSM362192","GSM362193","GSM362194","GSM362197","GSM362198")
# or simply
rownames(design) <- colnames(modified)

您的设计矩阵应该如下所示，其中包含一个1，其中一个样本属于特定的组类别。

design
          Groupp1 Groupp2 Groupp3 Groupp4
GSM362180       1       0       0       0
GSM362181       1       0       0       0
GSM362188       0       1       0       0
GSM362189       0       1       0       0
GSM362192       0       0       1       0
GSM362193       0       0       1       0
GSM362194       0       0       1       0
GSM362197       0       0       0       1
GSM362198       0       0       0       1
attr(,"assign")
[1] 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$Group
[1] "contr.treatment"