2012-10-17 23 views
0

我有以下列格式的RMA歸一化的矩陣:LIMMA包用於基因表達

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軟件包在幾個啓動子中找到差異表達的基因(帶有重複數據),並且我總是得到一個錯誤,因爲Iam是新的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) 

其中改性是如上述格式輸入的RMA標準化數據矩陣。 我收到以下錯誤:

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

錯誤lmfit(設計,T(M)):0(非NA)的情況下

+1

您應該在factor因子函數中用'labels'替換'levels'。否則,'Group'將只包含'NA'。 –

+0

Idid嘗試替換它們,但我得到了幾個錯誤,調整後的p值非常高。我使用RMA標準化值,並想知道在Limma分析之前是否需要記錄轉換它們。 –

+1

除非你提供你的對象'修改'或至少'str(修改)'的結果,否則很難回答你的問題。 –

回答

0

應該僅僅是:

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

組類別訂單應對應於樣本的順序colnames(modified)

在您的原始代碼中,您爲中的每個因子分配了不同的唯一所以基本上你沒有分組類別,沒有重複,也不能估計係數。

創建您的設計矩陣:

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

然後你要爲列代表分組分配樣本的名稱,設計矩陣的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"