Hello，大家好！此前，小编给大家介绍了：如何从GEO数据库下载芯片数据及相关样本处理信息等等（不知道的可以戳这里哦）。这些芯片数据下载下来干嘛呢？下载必然是为了挖数据啦！指不定什么有意思的东西我就发现了呢？指不定老天爷眼睛一闭让我发了文章了呢？
言归正传，拿到数据，分析的第一步往往是进行基因差异表达分析，所以，针对芯片数据，我们就来给大家介绍一款基因差异表达分析的常用方法——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)
GSM49939GSM49940GSM49941GSM49942GSM49943GSM49944
1007_s_at23.013.826.575.994.984.6
1053_at1449.91826.72242.81508.81523.02355.5
117_at109.271.5106.7128.884.179.6
针对表达矩阵，需要查看其整体分布情况，可以利用boxplot()绘制box分布图，GEO下载的表达矩阵数据基本上都是标准化的数据，可以由箱线图的分布特点看出这些样本的数据基本分布一致（中位数、上四分位数、下四分位数等等），如下图结果：
#获取样品数量，并设置图片颜色
n.sample = ncol(exprSet)
cols = rainbow(n.sample)
#利用boxplot()绘图
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
GSM4993910
GSM4994010
GSM4994110
GSM4994201
GSM4994301
GSM4994401
attr(,"assign")
[1]11
attr(,"contrasts")
attr(,"contrasts")$Treatment
[1]"contr.treatment"
3、差异比较矩阵
基于分组矩阵的信息构建差异比较矩阵（cont.matrix)，由差异比较矩阵显示结果可知，是进行INFECTED 与CONTROL之间的差异分析。
>cont.matrix = makeContrasts(INFECTED-CONTROL, levels=design)
>cont.matrix
Contrasts
LevelsINFECTED-CONTROL
CONTROL-1
INFECTED1
差异表达分析
差异表达分析主要是基于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.ValP.ValuelogFC
223020_at0.999642.196175e-05746.10000
1555758_a_at0.999646.467722e-05-540.53333
218676_s_at0.999641.352768e-04-280.86667
237249_at0.999642.669173e-04-93.53333
225100_at0.999642.836527e-04-124.96667
217825_s_at0.999642.903446e-04-143.73333
222099_s_at0.999643.425427e-04493.13333
212634_at0.999644.221452e-04-166.06667
211499_s_at0.999644.391776e-04-129.56667
221098_x_at0.999644.805746e-0495.16667
208974_x_at0.999645.060448e-04947.76667
209670_at0.999645.113338e-04374.20000
202088_at0.999645.262646e-04-594.40000
219394_at0.999645.307063e-04-117.56667
212221_x_at0.999645.393084e-04347.43333
以上，就完成了分析过程，之后可以根据结果设定合适的阈值筛选出差异表达基因，并基于结果进行图形化展示，或者进行富集分析、蛋白质互作网络分析等等，如此可以丰富分析结果，感兴趣的同学可以去探索一下哦！
更多生物信息课程：
1. 文章越来越难发？是你没发现新思路，基因家族分析发2-4分文章简单快速，学习链接：基因家族分析实操课程、基因家族文献思路解读
2. 转录组数据理解不深入？图表看不懂？点击链接学习深入解读数据结果文件，学习链接：转录组（有参）结果解读；转录组（无参）结果解读
3. 转录组数据深入挖掘技能-WGCNA，提升你的文章档次，学习链接：WGCNA-加权基因共表达网络分析
4. 转录组数据怎么挖掘？学习链接：转录组标准分析后的数据挖掘、转录组文献解读
5. 微生物16S/ITS/18S分析原理及结果解读、OTU网络图绘制、cytoscape与网络图绘制课程
6. 生物信息入门到精通必修基础课，学习链接：linux系统使用、perl入门到精通、perl语言高级、R语言画图
7. 医学相关数据挖掘课程，不用做实验也能发文章，学习链接：TCGA-差异基因分析、GEO芯片数据挖掘、GSEA富集分析课程、TCGA临床数据生存分析、TCGA-转录因子分析、TCGA-ceRNA调控网络分析
8.其他课程链接：二代测序转录组数据自主分析、NCBI数据上传、二代测序数据解读。
网友评论