怎么用R語言的limma方法進行芯片數據差異表達分析

這篇文章主要介紹“怎么用R語言的limma方法進行芯片數據差異表達分析”的相關知識，小編通過實際案例向大家展示操作過程，操作方法簡單快捷，實用性強，希望這篇“怎么用R語言的limma方法進行芯片數據差異表達分析”文章能幫助大家解決問題。

數據簡介與設置

為了方便演示，這里選擇了人的早幼粒細胞白血病細胞系NB4細胞的六個樣本數據（GSE2600），分析的輸入文件是下載的表達矩陣文件，而分析之前需要確保正確安裝和加載limma，同時需要對工作路徑進行設置。

library('limma')
workdir="F:/GEO/20180520"
setwd(workdir)

數據處理

1、表達矩陣
數據為六個樣本，讀取數據之后，大家可以利用head()簡單查看數據的情況等。

> expreSet=read.csv2("GSE2600expressionMatrix.csv", header = T, row.names = 1,check.names = F)

> head(exprSet,3)
          GSM49939 GSM49940 GSM49941 GSM49942 GSM49943 GSM49944
1007_s_at     23.0     13.8     26.5     75.9     94.9     84.6
1053_at     1449.9   1826.7   2242.8   1508.8   1523.0   2355.5
117_at       109.2     71.5    106.7    128.8     84.1     79.6

針對表達矩陣，需要查看其整體分布情況，可以利用boxplot()繪制box分布圖，GEO下載的表達矩陣數據基本上都是標準化的數據，可以由箱線圖的分布特點看出這些樣本的數據基本分布一致（中位數、上四分位數、下四分位數等等），如下圖結果：

n.sample = ncol(exprSet)
cols = rainbow(n.sample)
pdf(file=paste(workdir,"/","Probe_expressionDistribution.pdf",sep=""), width=24, height=18)
par(cex = 0.7)
if(n.sample>40) par(cex = 0.5)
boxplot(exprSet,col = cols, main = "expression", las = 2)
dev.off()

2、分組矩陣

確認表達矩陣之后，可以由下載保存的樣本處理信息進行分組，例如此處的樣本處理分組：CONTROL/INFECTED，經過整理，分組信息大致如下,并基于分組信息構建分組矩陣（design):

> group
         Treatment
GSM49939   CONTROL
GSM49940   CONTROL
GSM49941   CONTROL
GSM49942  INFECTED
GSM49943  INFECTED
GSM49944  INFECTED

> design = model.matrix(~ Treatment + 0, group)
> colnames(design) = levels(as.factor(c("CONTROL","INFECTED")))

> design
         CONTROL INFECTED
GSM49939       1        0
GSM49940       1        0
GSM49941       1        0
GSM49942       0        1
GSM49943       0        1
GSM49944       0        1
attr(,"assign")
[1] 1 1
attr(,"contrasts")
attr(,"contrasts")$Treatment
[1] "contr.treatment"

3、差異比較矩陣

基于分組矩陣的信息構建差異比較矩陣（cont.matrix)，由差異比較矩陣顯示結果可知，是進行INFECTED 與CONTROL之間的差異分析。

>cont.matrix = makeContrasts(INFECTED-CONTROL, levels=design)

> cont.matrix
          Contrasts
Levels     INFECTED - CONTROL
  CONTROL                  -1
  INFECTED                  1

差異表達分析

差異表達分析主要是基于lmFit()、eBayes()、topTable()完成分析過程，并提取了主要的結果(tT)。

> fit = lmFit(exprSet, design)
> fit2 = contrasts.fit(fit, cont.matrix)
> fit2 = eBayes(fit2, 0.01)
> tT = topTable(fit2, adjust="fdr", sort.by="logFC", resort.by = "P" ,n=Inf)

> tT = subset(tT, select=c("adj.P.Val","P.Value","logFC"))

> head(tT,15)
             adj.P.Val      P.Value      logFC
223020_at      0.99964 2.196175e-05  746.10000
1555758_a_at   0.99964 6.467722e-05 -540.53333
218676_s_at    0.99964 1.352768e-04 -280.86667
237249_at      0.99964 2.669173e-04  -93.53333
225100_at      0.99964 2.836527e-04 -124.96667
217825_s_at    0.99964 2.903446e-04 -143.73333
222099_s_at    0.99964 3.425427e-04  493.13333
212634_at      0.99964 4.221452e-04 -166.06667
211499_s_at    0.99964 4.391776e-04 -129.56667
221098_x_at    0.99964 4.805746e-04   95.16667
208974_x_at    0.99964 5.060448e-04  947.76667
209670_at      0.99964 5.113338e-04  374.20000
202088_at      0.99964 5.262646e-04 -594.40000
219394_at      0.99964 5.307063e-04 -117.56667
212221_x_at    0.99964 5.393084e-04  347.43333

關于“怎么用R語言的limma方法進行芯片數據差異表達分析”的內容就介紹到這里了，感謝大家的閱讀。如果想了解更多行業相關的知識，可以關注億速云行業資訊頻道，小編每天都會為大家更新不同的知識點。