怎么用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方法进行芯片数据差异表达分析”的内容就介绍到这里了，感谢大家的阅读。如果想了解更多行业相关的知识，可以关注亿速云行业资讯频道，小编每天都会为大家更新不同的知识点。