例:查看pData
15280238-a5686fcfdc47a321.png
其中title部分告訴了我們分組信息,2小時和18小時缔刹,每個時間段又有vehicle control, PE1.3 embolized, PE2.0 embolized,也就是2x2的雙因素試驗設計, 我們可以現(xiàn)在R語言里構建實驗設計的數(shù)據(jù)框沮明。
sample <- pData$geo_accession
treat_time <- rep(c("2h","18h"),each=11)
treat_type <- rep(rep(c("vehicle_control","PE1.3_embolized","PE2.0_embolized"), c(3,4,4)), times=2)
design_df <- data.frame(sample, treat_time, treat_type)
根據(jù)Limma的使用手冊的"9.5 Interaction Models: 2 X 2 Factorial Design"進行手續(xù)的分析本今。這里僅僅展示單個因素的分析過程闻蛀,多個因素看文檔依葫蘆畫瓢就行宰僧。
構建單因素試驗設計矩陣材彪,進行線性擬合:
TS <- paste(design_df$treat_time, design_df$treat_type, sep=".")
TS <- factor(TS, levels = unique(TS))
后續(xù)可有兩種方法:
design <- model.matrix(~0+TS)
fit <- lmFit(exprSet, design)
然后設置比較對象。比如不同時間段控制組哪些基因發(fā)生了差異表達琴儿、處理18小時后處理組相對于對照組有哪些基因發(fā)生差異表達段化,也就是做多次t檢驗。
cont.matrix <- makeContrasts(
vs1 = TS18.vehicle_control-TS2.vehicle_control,
vs2 = TS18.PE2.0_embolized-TS2.PE2.0_embolized,
vs3 = TS18.PE1.3_embolized-TS2.PE1.3_embolized,
diff = (TS18.PE2.0_embolized-TS18.vehicle_control)-(TS18.PE1.3_embolized-TS18.vehicle_control), levels = design)
# {{vs1:對照組在前后的差異表達基因;
#vs2: PE2.0處理前后的差異基因;
# vs3: PE1.3在處理前后差異基因;
# 處理18小時后造成,PE2.0相對于對照變化的基因再與PE1.3與對照的差異比較}}
fit2 <- contrasts.fit(fit, cont.matrix)
results <- decideTests(fit2)
design=model.matrix(~factor(TS))
fit=lmFit(exprSet,design)
fit=eBayes(fit)
options(digits = 4)
topTable(fit,coef=2,adjust='BH')
沒有實際操作之前显熏,覺得簡單的更清爽、適用晒屎,但使用后是第一種更可取喘蟆。差異基因fold change可能搞反!9穆场履肃!在沒有顯示指定的情況下,我們難以真正確定我們比對的結果是High-Low還是Low-High坐桩。另外,前一種方法更利于將差異的比較過程程序化封锉。