生物技术

基因芯片(Affymetrix)分析:差异表达基因的提取

在提取差异表达的基因时,人们总是会有这两种考虑,一是不可漏过一个,二是不能错杀过多(在英语里称为false discovery rate(FDR)错误发现率)。常见的手段是使用多种统计学方法来分析同样一个结果,尽可能多的得到差异表达的基因,而排除那些假的信号。然而学习和使 用多种统计分析手段并不一定对于每一个生物学工作者都是非常容易的,这需要付出时间和努力。在这里,我们尽量多介绍几种常用的统计分析手段,并给出实践中 人们常常使用的组合,来帮助你更好的分析自己的数据。

现在常用的分析手段主要有:significance analysis of microarrays(SAM),CyberT和Rank products(RP)三种手段。其中CyberT是bioconductor当中最为常用的分析手段,因为它的算法完整地被limma库实现。但有研 究指出,使用SAM和RP算法相结合可能是最佳的方案。其实任何一种算法都是有局限性的,我们需要从根本上对算法有所了解,然后才能有针对性地选择合适的 算法。

SAM:Tusher VG, Tibshirani R, Chu G. Significance analysis of microarrays applied to the ionizing radiation response. Proc Natl Acad Sci USA 2001; 98:5116-21

CyberT: Baldi P, Long AD. A Bayesian framework for the analysis of microarray expression data: regularized t-test and statistical inferences of gene changes. Bioinformatics 2001; 17:509-19

RP: 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. FEBS lett 2004;573:83-92

1.  读取数据

本程序使用的具体数据不能透露,但是实验方案是:共6个实验,1个处理组(3个重复),1个对照组(3个重复)。大家在按照教程时,务必要按照自己的实验方案处理。

其实affy包里面就已经提供了分析方法,即rma()函数。

library(affy)   
library(tcltk)   
data.raw <- ReadAffy()   
pData(data.raw)$treatment <- rep(c(“N”,“T”),each=3)   
pData(data.raw)   
#                       sample treatment   
N-2-2_(ATH1-121501).CEL      1         N   
N-2-3_(ATH1-121501).CEL      2         N   
N-2-5_(ATH1-121501).CEL      3         N   
Y-2-3_(ATH1-121501).CEL      4         T   
Y-2-4_(ATH1-121501).CEL      5         T   
Y-2-5_(ATH1-121501).CEL      6         T  

2. 然后对数据进行初步处理。
rma()其实是一个综合函数,他分为了三个步骤,在数据的初步处理中已经提到过了。当然也可以用其他函数处理。

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.
## |####################|

3. 计算基因表达量

很简单,用一个exprs函数就可以从eset数据中提取出表达量,得到的数据类型是矩阵。但是应该注意rma的eset结果是经过对数变换的,而mas5的eset结果是原始信号强度。虽然表达量是用对数变换的信号值表示的,但是有些计算过程要用到未经变换的原始值,应该把它们都计算出来:

emat.rma.log2 <- exprs(eset.rma)   
emat.mas5.nologs <- exprs(eset.mas5)  

这里一定要注意:前者得到的表达量是经过log处理了的,而后者没有。

emat.rma.log2 emat.mas5.nologs head(emat.rma.log2,1)
N-2-2_(ATH1-121501).CEL N-2-3_(ATH1-121501).CEL
244901_at 1.780563 1.767244
N-2-5_(ATH1-121501).CEL Y-2-3_(ATH1-121501).CEL
244901_at 1.784543 1.759243
Y-2-4_(ATH1-121501).CEL Y-2-5_(ATH1-121501).CEL
244901_at 1.683675 2.039616
&gt; head(emat.mas5.nologs,1)
N-2-2_(ATH1-121501).CEL N-2-3_(ATH1-121501).CEL
244901_at 88.73547 131.0235
N-2-5_(ATH1-121501).CEL Y-2-3_(ATH1-121501).CEL
244901_at 73.61482 127.7071
Y-2-4_(ATH1-121501).CEL Y-2-5_(ATH1-121501).CEL
244901_at 54.30507 118.1789

下面仅使用rma的结果做演示。计算平均表达量和差异表达倍数

results.rma <- data.frame((emat.rma.log2[,c(1,4)] + emat.rma.log2[,c(2,5)] + emat.rma.log2[,c(3,6)])/3)   
> head(results.rma,2)   
          N.2.2_.ATH1.121501..CEL Y.2.3_.ATH1.121501..CEL   
244901_at                1.777450                1.827511   
244902_at                1.818097                1.840394   
> results.rma$T <- results.rma[,2] - results.rma[,1]

4. 选取“表达”基因

很多人可能不做这一步,认为没必要。但是理论上说,芯片可以检测的基因不一定在样品中都有表达,对于样本量较小的非“pooled”样品来说这是肯定的。把芯片上所有基因当成样本中的表达基因去做统计分析显然不合适。

选取“表达”基因的方法常见的有两种,一是使用genefilter软件包,另外一种是调用affy包的mas5calls()函数。使用 genefilter需要设定筛选阈值,不同的人可能有不同的标准,稍嫌随机,不够自动化,不介绍了。mas5calls方法使用探针水平数据(AffyBatch类型数据)进行处理,一般使用没经过预处理的芯片数据通用性强些,其他参数用默认就可以:

data.mas5calls <- mas5calls(data.raw)  

继续用exprs计算“表达”量,得到的数据只有三个值P/M/A。对于这三个值的具体解释可以用mas5calls查看帮助。P为present,A为absent,M为marginal(临界值)。

> data.mas5calls <- mas5calls(data.raw)   
Getting probe level data…   
Computing p-values   
Making P/M/A Calls   
> eset.mas5calls <- exprs(data.mas5calls)   
> head(eset.mas5calls)   
          N-2-2_(ATH1-121501).CEL N-2-3_(ATH1-121501).CEL   
244901_at “P”                     “P”                       
244902_at “A”                     “P”                                      
          N-2-5_(ATH1-121501).CEL Y-2-3_(ATH1-121501).CEL   
244901_at “A”                     “P”                       
244902_at “A”                     “P”                                    
          Y-2-4_(ATH1-121501).CEL Y-2-5_(ATH1-121501).CEL   
244901_at “A”                     “P”                       
244902_at “A”                     “M”                                       

下面我们把至少一个芯片中有表达的基因选出来:从 22810个中 挑选出13537 。相当于取并集。

> AP <- apply(eset.mas5calls, 1, function(x)any(x==“P”))   
> present.probes <- names(AP[AP])   
> paste(length(present.probes),“/”,length(AP))   
[1] “13537 / 22810″  
> results.present <- results.rma[present.probes,] 

5. 获取差异表达基因

生物学数据分析时的”差异”应该有两个意思,一是统计学上的差异,另外一个是生物学上的差异。一个基因在两个条件下的表达量分别有3个测量值:99,100,101 和 102,103,104。统计上两种条件下的基因表达数值是有差异的,后者比前者表达量要大。但生物学上有意义吗?未必。按平均值计算表达变化上升了3%,能产生什么样的生物学效应?这得看是什么基因了。所以差异表达基因的选取一般设置至少两个阈值:基因表达变化量和统计显著性量度(p值、q值等)。

5.1 t-test 算法vs SAM 算法

简单t-测验
这种方法不用太多的统计学知识,生物专业的人很容易想到,而且确实有不少人在用。 经常使用的筛选阈值是表达量变化超过2倍,即|log2(fc)|>=log(2)。先简单看看有没有

apply(abs(results.present[,3]), 2, max)

或者 max(results.present$T)

apply是一个很有用的函数,它对数据按某个维度批量应用一个函数进行计算。第一个参数为向量或矩阵(或者是能转成向量或矩阵的数据,如数据框),第三个参数表示要使用的函数,第二个参数为应用的维度。上面语句的意思是对数据 abs(results.present[,5:7]) 按列(第二维)使用统计函数max(计算最大值)。 表达变化超过2倍的基因共有546个:

> sum(abs(results.present[,"T"])>=log2(2))   
[1] 546  

选出这546个基因:

t-test 测验,并选出p<0.05的差异表达基因:

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
> results.st <- results.present[abs(results.present$T)>=log2(2),]   
> sel.genes <- row.names(results.st)   
> head(sel.genes,5)   
[1]244939_at”   “244940_at”   “244979_at”   “244981_at”   “244987_s_at”  
> head(emat.rma.log2,5)   
          N-2-2_(ATH1-121501).CEL N-2-3_(ATH1-121501).CEL   
244901_at                1.780563                1.767244   
244902_at                1.812207                1.925442   
244903_at                4.079256                3.496978   
244904_at                2.845130                3.001135   
244905_at                1.699556                1.447394   
          N-2-5_(ATH1-121501).CEL Y-2-3_(ATH1-121501).CEL   
244901_at                1.784543                1.759243   
244902_at                1.716642                1.672951   
244903_at                3.709451                3.694250   
244904_at                3.011112                2.629638   
244905_at                1.469067                1.196991   
          Y-2-4_(ATH1-121501).CEL Y-2-5_(ATH1-121501).CEL   
244901_at                1.683675                2.039616   
244902_at                1.612765                2.235466   
244903_at                4.212147                3.296236   
244904_at                3.200534                2.997960   
244905_at                1.231213                1.685121   
> p.value <- apply(emat.rma.log2[sel.genes,], 1, function(x){t.test(x[1:3], x[4:6])$p.value})   
> results.st$p.value <- p.value   
> names(results.st)   
[1] “N.2.2_.ATH1.121501..CEL” “Y.2.3_.ATH1.121501..CEL[3] “T”                       “p.value> results.st <- results.st[p.value<0.05,]   
> nrow(results.st)   
[1] 208

通过计算,简单t测验方法得到208个表达倍数变化超过2倍的差异表达基因。

5.2 SAM 算法求差异表达基因

这种方法流行过一段时间,但由于FDR(错误检出率)控制太差,现在基本不用了。

要用也不复杂。但是注意SAM函数使用的emat表达数据是present.probes筛选出来的“表达”基因子集,如果你用没有经过筛选的数据,得到的结果会差别很大,不信可以自己试试(这点可能也是这种方法的毛病之一):

library(samr)   
samfit <- SAM(emat.rma.nologs[present.probes,c(1,2,3,4,5,6)], c(1,1,1,2,2,2),   
 resp.type=“Two class unpaired”, genenames=present.probes)  

SAM函数返回值一个列表结构,可以自己用?SAM看看。差异表达基因的数据在siggenes.table中,也是一个列表结构:

str(samfit$siggenes.table)

上调基因在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(i in 1:ncol(results.sam)) results.sam[,i] <- as.numeric(results.sam[,i])   
head(results.sam, 2)  

应用表达倍数进行筛选:

results.sam <- results.sam[abs(log2(results.sam$Fold.Change))>=log2(2), ] ; nrow(results.sam)

应用q值筛选,q<0.05只有23个,而q<0.1则有593个,选择筛选阈值也成了这种方法的一个问题:

1
2
3
4
5
#samr的q值表示方式为%,即5表示5%   
nrow(results.sam[results.sam$q.val<5,])   
## [1] 23   
nrow(results.sam[results.sam$q.val<10,])   
## [1] 593

5.3 CyberT 算法

CyberT模型将标准及信号强度的关系使用线性模型进一步强化,并推出了Bioconductor下完整的limma库,该库说明文件清晰,教程完善,故而使用率极高。CyberT模型使用的是基于贝叶斯估计调整后的t-test来确定差异表达的基因,准确率较高。并且有研究指出,它的计算结果要强于SAM算法。虽然它们思路很类似,但是CyberT和limma其实还是有一点点差别的。在t-test公式上,它们都是(区别于SAM的),其中代表倍数平均值,代表标准差的平均值,但是的计算方式不完全相同。不类CyberT,limma在计算标准差时考虑的是全部的基因,而非排序后相近的基因,所以,limma库不适合使用MAS5算法生成的数据,这就使得我们只能倾向于rma算法。但是也正因为如此,所以limma不象CyberT那样局限于两组数据,limma可以依照实验构建对比矩阵,计算多级复杂设计的实验。

用 limma 实现 寻找差异基因放在文章的最后。

5.4 Rank products(RP)算法

注意:对于RP来说,只能应对小于两组的实验分组,如果多组的话,得分别进行比对。

使用文档可以参看:http://www.bioconductor.org/packages/2.12/bioc/vignettes/RankProd/inst/doc/RankProd.pdf

方法文献为: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.

RankProd使用的是非参数结果排序法来判断两组对照实验之间基因是否上调或者下调。RP就是测量中某基因的几何平均值。具体RP决定差异表达基因的步骤如下:

1. 生成长度为n的k排序的p排列
2. 计算在p排列中n个基因的RP值
3. 计数在排列中有多少基因的RP值是小于等于对照组的RP值的。设定该计数为c。
4. 通过公式计算 RP期待值。
5. 通过公式计算假阳性的可能性,依照pfp值来决定是否为差异表达。其中rank(g)是RP升序排列的n个基因中基因g的排序。

source(“http://bioconductor.org/biocLite.R”)   
biocLite(“RankProd”)  
library(RankProd)

这里只以我的数据为例进行说明

数据承接之前的:

eset.rma <- rma(data.raw)   
eset<-exprs(eset.rma )  

对实验进行设计:

eset.cl <- c(0,0,0,1,1,1)   
eset.origin <- c(1,1,1,1,1,1)   

对比:

> eset<-exprs(eset.rma )    
> eset.cl <- c(0,0,0,1,1,1)   
> eset.origin <- c(1,1,1,1,1,1)     
> gnames<-geneNames(data.raw)   
> head(gnames,5)   
[1] “244901_at” “244902_at” “244903_at” “244904_at” “244905_at”  
> eset.sub<-eset[,which(eset.origin==1)]   
> eset.cl.sub<-eset.cl[which(eset.origin==1)]   
> eset.origin.sub<-eset.origin[which(eset.origin==1)]   
> library(RankProd)   
> RP.T<-RPadvance(eset.sub,eset.cl.sub,eset.origin.sub,num.perm=100,logged=TRUE, na.rm=FALSE, gene.names=gnames, plot=FALSE, rand=123)  

生成结果报表:

#生成结果报表   
>topGene(RP.T,cutoff=0.05,method=“pfp”,logged=TRUE,logbase=2,gene.names=gnames)   
Table1: Genes called significant under class1 < class2    
  
Table2: Genes called significant under class1 > class2    
  
$Table1  
            gene.index  RP/Rsum FC:(class1/class2)    pfp P.value   
247573_at         2673   6.4846             0.0880 0.0000   0e+00   
254447_at         9547  12.1243             0.1131 0.0000   0e+00   
249032_at         4132  17.0109             0.1356 0.0000   0e+00   
264202_at        19302  18.3599             0.1158 0.0025   0e+00   
252265_at         7365  19.7263             0.1343 0.0020   0e+00   
256787_at        11887  21.4240             0.1305 0.0017   0e+00   
$Table2  
                    gene.index  RP/Rsum FC:(class1/class2)    pfp P.value   
264146_at                19246  14.2418             6.2627 0.0000   0e+00   
256237_at                11337  28.5948             4.1473 0.0000   0e+00   
255822_at                10922  29.0446             4.6106 0.0000   0e+00   
261221_at                16321  30.0674             4.2118 0.0000   0e+00 

5.5 Wilcoxon’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(“T”, “N”))   
#T 是代表处理组;N 代表对照组;数据之前做过介绍  

pairwise.comparison返回的数据为simpleaffy自定义的”PairComp”类型,提取数据要用它专门的函数:平均值用means函数获得,变化倍数(log2)用fc函数获得,t测验的p值用tt函数获得:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
> library(simpleaffy)   
> sa.fit <- pairwise.comparison(eset.rma, “treatment”, c(“T”, “N”))   
> 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)   
                 T        N fc.sa.fit. tt.sa.fit.   
244901_at 1.835792 1.777469 0.05832259  0.6682087   
244902_at 1.868711 1.820625 0.04808520  0.9195176   
> colnames(results.sa) <- c(“T”, “N”, “fold.change, “p.val)   
> head(results.sa, 2)   
                 T        N fold.change     p.val   
244901_at 1.835792 1.777469  0.05832259 0.6682087   
244902_at 1.868711 1.820625  0.04808520 0.9195176   
> results.sa <- results.sa[abs(results.sa$fold.change)>=log2(2), ]; nrow(results.sa)   
[1] 638   
> results.sa <- results.sa[results.sa$p.val<0.05,]; nrow(results.sa)   
[1] 281

5.6 用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样品:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
> treatment <- factor(pData(data.raw)$treatment)   
> treatment   
[1] N N N T T T   
Levels: N T   
> design <- model.matrix(~0 + treatment)   
> colnames(design) <- c(“N”,“T”)   
>  #可以看到:矩阵的每.行代表一张芯片,每一列代表一种RNA来源(或处理)。此外,你可能还需要另外一个矩阵,用来说明你要进行哪些样品间的.确治   
> cntrast.matrix <- makeContrasts(T-N,levels=design)   
> cntrast.matrix   
      Contrasts   
Levels T - N   
     N    -1   
     T     1   
> #下一步建立线性模型,并进行分组比.?校   
> #下一步建立线性模型,并进行分组比较和p值校正:   
> fit <- lmFit(eset.rma[present.probes,], design)   
> fit2 <- contrasts.fit(fit, cntrast.matrix)   
> fit2 <- eBayes(fit2)   
># 参数解释: coef 代表第几组对比,在函数makeContrasts体现,我们这里就只有一组;   
># lfc 就是 fold change    
> nrow(topTable(fit2, coef=1, adjust.method=“fdr”, lfc=1, number=30000))   
[1] 546   
> nrow(topTable(fit2, coef=1, adjust.method=“fdr”, p.value=0.05, lfc=1, number=30000))   
[1] 220  #

为什么以上几种方法仅用表达倍数(2倍)筛选得到的数字不大一样?这是对eset信号值取对数和求平均值的先后导致的,limma先取对数再求平均值,而simpleaffy和SAM是先求平均值再取对数。

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
> sessionInfo()   
R version 2.15.2 (2012-10-26)   
Platform: x86_64-unknown-linux-gnu (64-bit)   
 
locale:   
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C                 
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8       
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8      
 [7] LC_PAPER=C                 LC_NAME=C                    
 [9] LC_ADDRESS=C               LC_TELEPHONE=C               
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C          
 
attached base packages:   
[1] tcltk     stats     graphics  grDevices utils     datasets  methods     
[8] base        
 
other attached packages:   
 [1] simpleaffy_2.34.0    gcrma_2.30.0         genefilter_1.40.0      
 [4] limma_3.14.4         ath1121501cdf_2.11.0 AnnotationDbi_1.20.7   
 [7] affy_1.36.1          Biobase_2.18.0       BiocGenerics_0.4.0     
[10] BiocInstaller_1.8.3    
 
loaded via a namespace (and not attached):   
 [1] affyio_1.26.0         annotate_1.36.0       Biostrings_2.26.3       
 [4] DBI_0.2-7             IRanges_1.16.6        parallel_2.15.2         
 [7] preprocessCore_1.20.0 RSQLite_0.11.4        splines_2.15.2          
[10] stats4_2.15.2         survival_2.36-14      tools_2.15.2            
[13] XML_3.98-1.1          xtable_1.7-1          zlibbioc_1.4.0

发表评论