怎么用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 GSM499441007_s_at     23.0     13.8     26.5     75.9     94.9     84.61053_at     1449.9   1826.7   2242.8   1508.8   1523.0   2355.5117_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         TreatmentGSM49939   CONTROLGSM49940   CONTROLGSM49941   CONTROLGSM49942  INFECTEDGSM49943  INFECTEDGSM49944  INFECTED

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

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

3、差异比较矩阵

基于分组矩阵的信息构建差异比较矩阵（cont.matrix)，由差异比较矩阵显示结果可知，是进行INFECTED 与CONTROL之间的差异分析。

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

> cont.matrix          ContrastsLevels     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      logFC223020_at      0.99964 2.196175e-05  746.100001555758_a_at   0.99964 6.467722e-05 -540.53333218676_s_at    0.99964 1.352768e-04 -280.86667237249_at      0.99964 2.669173e-04  -93.53333225100_at      0.99964 2.836527e-04 -124.96667217825_s_at    0.99964 2.903446e-04 -143.73333222099_s_at    0.99964 3.425427e-04  493.13333212634_at      0.99964 4.221452e-04 -166.06667211499_s_at    0.99964 4.391776e-04 -129.56667221098_x_at    0.99964 4.805746e-04   95.16667208974_x_at    0.99964 5.060448e-04  947.76667209670_at      0.99964 5.113338e-04  374.20000202088_at      0.99964 5.262646e-04 -594.40000219394_at      0.99964 5.307063e-04 -117.56667212221_x_at    0.99964 5.393084e-04  347.43333

关于"怎么用R语言的limma方法进行芯片数据差异表达分析"的内容就介绍到这里了，感谢大家的阅读。如果想了解更多行业相关的知识，可以关注行业资讯频道，小编每天都会为大家更新不同的知识点。