科研代码分享|R的Limma包实现表达谱数据标准化和差异基因提取
科研有捷径，输入代码，一键获取科研成果！就是这么省事，来具体看下有多方便！
搜索http://985.so/a9kb查看全部代码（目前共计50+持续新增中），也可以点击右侧【目录】，可以看到更多有趣的代码
真香提示：文末可以知道如何获取代码~
在处理高通量数据过程中，我们往往需要借助一些表达谱信息进行比对分析，通过不同类型样本之间的比较，获得与疾病表型存在显著相关性的基因标识。
GEO数据库http://www.ncbi.nlm.nih.gov/geo/为我们提供了大量和疾病相关的表达谱信息，其包含物种丰富（人，鼠， 植物等），表达谱类型多样（芯片数据，甲基化数据， RNA-seq数据等），并且为使用者提供了各种可下载的数据资源。其中比较常用的一种下载资源是CEL格式文件，即原始数据，cel文件就是Affymetrix扫描仪产出的原始数据文件，记录了每个probe 的signal intensity。
这里我们通过从GEO数据库下载的一个样例，逐步介绍RMA进行CEL表达谱数据标准化，以及limma算法提取差异表达基因的过程。
首先我们从GEO下载瘢痕瘤的表达谱数据GSE7890，下载地址http://www.ncbi.nlm.nih.gov/geo/query/acc.cgiacc=GSE7890，下载CEL文件到本地，解压，这里默认目录为F:\。之后调用R语言的RMA包，标准化过程如下：
直方图为标准化前后数据分布的比较，可以观察出基因的以2为底对数转换后的信号值集中在什么区间
箱式图为RMA方法消除样本间差异前后的数据比较，红色图代表处理前数据分布，蓝色为处理后数据分布，可见RMA方法可以有效去除样本间差异。
实现代码如下：
source("http://bioconductor.org/biocLite.R")
biocLite("affy")
biocLite("Biobase")
biocLite("tkWidgets")
library(affy)
AffyData<-ReadAffy(widget=TRUE)
exprsSet.RMA<-rma(AffyData) ##用rma方法处理数据
exp.RMA<-exprs(exprsSet.RMA) ## exprsSet.MAS5<-mas5(AffyData) ##用mas5方法处理数据
write.exprs(exprsSet.RMA, file="mydata.txt") ##结果保存
hist(log2(intensity(AffyData[,4])),breaks=100,col="blue") ##画芯片数据的直方图，因为实验数据变化范围大所以取log2，结果集中在6~8，大于10 的很少，说明样本仅少部分高表达，与实际吻合。
par(mfrow=c(1,2))
boxplot(AffyData,col="red")
boxplot(data.frame(exp.RMA),col="blue") ## 画箱式图，比较数据分布情况
plot(exp.RMA[,1:2],log="xy",pch=".",main="all")
plot(pm(AffyData)[,1:2],log="xy",pch=".",main="pm")
plot(mm(AffyData)[,1:2],log="xy",pch=".",main="mm")
plot(exp.RMA[,1:2],pch=".",main="exp.RMA") ## 散点图评价两个样本之间的数据差异
library(genefilter) ## 函数rowSds在genenfliter中
rsd<-rowSds(exp.RMA) ## 计算行的标准差
sel<-order(rsd,decreasing=TRUE)[1:20] ## 选择表达差异高的top20，order返回行号
heatmap(exp.RMA[sel,],col=rainbow(256)) ## 画热图
之后我们获得了标准化后的表达谱数据，接下来我们可以利用limma package的算法对表达谱数据进行差异基因提取。
需要设置的参数包括表达谱数据的文件目录，case组和control组的数据，以及输出文件的文件名，这里默认为case_control_limma.txt。那么对于我们的样例文件，我们将其放置于F：\盘下，case组（treated with**）为9个，control组（treated without**）为10个，实现的代码如下：
library(limma)
setwd("F:\\")
exp<-read.table("profile_normalized.txt",header=T,sep="\t")
rownames(exp)<-exp[,1]
exp<-exp[2:20]
exp<-as.matrix(exp)
samps<-factor(c(rep("case",9),rep("control",10)))
design <- model.matrix(~0+samps) ;
colnames(design) <- c("case","control")
fit <- lmFit(exp, design)
cont.matrix<-makeContrasts(case-control,levels=design)
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)
final<-topTable(fit2, coef=1, number=dim(exp)[1], adjust.method="BH", sort.by="B", resort.by="M")
write.table(final,paste("case_control","_limma.txt"),quote=FALSE,sep="\t")
最后输出的结果如下图所示：
表中共包含8列，分别代表探针ID，校正后P值，校正前P值，t/B统计量，对数转换的fold change值，gene symbol以及基因名称。接下来使用者可根据需要，设定P值和logFC的阈值，从而提取出显著差异表达的基因。
当然数据标准化和提取差异基因的方法很多，例如Z-score标准化，通过计算均值和标准差进行零均值校正，提出基因固有表达差异的影响也是一种常用的标准化方法，另外还有类似Python语言中preprocessing包中的scale，normalize等方法。
差异基因的提取方法还包括student T test，秩和检验，SAM package，甚至更高级一些的机器学习中的feature selection也可以实现这一目的。使用者可根据需要或数据类型进行选择。
