怎么用R語(yǔ)言的limma方法進(jìn)行芯片數(shù)據(jù)差異表達(dá)分析

這篇文章主要介紹“怎么用R語(yǔ)言的limma方法進(jìn)行芯片數(shù)據(jù)差異表達(dá)分析”的相關(guān)知識(shí)，小編通過(guò)實(shí)際案例向大家展示操作過(guò)程，操作方法簡(jiǎn)單快捷，實(shí)用性強(qiáng)，希望這篇“怎么用R語(yǔ)言的limma方法進(jìn)行芯片數(shù)據(jù)差異表達(dá)分析”文章能幫助大家解決問(wèn)題。

數(shù)據(jù)簡(jiǎn)介與設(shè)置

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

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

數(shù)據(jù)處理

1、表達(dá)矩陣
數(shù)據(jù)為六個(gè)樣本，讀取數(shù)據(jù)之后，大家可以利用head()簡(jiǎn)單查看數(shù)據(jù)的情況等。

> 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

針對(duì)表達(dá)矩陣，需要查看其整體分布情況，可以利用boxplot()繪制box分布圖，GEO下載的表達(dá)矩陣數(shù)據(jù)基本上都是標(biāo)準(zhǔn)化的數(shù)據(jù)，可以由箱線(xiàn)圖的分布特點(diǎn)看出這些樣本的數(shù)據(jù)基本分布一致（中位數(shù)、上四分位數(shù)、下四分位數(shù)等等），如下圖結(jié)果：

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、分組矩陣

確認(rèn)表達(dá)矩陣之后，可以由下載保存的樣本處理信息進(jìn)行分組，例如此處的樣本處理分組：CONTROL/INFECTED，經(jīng)過(guò)整理，分組信息大致如下,并基于分組信息構(gòu)建分組矩陣（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、差異比較矩陣

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

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

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

差異表達(dá)分析

差異表達(dá)分析主要是基于lmFit()、eBayes()、topTable()完成分析過(guò)程，并提取了主要的結(jié)果(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

關(guān)于“怎么用R語(yǔ)言的limma方法進(jìn)行芯片數(shù)據(jù)差異表達(dá)分析”的內(nèi)容就介紹到這里了，感謝大家的閱讀。如果想了解更多行業(yè)相關(guān)的知識(shí)，可以關(guān)注億速云行業(yè)資訊頻道，小編每天都會(huì)為大家更新不同的知識(shí)點(diǎn)。