【图文】基因芯片分析_百度文库
两大类热门资源免费畅读
续费一年阅读会员，立省24元！
评价文档：
基因芯片分析
上传于|0|0|文档简介
&&基因芯片分析
大小：9.66MB
登录百度文库，专享文档复制特权，财富值每天免费拿！
你可能喜欢13459人阅读
基因芯片分析（5）
“差异”是个统计学概念，获取差异表达基因就要用统计方法，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&
# 简化一下名称，设置pData
sampleNames(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
使用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
## |####################|
1 计算基因表达量
很简单，用一个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)
## 244901_at 4.04 4.348 4.048 4.052 4.019 3.962 4.03 4.062
head(emat.mas5.nologs, 1)
## 244901_at 39.79 94.94 57.35 60.23 61.08 55.57 53.95 85.98
emat.rma.nologs &- 2^emat.rma.log2
emat.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)
## 244901_at 4.194 4.050 3.991 4.046 -0.7 -0.1481
## 244902_at 4.293 4.159 4.061 3.937 -0.6 -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&0
subset.data &- results.rma[subset.logic,]
要注意的是逻辑向量的长度要和相应维度的数据长度一致，逻辑向量中为TRUE的就保留，FALSE的就丢弃：
length(subset.logic); nrow(results.rma)
## [1] 22810
head(subset.logic)
## [1] FALSE FALSE FALSE FALSE
2 选取“表达”基因
很多人可能不做这一步，认为没必要。但是理论上说，芯片可以检测的基因不一定在样品中都有表达，对于样本量较小的非“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)
## 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 获取差异表达基因
生物学数据分析时的&差异&应该有两个意思，一是统计学上的差异，另外一个是生物学上的差异。一个基因在两个条件下的表达量分别有3个测量值：99,100,101 和 102,103,104。统计上两种条件下的基因表达数值是有差异的，后者比前者表达量要大。但生物学上有意义吗？未必。按平均值计算表达变化上升了3%，能产生什么样的生物学效应？这得看是什么基因了。所以差异表达基因的选取一般设置至少两个阈值：基因表达变化量和统计显著性量度（p值、q值等）。
3.1 简单t-测验
这种方法不用太多的统计学知识，生物专业的人很容易想到，而且确实有不少人在用。经常使用的筛选阈值是表达量变化超过2倍，即|log2(fc)|&=log(2)。先简单看看有没有：
apply(abs(results.present[,5:7]), 2, max)
fc.1h fc.24h
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.value
names(results.st)
## [1] &S1&
## [8] &p.value&
results.st &- results.st[, c(1,4,7,8)]
results.st &- results.st[p.value&0.05,]
head(results.st, 2)
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.2 SAM（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
SAM函数返回值一个列表结构，可以自己用?SAM看看。差异表达基因的数据在siggenes.table中，也是一个列表结构：
str(samfit$siggenes.table)
## List of 5
$ genes.up
: chr [1:] &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:] &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(i in 1: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
## 248583_at
q.value...
## 265483_at
## 248583_at
应用表达倍数进行筛选，有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,])
nrow(results.sam[results.sam$q.val&10,])
## [1] 685
3.3 Wilcoxon's signed-rank test
这个方法发表在 Liu, W.-m. et al, Analysis of high density expression microarrays with signed-rank call algorithms. Bioinformatics, 93-1599。R软件包simpleaffy的detection.p.val函数有实现，可以通过parison函数调用：
library(simpleaffy)
#注意下面语句中的数据顺序
sa.fit &- parison(eset.rma, &treatment&, c(&7d&, &0h&))
<parison返回的数据为simpleaffy自定义的&PairComp&类型，提取数据要用它专门的函数：平均&#20540;用means函数获得，变化倍数（log2）用fc函数获得，t测验的p&#20540;用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)
X0h fc.sa.fit. tt.sa.fit.
## 244901_at 4.047 4.203
## 244902_at 3.938 4.295
colnames(results.sa) &- c(&7d&, &0h&, &fold.change&, &p.val&)
head(results.sa, 2)
0h fold.change
## 244901_at 4.047 4.203
## 244902_at 3.938 4.295
应用表达倍数筛选得到表达倍数超过2倍的基因数量有862个，应用p&#20540;筛选后得到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
3.4 Moderated T statistic
这种方法在R软件包limma里面实现得最好。limma最初主要用于双色（双通道）芯片的处理，现在不仅支持单色芯片处理，新版还添加了对RNAseq数据的支持，很&#20540;得学习使用。安装方法同前面其他Bioconductor软件包的安装。载入limm软件包后可以用limmaUsersGuide()函数获取pdf&#26684;式的帮助文档。
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, ): 3.
limma需要先产生一个design矩阵，用于描述RNA样品：
library(limma)
treatment &- factor(pData(eset.rma)$treatment)
design &- model.matrix(~ 0 &#43; treatment)
colnames(design) &- c(&C0h&, &T1h&, &T24h&, &T7d&)
C0h T1h T24h T7d
## 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
## Levels T1h - C0h T24h - C0h T7d - C0h
下一步建立线性模型，并进行分组比较和p&#20540;校正：
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
P.Value adj.P.Val
## 254818_at
41.38 7.304e-10 1.169e-05 11.538
## 254805_at
30.81 6.095e-09 4.878e-05 10.431
## 245998_at
25.44 2.411e-08 9.545e-05
## 265119_at
24.07 3.588e-08 9.545e-05
## 256114_at
23.92 3.745e-08 9.545e-05
## 265722_at -2.913
9.276 -23.91 3.760e-08 9.545e-05
为什么以上几种方法仅用表达倍数（2倍）筛选得到的数字不大一样？limma和直接计算的结果都是842个，而simpleaffy和SAM为862/861个。这是对eset信号&#20540;取对数和求平均&#20540;的先后导致的，limma先取对数再求平均&#20540;，而simpleaffy和SAM是先求平均&#20540;再取对数。
3.5 其他方法：
如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, ): 83-92.
最后我们保存部分数据备下次使用：
#所有表达基因的名称
write(present.probes, &genes.expressed.txt&)
#处理7天的差异表达基因
write.csv(results.lim, &results.lim.7d.csv&)
#emat.rma.log2
write.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
## 254818_at
10.363 .704e-10
## 245998_at
630.3 4.192e-09
## 265119_at
579.6 5.671e-09
## 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)
## 254818_at
## 245998_at -0. 2.778
## 265119_at -0. 4.380
write.csv(results.lim.all, 'results.lim.all.csv')
4 Session Info
sessionInfo()
## R version 3.1.0 ()
## 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
[9] LC_ADDRESS=C
LC_TELEPHONE=C
## [11] LC_MEASUREMENT=zh_CN.UTF-8 LC_IDENTIFICATION=C
## attached base packages:
## [1] parallel
grDevices utils
## [8] methods
## other attached packages:
[1] limma_3.21.1
simpleaffy_2.41.0
gcrma_2.37.0
[4] genefilter_1.47.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
[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
&&相关文章推荐
参考知识库
* 以上用户言论只代表其个人观点，不代表CSDN网站的观点或立场
访问：187588次
积分：2684
积分：2684
排名：第12380名
原创：81篇
评论：34条基因表达谱芯片数据分析 - 广州济恒医药科技有限公司
基因表达谱芯片数据分析
&基因表达谱芯片数据分析及其Bioconductor实现&1.表达谱芯片及其应用表达谱DNA芯片（DNA microarrays for gene expression profiles）是指将大量DNA片段或寡核苷酸固定在玻璃、硅、塑料等硬质载体上制备成基因芯片，待测样品中的mRNA被提取后，通过逆转录获得cDNA，并在此过程中标记荧光，然后与包含上千个基因的DNA芯片进行杂交反应30min~20h后，将芯片上未发生结合反应的片段洗去，再对玻片进行激光共聚焦扫描，测定芯片上个点的荧光强度，从而推算出待测样品中各种基因的表达水平。用于研究基因表达的芯片可以有两种：① cDNA芯片；②寡核苷酸芯片。cDNA芯片技术及载有较长片段的寡核苷酸芯片采用双色荧光系统：目前常用 Cy3一dUTP（绿色）标记对照组mRNA，Cy5一dUTP（红色）标记样品组 mRNA[1]。用不同波长的荧光扫描芯片，将扫描所得每一点荧光信号值自动输入计算机并进行信息处理，给出每个点在不同波长下的荧光强度值及其比值（ratio值），同时计算机还给出直观的显色图。在样品中呈高表达的基因其杂交点呈红色，相反，在对照组中高表达的基因其杂交点呈绿色，在两组中表达水平相当的显黄色，这些信号就代表了样品中基因的转录表达情况[2]。基因芯片因具有高效率，高通量、高精度以及能平行对照研究等特点，被迅速应用于动、植物和人类基因的研究领域，如病原微生物毒力相关基因的。基因表达谱可直接检测mRNA的种类及丰度，可以同时分析上万个基因的表达变化，来揭示基因之间表达变化的相互关系。表达谱芯片可用于研究：①同一个体在同一时间里，不同基因的表达差异。芯片上固定的已知序列的cDNA或寡聚核苷酸最多可以达到30 000多个序列，与人类全基因组基因数相当，所以基因芯片一次反应几乎就能够分析整个人的基因[3]。②同一个体在不同时间里，相同基因的表达差异。③不同个体的相同基因表达上的差异。利用基因芯片可以分析多个样本，同时筛选不同样本（如肿瘤组织、癌前病变和正常组织）之间差异表达的基因，这样可以避免了芯片间的变异造成的误差[4]。张辛燕[5]等将 512个人癌基因和抑癌基因的cDNA用点样仪点在特制玻片上制成表达谱芯片，对正常人卵巢组织及卵巢癌组织基因表达的差异性进行比较研究，结果发现在卵巢癌组织中下调的基因有23个，上调的基因有15个，初步筛选出了卵巢癌相关基因。Lowe[6]等利用胰腺癌、问充质细胞癌等组织的cDNA制备基因芯片，筛选到胰腺癌细胞中高表达的基因，为医疗诊断、病理研究及新药设计奠定基础。&2. 表达谱芯片的数据处理技术2.1 探针水平数据（probe-level data）的获得提取生物样品的mRNA并反转录成cDNA，同时用荧光素或同位素标记。在液相中与基因芯片上的探针杂交，经洗膜后用图像扫描仪捕获芯片上的荧光或同位素信号[7]，由此获得的图像就是基因芯片的原始数据（raw data），也叫探针水平数据。获取探针水平的数据是芯片数据处理的第一步，然后需要对其进行预处理（pre-processing），以获得基因表达数据（gene expression data）。基因表达数据是芯片数据处理的基础。2.2 预处理2.2.1背景（background）处理背景处理即过滤芯片杂交信号中属于非特异性的背景噪音部分。一般以图像处理软件对芯片划格后，每个杂交点周围区域各像素吸光度的平均值作为背景。但此法存在芯片不同区域背景扣减不均匀的缺点，同时会使1％～5％[7]的点产生无意义的负值。也可利用芯片最低信号强度的点（代表非特异性的样本与探针结合值）或综合整个芯片非杂交点背景所得的平均值做为背景[8]。Brown[8]等提出利用整个芯片杂交点外的平均吸光度值作为背景的best-fit方法，使该问题得到较好的解决，并有效地提高了处理数据的质量。背景处理之后，我们可以将芯片数据放入一个矩阵中：&M = &&其中，各字母的意义如下：N：条件数；G：基因数目（一般情况下，G&&N）；行向量mi=(mi1,mi2,&,miN)表示基因i在N个条件下的表达水平（这里指绝对表达水平，亦即荧光强度值）；列向量mj=(m1j,m2j,&,mGj)表示在第j个条件下各基因的表达水平（即一张芯片的数据）；元素mij表示第基因i在第j个条件下（绝对）基因表达数据。m可以是R（红色，Cy5，代表样品组）。也可以是G（绿色，Cy3,代表对照组）。2.2.2数据清洗（data cleaning）经过背景校正后的芯片数据中可能会产生负值，显然负值是没有生物学意义的。数据集中还可能包括一些单个异常大（或小）的峰（谷）信号，它们被认为是随机噪声。另外，对于负值和噪声信号，通常的处理方法就是将其去除。然而，数据的缺失（除了上述原因会造成数据缺失以外，扫描的过程中也可能会产生缺失）对后续的统计分析（尤其是层式聚类和主成分分析）有致命的影响。所以对数据的删除，通常是删去所在的列向量或行向量。一个比较常用的做法是，事先定义个阈值M。若行（列）向量中的缺失数据量达到阈值M，则删去该向量。若未达到M，有两种方法处理，一是以0或者用基因表达谱中的平均值或中值代替，另一个是分析基因表达谱的模式，从中得到相邻数据点之间的关系，据此利用相邻数据点估算得到缺失值（类似于插值）。2.2.3归一化（normalization）经过背景处理和数据清洗处理后的修正值反映了基因表达的水平[9]。然而在芯片试验中，各个芯片的绝对光密度值是不一样的，在比较各个试验结果之前必需将其归一化（normalization，也称作标准化）。在同一块芯片上杂交的、由不同荧光分子标记的两个样品间的数据，也需归一化。常用的标准化方法有&看家基因法&、基于总光密度的方法、回归方法、比率统计法 [10]等。⑴ &看家基因（house-keeping gene）&法此法最为常用，可以用于几张芯片的数据归一化。它预先选择一组表达水平不变的看家基因，计算出这组基因平均ratio值为1时标准化系数，然后将其应用于全部的数据以达到归一化的目的。但是目前尚未找到理想的看家基因[11]，另外此前有研究表明，所谓&看家基因&在不同实验条件下其表达水平同样发生变化[12]。⑵基于总光密度的方法[13]此方法用于标准化同一块芯片上杂交的两种样品，它假设两批待标记的mRNA的量相同；相对于对照组样品，实验组的表达应既有上调也有下调。而且，扫描所得的所有Cy5和Cy3荧光分子的光密度值是相同的。据此计算出一个标准化系数，用以重新计算芯片上每个基因的光密度。⑶回归的方法[13]此方法用于标准化同一块芯片上杂交的两种样品。如果mRNA来自紧密相关的样品，那么大部分基因的表达水平是相近的。这样，在以Cy5和Cy3为坐标的散点图上，这些基因应呈一直线。如果两批样品的标记和检测效率相同，则直线的斜率也是惟一的。那么，标准化这些数据就等同于用回归的方法计算其最适斜率。但在实际试验中，光密度值常为非线性，此时应该使用局部回归方法，如LOWESS（1ocally weighted scatterplot smoothing）回归法。⑷比率统计法[13]此方法用于标准化同一块芯片上杂交的两种样品，并且建立于以下的假设之上：在近似的两个样品中，虽然基因有上调和下调，但一些基本的基因（如管家基因）的表达量是近似相同的。由此得出一个近似概率密度公式：比率T =R /G（R 和G分别是芯片上第K个点的红光和绿光的强度），经过迭代算法处理得到一个平均表达比率及其可信限，用于数据的标准化计算。2.3基因表达数据经过预处理，探针水平数据转变为基因表达数据。为了便于应用一些统计和数学术语，基因表达数据仍采用矩阵形式。随着生物学进入后基因组学时代，类似芯片数据这样的的非序列生物数据几乎呈指数形式膨胀。这些生物数据往往维数高，具有异质性和网络性，传统的分析方法已不能胜任，发展优秀的算法分析生物数据成为生物学研究的瓶颈。数据挖掘等决策支持技术因其在大规模数据处理方面的卓越能力而在其中占据越来越重要的地位[14]。数据挖掘也被称为数据库知识发现（knowledge discovery in database，KDD），是从数据库中识别出有效的、新颖的、潜在有用的并且最终可理解的、模式的非平凡过程[15]。迄今还没有一套完整、统一的数据挖掘理论体系来指导如何获取有用信息[16]。2.4差异表达基因（differentially expressed genes,DEGs）筛选用于检测基因表达水平的DNA微阵列的应用之一是比较实验，目的是比较两个条件下的基因差异表达，从中识别出与条件相关的特异性基因或显著差异表达的基因。2.4.1倍数变化法（fold chang, FC）比较两个各不同生物样本时，可根据ratio值来筛选，一般认为ratio值在0.5&2.0范围内的基因不存在显著表达差异，该范围之外则认为表达有显著差异。当然，上述范围需要根据不同实验条件作调整。FC法的优点是简单直观，需要的芯片量少，节约研究成本[16]；缺点是结论过于简单，其阈值的划分主观性较强、缺乏生物学和统计学支持，尤其对于分析样本中的低拷贝或高拷贝转录子，容易产生假阳性和假阴性问题[12]。一般而言，FC法可用于对于预实验和实验初筛。2.4.2&参数分析（parameter analysis）⑴ t-检验（t-test）t 检验可用于两个生物条件下多个重复样本的差异表达基因的筛选。当 t 超过根据可信度选择得标准时，比较的两样本被认为存在着差异。受样本量和成本的限制，研究者提出了调节性 t 检验（regulated t-test）。它根据在基因表达水平和变异之间存在着相互关系，相似的基因表达水平有着相似的变异这个经验，应用贝叶斯条件概率统计方法，通过检测同一张芯片上其他临近基因表达水平，理论上可对任何基因的变异程度估计进行弥补。调节性 t 检验法对基因表达的标准差估计优于一般t检验和FC法[20]。⑵ F 检验F 检验又称变异数分析或方差分析（analysis of variance, ANOVA）。F 检验适用于多个生物条件下DEGs的检测，它检验两个或多个样本均数的差异是否有统计学意义。方差分析需要参照实验设计，参照样本常用多种细胞的mRNA混合而成，由于所有的细胞同时表达基因众多，结果低表达基因在样本混合后就被稀释而减少了参照样本的代表性，因此，增加参照样本的细胞不会提高参照样本的代表性。方差分析的缺点在于虽然能计算出那些基因有统计差异，但是他没有对那些组之间有统计差异进行区分。如果相区分组间的统计差异，则需要使用均值间的两两比较（post-hoc comparisons）检验，该检验是对方差分析后的基因进行下一水平更细节的分析[15]。⑶回归分析（regression analysis）基因表达谱的回归分析可以处理多个基因变量间线性依存关系，研究者提出了&使用回归分析的基因表达谱数据&。Li 等[22]使用互变量（Cox）回归方法分析基因表达谱数据，用于患者的生存率预判；Huang 等[23]将线性回归方法应用于肿瘤的分类研究中。2.4.3非参数分析（nonparameter analysis）由于噪声的存在，通过数据转换后微阵列数据可能仍然不呈正态分布，因此使用参数分析法可能有风险。非参数检验的优点在于不必假设数据满足特殊的正态分布，尽管其对数据进行筛选有些粗放，而且其对表达数据分析的敏感性不如参数分，但是仍然可行。常用的基因表达谱数据分析的非参数检验方法有：传统的非参数t-检验（nonparametric t-test）[24]、Wilcoxon秩和检验（Wilcoxon rank sun test）[24]和新的非参数法如检验贝叶斯法（empirical Bayes method）[25]、芯片显著性分析法（signifcance analysis of micorarray,SAM）[26]、混合模型法（the mixture model method,MMM）[27]等。2.4.4假表达谱（pseudo profile）假表达谱常用于鉴别基因的某一特定性为。比如要鉴别在肺癌中高表达而在正常肺组织中和其他肿瘤组织中低表达的基因，就可以先假设具有这样一个假表达谱，然后在实际芯片数据中去寻找与其相吻合的基因[7]。关于DEGs的检测，目前尚无统一性标准，芯片后验证性实验（RT-PCR、荧光定量 RT-PCR、Northern等）是确定样本基因差异表达的黄金标准。2.5基因芯片数据分析的非监督方法在基因表达谱中找出差异表达基因只是对表达谱数据进行统计学分析第一步，通过建立共调控网络，发掘未知和已知基因功能才是芯片实验的最终目的。前者可以看成是基因表达的单基因水平分析，后者则为基因与蛋白质网络分析。根据对所研究的基因表达规律和实验分组是否了解，可将分析方法分为监督的（supervised）和非监督导的（unsupervised）。前者根据特定样本或基因的已知生物学信息对表达谱建立分类器，进而对各基因进行功能分类和预测，后者则通过计算和比较表达谱各基因统计学距离，聚类&相似性&样本或基因。两者都假设功能相似的基因其表达谱也是相似的，但Zhou[19]等认为，一些相似功能的基因并不总是表现相似的表达谱，针对此他们提出了&过渡共表达基因&概念及相应的数学模型鉴定表达谱中此类基因。2.5.1非监督的分析方法概述芯片数据统计分析的非指导的方法即聚类分析（cluster analysis），在目前最为常用。聚类分析是研究事物分类的一种方法，是在事物分类面貌尚不清楚的情况下研究事物的分类，其原理是直接比较样本中各指标之间的性质，将性质相近的归为一类，性质差别较大的归在另一类。统计学上通过计算相似距离（similarity distances）来比较数据，常用相关系数或欧氏距离表示。2.5.2&非监督分析中的数据降维（dimension reduction）在芯片数据中，有些数据并未提供有显著意义的信息，反而会给数据分析带来不必要的复杂。理想情况下，经过数据降维处理后，剩余数据即为非冗余数据（non-redundant data），不同组间的数据提供的信息是互相独立的。数据降维技术也分为监督的方法和非监督的方法。非监督分析中的数据降维主要是指删除不提供信息的数据。如果某一基因在不同条件下的表达水平相同，则它对区分这些不同条件没有任何作用，该基因所提供的数据即为冗余数据。为了去除冗余数据，可将冗余的数据整合到一个新的杂合分组中。主成分分析可很好的完成这一任务。2.5.3&非监督分析的各种技术简介⑴系统聚类（hierarchical clustering）[10]系统聚类根据聚类的方式分为凝聚法（agglomerative approach）和分裂法（divisive approach）。a. 凝聚法按照从下到上的方式对个体进行聚类，初始每个个体从各为一类、按照一定的规则进行逐步合并，直到所有个体都归为一类或达到预定的终止条件。凝聚法因类问相似性的度量方法的不同而又有所差异。b. 分裂法按照从上到下的方式对个体进行聚类，初始所有个体为一类，然后按照一定规则逐渐分裂，直到每个个体形成一类或满足某个特定的结束条件，如达到预定的类数或两个最邻近的类之间的距离超过某预定值。系统聚类方法简单，但有时在选择分裂点或合并点时存在困难。一旦将一组个体分裂或合并，后续的类将在新类的基础上产生，而不能取消己经完成的分裂或合并，也不能在类问对个体进行调整。系统聚类不适于分析基因表达谱复杂的数据[9]。⑵分割聚类(partitioning methods)[10][11]对于一个给定的基因芯片矩阵，分割算法将把观察个体分为预定的几部分，使得对个体的分割达到最优的客观标准，即类内个体间的相似性达到最大，而类间个体间的相似性达到最小。最常用的分割算法为k-means法和k-medoids法。a. k-means法把n个观察个体分成k个类，使类内的相似性高，而类间的相似性低。类的相似性用类内观察个体的均值来度量，此均值被视为类的重心。通过计算新形成的k类的类均数，达到目标函数收敛。具体步骤如下：所有数据随机分入k个簇中，每个簇的平均向量用于计算各簇间的距离。然后用迭代方法计算簇间数据移动后的距离，某个数据只有在比原先所在的簇更为接近现在所在的簇时，才能留在目前所在的簇，每次移动后簇的平均向量都重新计算，如此不断重复，直至一旦有任何移动，都会增加簇内的距离或减小簇间的非相似性为止。该法的局限性在于：①此方法在较大数据量时的扩展性和效率都较理想，但可能陷入局部最优。②只能用于类均数确定的情况下，若包含分类变量时就不适用。③必须提前确定类数。④受噪声和异常值的影响较大。目前常先使用凝聚算法确定类数和初始的类，再利用迭代重定位技术提高聚类的效果。k-medoids算法中用模式代替类均数，使用新的非相似性指标处理分类资料，用以频数为基础的方法对类的模式进行更替，而k-prototypes算法（k-means和k-medoids的结合）可以处理数值变量和分类变量的混合资料。EM（expectation maximization）算法是k-means算法的另一种扩展，把每一个体不是划为具体的某种类别，而是赋予其属于各类的概率。b. k-medoids法k-means算法对于异常值敏感，因为极端值可能歪曲资料的分布。k-medoids算法选择类的最中心的一点作为参照点，而不是类中所有个体的均数。当数据中存在噪声和奇异值多时，k-medoids算法比k-means算法具有更高的稳健性，因为一个类中具有代表性的中心点比该类中所有个体的均数更不易受异常值的影响。但k-medoids方法同样需要预先确定类数[10]。分割聚类分析适合于对具有相似性的基因进行分类。系统聚类和分割聚类是基因芯片数据分析中最传统、应用最广泛的方法，对于一般资料具有较理想的分类效果，但在处理复杂非线性及变量问的交互作用时效果较差。⑶主成分分析（principal component analysis，PCA）[28]在大规模基因表达数据的分析工作中，由于组织样本例数远远小于所观察基因个数（G&&N），如果直接采用前述聚类分析可能产生较大误差，故需要对聚类算法进行改进。目前已经提出很多改进的聚类方法，其中较为流行的方法是应用主成分分析方法对数据进行分析。主成分分析的目的是要对多变量数据矩阵进行最佳综合简化。使用的方法是寻找这些变量的线性组合&&称之为&主成分&（principal component），使这些主成分间不相关。为了能用尽量少的主成分个数去反映原始变量间提供的变异信息，要求各主成分的方差从大到小排列，第一主成分最能反映数据间的差异。主成分分析通过合并原来的维数得到更少的维数来表示对象，同时要求新的维数必须尽可能地反映原有维数所反映的信息，它有较少的信息丢失．主成分分析有助于简化分析和多维数据的可视化[17]。⑷自组织映射图网络（self-organizing map clustering，SOM）[10]所谓自组织特征映射是指神经网络中邻近的各个神经元通过侧向交互作用彼此相互竞争，自适应地发展成检测不同信号的特殊检测器。自组织映射网络图的基本原理是：将多为数据输入成几何学节点，相似的数据模式聚成节点，相隔较近的节点组成相邻的类，从而使多维的数据模式聚成2维节点的自组织映射图。SOM适合于复杂的多维数据的模式识别和特征分类等探索性分析，它允许对聚类的部分结构施加干预。相对于系统聚类中的严格结构和k-means聚类的无结构，SOM更灵活。与主成分分析（PCA）类似，SOM可以对数据集中的不同表达模式实现可视化，从而判断某种模式是否为另外一种模式的变异。SOM同样需要实现确定类数。⑸模糊聚类法（fuzz clustering）[10]在真实情况下，基因各功能类间的边界经常是不能截然分开的，模糊聚类适合于解决此类问题。该方法首先由Bezdek提出，后被Guthke用于基因芯片数据中的基因的分类。它给出向量（代表观察个体或基因）隶属于各类的隶属度，亦即该向量属于各类的概率。非监督模糊聚类的应用包括模糊c-means法、概率SOM和Gustafson-plaid法。⑹双向聚类（two-way clustering，TWC）基因表达谱常采用单向聚类法（one-way clustering），即要么以整个样本中特性相似的基因进聚类，或者以基因表达相似的样本进行聚类。对样本和基因同时进行聚类就是双向聚类法（two-wayclustering），目前基因表达谱的数据分析常用的双向聚类有基因剃须（gene shaving，GS）和格子模型（plaid models）。基因剃须是通过基因的共同表达值或表达量来鉴定基因的亚类，基因表达谱分析方法常用监督进行聚类，没有考虑一个基因可能属于多个类。基因剃须对基因或样本进行分类既可以是监督的，也可以是非监督的。基因剃须近年逐渐被应用于基因表达谱的分析中，Hastie[23]使用基因剃须方法分析了B细胞淋巴瘤患者的基因表达谱，鉴定了一小类可用于生存率预判的基因。2.6基因芯片数据分析的监督方法监督的方法又称判别分析（discriminant analysis）,以判别样本所属的类型。判别分析在已有数据的基础上建立分类器，并利用所建立的分类器对未知样品的功能或状态进行预测。与聚类分析不同，判别分析使用某种方法将研究对象分成若干类的前提下，建立判别函数，用以判定未知对象属于已知分类中的哪一类[16]。2.6.1&监督方法的数据降维监督的方法中的数据降维主要指数据选择，其目的有二：⑴挑选相对基本的数据了；⑵减少同济分析所必须的数据量。最简单的降维方法是，不断重复为每个数据加权的分类算法。首先用分类算法去除加权最小的数据，然后在剩下的数据中再用分类算法去除加权最小的数据，如此不断重复，直到这种处理已经失去统计显著性的时候，表明有重要信息已经被错误删去了。此时，立即停止计算，然后找回被误删的数据，这样剩下的数据则为非冗余的。当然，这个方法的缺点在于难以确定统计显著性的大小。2.6.2&各类监督方法技术简介⑴线性判别分析（1inear discriminant analysis，LDA）线性判别分析是指在输入变量上构造线性判别函数的方法。即寻找一种变换，使得在某种意义下类间分离性最大，类内相异性最小。它是一种有监督的维数约简方法[30][31]。线性判别分析的特点是计算简单，易于应用，一般具有较低的误差率，但不能处理基因（或个体）间的交互作用。因此，当基因（或个体）间存在复杂的交互作用时，线性判别分析不易发现数据中的规律性[10]。有研究指出，在基因芯片的分类中，Diagonal线性判别分析具有与最临近分类相接近的较高的判别性能，而Fisher线性判别分析的判别性能比其他方法要差[10]。另外，与LDA接近的还包括二次方判别分析等。Cho等[32]应用Fisher判别方法分析肿瘤患者的基因表达谱资料以判别肿瘤的分型；Dangond[33]等将Fisher判别方法应用于计算肌萎缩侧索硬化病的基因表达谱研究中。⑵k最临近分类法（k-nearest neighbor classfiers）k最临近分类法建立在通过类比进行学习的基础上，训练样本由n维计量变量描述，而每个观察个体由n维空间中的一点来描述。当给定一个未知样本，k-最临近分类算法将在模式空间中搜寻与此样本最临近的k个观察个体，这k个个体就是该位置观察个体的k个最临近点。一般采用欧氏距离来衡量临近程度。未知样本被赋予k个租赁金的个体中类数最多的类。与复杂的分类算法相比，k-最邻近算法具有简单、直观、误差率较低等特点，能够以&黑箱&的方式处理基因间的交互作用，但不能洞悉数据的结构。⑶决策树（decision trees）[16]决策树是一种常用于预测模型的算法，它通过将大量数据有目的的分类，从中找到一些有价值的，潜在的信息。它的主要优点是描述简单，分类速度快，特别适合大规模的数据处理[34]。⑷人工神经网络法（artificial neural network,ANN）ANN是一种应用类似于大脑神经突触联接的结构进行信息处理的数学模型。在ANN中，大量节点（&神经元&或&单元&）之间相互联接构成网络，即&神经网络&，以达到处理信息的目的。其优势是运行分析师无需在心中有特定模型，而且神经网络可以发现交互作用效果。Sawa等对酵母属基因表达谱数据进行欧式距离、相关系数、相互信息和基于神经网络的聚类分析，发现基于神经网络的聚类结果较前3种更为合理。⑸支持向量机（support vector machine,SVMs）[11][35]支持向量机是数据挖掘中的一个新方法。支持向量机能非常成功地处理回归问题（时间序列分析）和模式识别（分类问题、判别分析）等诸多问题，它通过训练一种&分类器&来辨识与已知的共调控基因表达类型相似的新基因。它起源于统计学习理论，研究如何构造学习机，实现模式分类问题。支持向量机使用结构风险最小化，使每一类数据之间的分类间隔最大。SVM用构建训练组的方法来学习如何区分不同的类别，它可以利用生物学的信息决定如何分组，也可以找出已分类组中的错误值。例如Williams [36]为了鉴定出肾母细胞瘤复发的基因表达谱模型，研究了27例肾母细胞瘤患者的肿瘤组织，其中13例2年内复发，对复发和未复发的肿瘤组织进行基因芯片实验，并应用支持向量机对基因表达谱数据进行分析，结果发现了一小类可能用于肿瘤预诊的基因。2.7数据的可视化方法[14]生物信息数据量大，形式复杂，直观地显示数据挖掘结果使其易于理解甚至关系到数据挖掘的成功与发展。目前已有一些可视化方法和工具，如Stanfold大学的TreeView软件采用色彩图与树图显示聚类的结果，树图能够清楚地显示层次聚类的每一步骤。另外，Bioconductor也提供芯片数据可视化的方法。&3. R & Bioconductor应用于表达谱芯片的数据处理3.1 R & Bioconductor简介R语言是一种计算机程序设计语言，也是一个开放式的软件开发平台，它具有强大的数学统计分析和科学数据可视化功能，能提供给各种数据处理和统计分析工具，如线性和非线性建模、经典的统计测试、时间序列分析、分类和聚类，同时也提供各种图形显示和分析工具。由于R语言是一个开放式的软件开放平台，软件开发人员可以再这个平台上不断扩充R语言的功能，并开发出面向特定应用的软件，如Bioconductor。Bioconductor实际上是一个开源和开放式的软件开发项目，该项目起始于2001年秋季，项目核心成员主要是哈佛医学院/哈佛公共卫生学院的Dana Farber癌症研究所生物统计组，还有来自美国和国际上的其他研究机构的一些研究人员。该项目的目标是建立多方面的、强有力的基因组数据的统计与图形分析方法，促进各种生物数据的集成，推动数据的综合分析和利用，促进各种生物数据的集成，推动数据的综合分析和利用，促进形成?%上一篇：下一篇：
技术部联系人：李小姐
移动电话：
邮政编号：510463
公司地址：广州萝岗科学城创新基地揽月路80号C区406-410
公司邮箱：
&联系人：石小姐
联系电话：020-&
邮政编号：510463
公司地址：广州萝岗科学城创新基地揽月路80号C区406-410
Copyright@ &
广州济恒医药科技有限公司版权所有&&&&&
电话：020-
地址：广州萝岗科学城创新基地揽月路80号C区406-410
邮编：510463
&&&&&&&&&&}