“差异”是个统计学概念，获取差异表达基因就要用统计方法，R的统计功能很强大，适合做这样的事情。用前面的方法读取数据：
library(affy)library(tcltk)filters<-matrix(c("CEL file",".[Cc][Ee][Ll]","All",".*"),ncol=2,byrow= T)cel.files<-tk_choose.files(caption="Select CELs",multi=TRUE,filters= filters,index=1)data.raw<-ReadAffy(filenames= cel.files)n.cel<-length(cel.files)# 查看样品名称sampleNames(data.raw)
## [1] "NRID9780_Zarka_2-1_MT-0HCA(SOIL)_Rep1_ATH1.CEL"
## [2] "NRID9781_Zarka_2-2_MT-0HCB(SOIL)_Rep2_ATH1.CEL"
## [3] "NRID9782_Zarka_2-3_MT-1HCA(SOIL)_Rep1_ATH1.CEL"
## [4] "NRID9783_Zarka_2-4_MT-1HCB(SOIL)_Rep2_ATH1.CEL"
## [5] "NRID9784_Zarka_2-5_MT-24HCA(SOIL)_Rep1_ATH1.CEL"
## [6] "NRID9785_Zarka_2-6_MT-24HCB(SOIL)_Rep2_ATH1.CEL"
## [7] "NRID9786_Zarka_2-7_MT-7DCA(SOIL)_Rep1_ATH1.CEL"
## [8] "NRID9787_Zarka_2-8_MT-7DCB(SOIL)_Rep2_ATH1.CEL"
# 简化一下名称，设置pDatasampleNames(data.raw)<-paste("S",1:n.cel,sep='')pData(data.raw)$treatment<-rep(c("0h","1h","24h","7d"),each=2)pData(data.raw)
## sample treatment
## S1 1 0h
## S2 2 0h
## S3 3 1h
## S4 4 1h
## S5 5 24h
## S6 6 24h
## S7 7 7d
## S8 8 7d
使用rma和mas5方法进行预处理：
eset.rma<-rma(data.raw)
## Background correcting
## Normalizing
## Calculating Expression
eset.mas5<-mas5(data.raw)
## background correction: mas
## PM/MM correction : mas
## expression values: mas
## background correcting...done.
## 22810 ids to be processed
## | |
## |####################|
很简单，用一个exprs函数就可以从eset数据中提取出表达量，得到的数据类型是矩阵。但是应该注意rma的eset结果是经过对数变换的，而mas5的eset结果是原始信号强度。虽然表达量是用对数变换的信号值表示的，但是有些计算过程要用到未经变换的原始值，应该把它们都计算出来：
emat.rma.log2<-exprs(eset.rma)emat.mas5.nologs<-exprs(eset.mas5)class(emat.rma.log2)
## [1] "matrix"
head(emat.rma.log2,1)
## S1 S2 S3 S4 S5 S6 S7 S8
## 244901_at 4.04 4.348 4.048 4.052 4.019 3.962 4.03 4.062
head(emat.mas5.nologs,1)
## S1 S2 S3 S4 S5 S6 S7 S8
## 244901_at 39.79 94.94 57.35 60.23 61.08 55.57 53.95 85.98
emat.rma.nologs<-2^emat.rma.log2emat.mas5.log2<-log2(emat.mas5.nologs)
下面我们仅使用rma的结果做演示。计算平均表达量和差异表达倍数（和0h对照比）：
rm(eset.mas5)rm(emat.mas5.nologs)rm(emat.mas5.log2)#计算平均值，并做对数转换results.rma<-data.frame((emat.rma.log2[,c(1,3,5,7)]+emat.rma.log2[,c(2,4,6,8)])/2)#计算表达量差异倍数results.rma$fc.1h<-results.rma[,2]-results.rma[,1]results.rma$fc.24h<-results.rma[,3]-results.rma[,1]results.rma$fc.7d<-results.rma[,4]-results.rma[,1]head(results.rma,2)
## S1 S3 S5 S7 fc.1h fc.24h fc.7d
## 244901_at 4.194 4.050 3.991 4.046 -0.1448 -0.2037 -0.1481
## 244902_at 4.293 4.159 4.061 3.937 -0.1340 -0.2316 -0.3557
简单补充介绍一下R语言中取数据子集的三种方法，主要是矩阵和数据框：
用下标子集取数据子集::比如上面用到的eset.rma[, c(1,3,5,7)]。由于eset.rma是2维矩阵，eset.rma[, c(1,3,5,7)]的第一维留空（逗号前不写东西）表示取全部的行，第二维下标的取值为向量c(1,3,5,7)，表示取1,3,5,7共4列。
用行、列名称取子集::eset.rma[c("244901_at", "244902_at"), ]的第一维（行）是名称向量为c("244901_at", "244902_at")，第二维留空，表示取数据中行名称为c("244901_at", "244902_at")的所有列。同样方法可应用在列选取上。
用逻辑向量取子集::比如我们要选取results.rma中fc.7d大于0的所有行，分两步：先产生一个逻辑向量，然后用这个逻辑向量取子集，也可以一步完成。
subset.logic<-results.rma$fc.7d>0subset.data<-results.rma[subset.logic,]
要注意的是逻辑向量的长度要和相应维度的数据长度一致，逻辑向量中为TRUE的就保留，FALSE的就丢弃：
length(subset.logic);nrow(results.rma)
## [1] 22810
head(subset.logic)
## [1] FALSE FALSE FALSE FALSE TRUE TRUE
很多人可能不做这一步，认为没必要。但是理论上说，芯片可以检测的基因不一定在样品中都有表达，对于样本量较小的非“pooled”样品来说这是肯定的。把芯片上所有基因当成样本中的表达基因去做统计分析显然不合适。
选取“表达”基因的方法常见的有两种，一是使用genefilter软件包，另外一种是调用affy包的mas5calls()函数。使用 genefilter需要设定筛选阈值，不同的人可能有不同的标准，稍嫌随机，不够自动化，不介绍了。mas5calls方法使用探针水平数据（AffyBatch类型数据）进行处理，一般使用没经过预处理的芯片数据通用性强些，其他参数用默认就可以：
data.mas5calls<-mas5calls(data.raw)
## Getting probe level data...
## Computing p-values
## Making P/M/A Calls
继续用exprs计算“表达”量，得到的数据只有三个值P/M/A。对于这三个值的具体解释可以用?mas5calls查看帮助。P为present，A为absent，M为marginal（临界值）。
eset.mas5calls<-exprs(data.mas5calls)head(eset.mas5calls)
## S1 S2 S3 S4 S5 S6 S7 S8
## 244901_at "A" "P" "P" "A" "P" "P" "A" "P"
## 244902_at "A" "P" "P" "M" "A" "A" "P" "A"
## 244903_at "P" "P" "P" "P" "P" "P" "P" "P"
## 244904_at "A" "A" "A" "A" "A" "A" "A" "A"
## 244905_at "A" "A" "A" "A" "A" "A" "A" "A"
## 244906_at "A" "A" "A" "A" "A" "A" "A" "A"
下面我们把至少一个芯片中有表达的基因选出来：从22810中选出了16005个。
AP<-apply(eset.mas5calls,1,function(x)any(x=="P"))present.probes<-names(AP[AP])paste(length(present.probes),"/",length(AP))
## [1] "16005 / 22810"
删掉一些中间数据很必要：
rm(data.mas5calls)rm(eset.mas5calls)
present.probes是名称向量，用它进行数据子集提取：
results.present<-results.rma[present.probes,]
生物学数据分析时的"差异"应该有两个意思，一是统计学上的差异，另外一个是生物学上的差异。一个基因在两个条件下的表达量分别有3个测量值：99,100,101 和 102,103,104。统计上两种条件下的基因表达数值是有差异的，后者比前者表达量要大。但生物学上有意义吗？未必。按平均值计算表达变化上升了3%，能产生什么样的生物学效应？这得看是什么基因了。所以差异表达基因的选取一般设置至少两个阈值：基因表达变化量和统计显著性量度（p值、q值等）。
这种方法不用太多的统计学知识，生物专业的人很容易想到，而且确实有不少人在用。经常使用的筛选阈值是表达量变化超过2倍，即|log2(fc)|>=log(2)。先简单看看有没有：
apply(abs(results.present[,5:7]),2, max)
## fc.1h fc.24h fc.7d
## 5.309 6.688 6.844
apply是一个很有用的函数，它对数据按某个维度批量应用一个函数进行计算。第一个参数为向量或矩阵（或者是能转成向量或矩阵的数据，如数据框），第三个参数表示要使用的函数，第二个参数为应用的维度。上面语句的意思是对数据 abs(results.present[,5:7]) 按列（第二维）使用统计函数max（计算最大值）。表达变化超过2倍的基因共有842个：
sum(abs(results.present[,"fc.7d"])>=log2(2))
## [1] 842
选出这842个基因：
results.st<-results.present[abs(results.present$fc.7d)>=log2(2),]sel.genes<-row.names(results.st)
t测验，并选出p<0.05的差异表达基因：
p.value<-apply(emat.rma.log2[sel.genes,],1,function(x){t.test(x[1:2], x[7:8])$p.value})results.st$p.value<-p.valuenames(results.st)
## [1] "S1" "S3" "S5" "S7" "fc.1h" "fc.24h" "fc.7d"
## [8] "p.value"
results.st<-results.st[,c(1,4,7,8)]results.st<-results.st[p.value<0.05,]head(results.st,2)
## S1 S7 fc.7d p.value
## 245042_at 8.153 7.021 -1.133 0.01004
## 245088_at 7.041 5.419 -1.622 0.03381
nrow(results.st)
## [1] 347
通过简单t测验方法得到347个表达倍数变化超过2倍的差异表达基因。
3.2SAM（Significance Analysis of Microarrays）
R软件包samr可以做这个。得先安装：
library(BiocInstaller)biocLite("samr")
这种方法流行过一段时间，但由于FDR（错误检出率）控制太差，现在基本不用了。
要用也不复杂。但是注意SAM函数使用的emat表达数据是present.probes筛选出来的“表达”基因子集，如果你用没有经过筛选的数据，得到的结果会差别很大，不信可以自己试试（这点可能也是这种方法的毛病之一）：
library(samr)samfit<-SAM(emat.rma.nologs[present.probes,c(1,2,7,8)],c(1,1,2,2),resp.type="Two class unpaired",genenames=present.probes)
## perm= 1
## perm= 2
## perm= 3
## perm= 4
## perm= 5
## perm= 6
## perm= 7
## perm= 8
## perm= 9
## perm= 10
## perm= 11
## perm= 12
## perm= 13
## perm= 14
## perm= 15
## perm= 16
## perm= 17
## perm= 18
## perm= 19
## perm= 20
## perm= 21
## perm= 22
## perm= 23
## perm= 24
##
## Computing delta table
## 1
## 2
## 3
## 4
## 5
## 6
## 7
## 8
## 9
## 10
## 11
## 12
## 13
## 14
## 15
## 16
## 17
## 18
## 19
## 20
## 21
## 22
## 23
## 24
## 25
## 26
## 27
## 28
## 29
## 30
## 31
## 32
## 33
## 34
## 35
## 36
## 37
## 38
## 39
## 40
## 41
## 42
## 43
## 44
## 45
## 46
## 47
## 48
## 49
## 50
SAM函数返回值一个列表结构，可以自己用?SAM看看。差异表达基因的数据在siggenes.table中，也是一个列表结构：
str(samfit$siggenes.table)
## List of 5
## $ genes.up : chr [1:6748, 1:7] "265483_at" "248583_at" "253183_at" "252409_at" ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:7] "Gene ID" "Gene Name" "Score(d)" "Numerator(r)" ...
## $ genes.lo : chr [1:5341, 1:7] "259382_s_at" "248748_at" "249658_s_at" "266863_at" ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:7] "Gene ID" "Gene Name" "Score(d)" "Numerator(r)" ...
## $ color.ind.for.multi: NULL
## $ ngenes.up : int 6748
## $ ngenes.lo : int 5341
上调基因在siggenes.table的genes.up中，下调基因在genes.lo中。从上面的数据结构显示还可以看到差异表达基因的数量： ngenes.up和ngenes.lo。提取差异表达基因数据：
results.sam<-data.frame(rbind(samfit$siggenes.table$genes.up,samfit$siggenes.table$genes.lo),row.names=1,stringsAsFactors=FALSE)for(iin1:ncol(results.sam)) results.sam[,i]<-as.numeric(results.sam[,i])head(results.sam,2)
## Gene.Name Score.d. Numerator.r. Denominator.s.s0. Fold.Change
## 265483_at 14534 222.3 22.48 0.101 1.989
## 248583_at 2667 186.8 88.45 0.473 1.670
## q.value...
## 265483_at 0
## 248583_at 0
应用表达倍数进行筛选，有861个基因表达变化超过2倍（和前面简单t测验结果仅差1个，说明t测验还是可以的嘛！）：
results.sam<-results.sam[abs(log2(results.sam$Fold.Change))>=log2(2), ] ;nrow(results.sam)
## [1] 861
应用q值筛选，q<0.05只有10个，而q<0.1则有685个，选择筛选阈值也成了这种方法的一个问题：
#samr的q值表示方式为%，即5表示5%nrow(results.sam[results.sam$q.val<5,])
## [1] 10
nrow(results.sam[results.sam$q.val<10,])
## [1] 685
3.3Wilcoxon's signed-rank test
这个方法发表在 Liu, W.-m. et al, Analysis of high density expression microarrays with signed-rank call algorithms. Bioinformatics, 2002, 18, 1593-1599。R软件包simpleaffy的detection.p.val函数有实现，可以通过pairwise.comparison函数调用：
library(simpleaffy)#注意下面语句中的数据顺序sa.fit<-pairwise.comparison(eset.rma,"treatment",c("7d","0h"))
pairwise.comparison返回的数据为simpleaffy自定义的"PairComp"类型，提取数据要用它专门的函数：平均值用means函数获得，变化倍数（log2）用fc函数获得，t测验的p值用tt函数获得：
class(sa.fit)
## [1] "PairComp"
## attr(,"package")
## [1] "simpleaffy"
results.sa<-data.frame(means(sa.fit),fc(sa.fit),tt(sa.fit))#选择有表达的基因results.sa<-results.sa[present.probes,]head(results.sa,2)
## X7d X0h fc.sa.fit. tt.sa.fit.
## 244901_at 4.047 4.203 -0.1562 0.43982
## 244902_at 3.938 4.295 -0.3570 0.05824
colnames(results.sa)<-c("7d","0h","fold.change","p.val")head(results.sa,2)
## 7d 0h fold.change p.val
## 244901_at 4.047 4.203 -0.1562 0.43982
## 244902_at 3.938 4.295 -0.3570 0.05824
应用表达倍数筛选得到表达倍数超过2倍的基因数量有862个，应用p值筛选后得到562个差异表达基因：
results.sa<-results.sa[abs(results.sa$fold.change)>=log2(2), ];nrow(results.sa)
## [1] 862
results.sa<-results.sa[results.sa$p.val<0.05,];nrow(results.sa)
## [1] 562
这种方法在R软件包limma里面实现得最好。limma最初主要用于双色（双通道）芯片的处理，现在不仅支持单色芯片处理，新版还添加了对RNAseq数据的支持，很值得学习使用。安装方法同前面其他Bioconductor软件包的安装。载入limm软件包后可以用limmaUsersGuide()函数获取pdf格式的帮助文档。
limma的功能很多，这里只看看差异表达基因的分析流程，具体算法原理请参考limma在线帮助和这篇文献：Smyth G K. Linear models and empirical bayes methods for assessing differential expression in microarray experiments[J]. Statistical applications in genetics and molecular biology, 2004, 3(1): 3.
limma需要先产生一个design矩阵，用于描述RNA样品：
library(limma)treatment<-factor(pData(eset.rma)$treatment)design<-model.matrix(~0+treatment)colnames(design)<-c("C0h","T1h","T24h","T7d")design
## C0h T1h T24h T7d
## 1 1 0 0 0
## 2 1 0 0 0
## 3 0 1 0 0
## 4 0 1 0 0
## 5 0 0 1 0
## 6 0 0 1 0
## 7 0 0 0 1
## 8 0 0 0 1
## attr(,"assign")
## [1] 1 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$treatment
## [1] "contr.treatment"
可以看到：矩阵的每一行代表一张芯片，每一列代表一种RNA来源（或处理）。此外，你可能还需要另外一个矩阵，用来说明你要进行哪些样品间的对比分析：
contrast.matrix<-makeContrasts(T1h-C0h, T24h-C0h, T7d-C0h,levels=design)contrast.matrix
## Contrasts
## Levels T1h - C0h T24h - C0h T7d - C0h
## C0h -1 -1 -1
## T1h 1 0 0
## T24h 0 1 0
## T7d 0 0 1
下一步建立线性模型，并进行分组比较和p值校正：
fit<-lmFit(eset.rma[present.probes,], design)fit2<-contrasts.fit(fit, contrast.matrix)fit2<-eBayes(fit2)
先统计一下数量。可以看到：对于T7d-C0h比较组（coef=3），表达变化超过2倍（lfc参数）的基因数为842个，而变化超过2倍、p<0.05的基因总数为740个：
nrow(topTable(fit2,coef=3,adjust.method="fdr",lfc=1,number=30000))
## [1] 842
nrow(topTable(fit2,coef=3,adjust.method="fdr",p.value=0.05,lfc=1,number=30000))
## [1] 740
把toTable函数的返回结果存到其他变量就可以了，它是数据框类型数据，可以用write或write.csv函数保存到文件：
results.lim<-topTable(fit2,coef=3,adjust.method="fdr",p.value=0.05,lfc=1,number=30000)class(results.lim)
## [1] "data.frame"
head(results.lim)
## logFC AveExpr t P.Value adj.P.Val B
## 254818_at 6.215 10.363 41.38 7.304e-10 1.169e-05 11.538
## 254805_at 6.844 7.280 30.81 6.095e-09 4.878e-05 10.431
## 245998_at 2.778 10.011 25.44 2.411e-08 9.545e-05 9.528
## 265119_at 4.380 8.282 24.07 3.588e-08 9.545e-05 9.240
## 256114_at 4.461 7.668 23.92 3.745e-08 9.545e-05 9.208
## 265722_at -2.913 9.276 -23.91 3.760e-08 9.545e-05 9.205
为什么以上几种方法仅用表达倍数（2倍）筛选得到的数字不大一样？limma和直接计算的结果都是842个，而simpleaffy和SAM为862/861个。这是对eset信号值取对数和求平均值的先后导致的，limma先取对数再求平均值，而simpleaffy和SAM是先求平均值再取对数。
如Rank products方法，在R软件包RankProd里实现，方法文献为：Breitling R, Armengaud P, Amtmann A, et al. Rank products: a simple, yet powerful, new method to detect differentially regulated genes in replicated microarray experiments[J]. FEBS letters, 2004, 573(1): 83-92.
最后我们保存部分数据备下次使用：
#所有表达基因的名称write(present.probes,"genes.expressed.txt")#处理7天的差异表达基因write.csv(results.lim,"results.lim.7d.csv")#emat.rma.log2write.csv(emat.rma.log2[present.probes,],"emat.rma.log2.csv")
如果要全部结果：
results.lim.all<-topTable(fit2,coef=1:3,adjust.method="fdr",p.value=1,lfc=0,number=30000)head(results.lim.all,3)
## T1h...C0h T24h...C0h T7d...C0h AveExpr F P.Value
## 254818_at 0.34085 6.024 6.215 10.363 1048.0 6.704e-10
## 245998_at -0.13675 3.676 2.778 10.011 630.3 4.192e-09
## 265119_at -0.02536 6.061 4.380 8.282 579.6 5.671e-09
## adj.P.Val
## 254818_at 1.073e-05
## 245998_at 3.025e-05
## 265119_at 3.025e-05
## 仅保存logFC供后面的分析使用results.lim.all<-results.lim.all[,1:3]colnames(results.lim.all)<-c('T1h','T24h','T7d')head(results.lim.all,3)
## T1h T24h T7d
## 254818_at 0.34085 6.024 6.215
## 245998_at -0.13675 3.676 2.778
## 265119_at -0.02536 6.061 4.380
write.csv(results.lim.all,'results.lim.all.csv')
sessionInfo()
## R version 3.1.0 (2014-04-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
##
## locale:
## [1] LC_CTYPE=zh_CN.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=zh_CN.UTF-8 LC_COLLATE=zh_CN.UTF-8
## [5] LC_MONETARY=zh_CN.UTF-8 LC_MESSAGES=zh_CN.UTF-8
## [7] LC_PAPER=zh_CN.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=zh_CN.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel tcltk stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] limma_3.21.1 simpleaffy_2.41.0 gcrma_2.37.0
## [4] genefilter_1.47.0 samr_2.0 matrixStats_0.8.14
## [7] impute_1.39.0 ath1121501cdf_2.14.0 affy_1.43.0
## [10] Biobase_2.25.0 BiocGenerics_0.11.0 zblog_0.1.0
## [13] knitr_1.5
##
## loaded via a namespace (and not attached):
## [1] affyio_1.33.0 annotate_1.43.2 AnnotationDbi_1.27.3
## [4] BiocInstaller_1.15.2 Biostrings_2.33.3 DBI_0.2-7
## [7] evaluate_0.5.3 formatR_0.10 GenomeInfoDb_1.1.2
## [10] highr_0.3 IRanges_1.99.2 preprocessCore_1.27.0
## [13] R.methodsS3_1.6.1 RSQLite_0.11.4 S4Vectors_0.0.2
## [16] splines_3.1.0 stats4_3.1.0 stringr_0.6.2
## [19] survival_2.37-7 tools_3.1.0 XML_3.98-1.1
## [22] xtable_1.7-3 XVector_0.5.3 zlibbioc_1.11.1
版权声明：本文为博主原创文章，未经博主允许不得转载。
