基因芯片（Affymetrix）分析2：芯片数据预处理

• 芯片质量评估
• 芯片数据前处理
• 识别关键基因变化
• 功能富集分析
• 基因表达模式识别

(本文于2013.09.04更新）

基因芯片技术的一个显著特点是通过寡聚核苷酸探针来检测基因。在上一节中使用ReadAffy函数读取CEL文件所得数据属于探针水平（probe level），即杂交信号。然而，在这一部分我们关注的重点则是将这些杂交信号转化为表达相关的度量值（亦即所谓的表达水平数据或expression level data）。具体而言，在探针水平的数据存储于AffyBatch类对象中，则属于ExpressionSet类对象所管理的范畴。为了满足这一需求，在R语言环境中可以调用多个用于处理基因芯片探针水平数据的软件包包括affy、affyPLM、affycomp及gcrma等工具包。这些工具包都经过精心设计以提高分析效率与准确性。如果尚未安装相关软件，请运行以下R语句进行安装：

    library(BiocInstaller)biocLite(c("affy","gcrma","affyPLM","affycomp"))

Affy芯片数据的预处理一般有三个步骤：

• 背景调整（background adjustment）
• 归一化过程（normalization process, 也可称为'标准化方法'）
• 总结（summarization process）

最后阶段完成对表达水平数据的获取工作。需要注意的是，在每一步骤中都涉及多种处理途径（算法）。采用何种处理手段往往因人而异，并会对其最终结果产生显著影响。由于不同级别的刊物或编辑可能会有不同的偏好，在实际应用中需要根据具体情况灵活选择适合自己的方式

1 需要了解的一点Affy芯片基础知识

Affy基因芯片采用25个碱基长的 probe设计。对于每个 mRNA样本而言,研究人员通常会使用 11 到 20 个独特的 probe 来进行检测,这些用于同一 mRNA 分析的一组 probe 被定义为其 "probe set"。鉴于所使用的 probe 长度较短,为了确保杂交反应的独特性, Affymetrix 公司针对每个基因设计了两种类型的 probe:其中一类被称为 perfect match (PM) 的 probe 阵列,另一类则被命名为 mismatch (MM) 的 probe 阵列。值得注意的是,在 MM 阵列中将 PM 阵列中的第 13 位碱基替换为其互补配对碱基,这样形成 PM 和 MM 探测器对

我们先使用前一节的方法载入数据并修改芯片名称：

    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)# 最好查看一下文件名称basename(cel.files)

    data.raw<-ReadAffy(filenames= cel.files)sampleNames(data.raw)<-paste("CHIP",1:length(cel.files),sep="")

用pm和mm函数可查看每个探针的检测情况：

    pm.data.raw<-pm(data.raw)head(pm.data.raw,2)

    mm.data.raw<-mm(data.raw)head(mm.data.raw,2)

上面显示的列名称就是探针的名称。而基因名称用probeset名称表示：

    head(geneNames(data.raw))

probeset不一定是和实际基因一一对应的，在处理过程中会对探针集进行基因名称映射时会看到相关部分

2 背景处理（background adjustment）

虽然自称是背景处理, 但这一步不仅涉及处理 background 值, 还涉及处理 noise signals. 需要注意的是 backgrounds 和 noise 是两个不同的概念, 例如在乡村夜晚时分的蛙鸣声尽管显得嘈杂但却十分稳定, 可以被视为 background, 如果突然来一声狗叫, 那么这就是 noise. 它们在统计上有明显的区别.

从理论上讲，在芯片设计中进行背景处理相对较为简单。Affy公司的MM设计初衷即在于检测非特异交叉信号。然而研究发现约有30%左右的MM探针所测得的结果反而比对应的PM探针更为显著（有趣的是……）。由此可见，在那些看似合理的方法背后或许隐藏着一些值得深入探讨的问题

该软件包提供了芯片背景噪声消减的相关函数，并且其中最为常用的是bg.correct()这一功能。此外，在分析RNA微阵列数据时，默认采用的依然是最为经典的MAS与RMA分析方法。

MAS方法采用网格划分策略将芯片划分为k个网格区域，默认值k设为16；通过每个网格区域使用信号强度最低的2%探针来计算背景值与噪声。RMA方法的工作原理较为复杂，请参考相关文献。

Ricardo A. Irizarry et Brazelton Hobbs et Felix Collin et al., were conducted to explore the techniques for normalization and summarization of data derived from high-density oligonucleotide array probes at the level of individual probes during the process of microarray analysis and interpretation in biostatistics research., Biostatistics (Oxford), Volume: 4 (Issue): pages 249–64 (Year): October (Month) 2003 (Year). Numbers: [11;18;27;232;241;432;443].

    data.rma<-bg.correct(data.raw,method="rma")data.mas<-bg.correct(data.raw,method="mas")class(data.rma)

    class(data.mas)

    class(data.raw)

通过该函数ReadAffy()解析得到的CEL芯片数据采用AffyBatch类数据形式存储；经背景校正后的数据依旧保持在AffyBatch类数据形式中。

MAS方法施加于实验样本后导致PM与MM指标数值均被重新评估。在实施RMA方法时主要基于PM探针数据进行分析，在经过背景校正处理之后, MM探针的数据保持不变

    head(pm(data.raw)-pm(data.mas),2)

    head(pm(data.raw)-pm(data.rma),2)

    head(mm(data.raw)-mm(data.mas),2)

    #差值为全部为0，说明rma方法没有对mm数值进行处理head(mm(data.raw)-mm(data.rma),2)

    identical(mm(data.raw),mm(data.rma))

3 归一化处理（normalization）

同一份RNA样品通过多块类型一致的探针进行杂交实验所得的结果（如信号强度等指标）可能会有明显的差异甚至显著差别为此必须采用统一化的数据处理方法以确保来自不同探针组的数据具有可比性这些统一化方法种类繁多其中由Affymetrix公司提供的normalize函数是一个常用工具该函数采用Affybatch对象作为输入并通过指定处理方法实现数据标准化

3.1 线性缩放方法

该方法在Affymetrix公司的软件版本4.0及5.0中被采用。该方法首先选择一个芯片作为基准。随后将其他所有芯片与基准芯片进行线性回归分析。利用无截距的一元线性回归模型对各芯片的信号值进行标准化处理。Affymetrix公司的软件在回归分析前剔除了强度最高的2%和最低的2%的数据点。使用相对简便的方法，在完成背景校正后即可执行归一化分析，请问您是否知道具体的R代码？

    data.mas.ln<-normalize(data.mas,method="constant")head(pm(data.mas.ln)/pm(data.mas),5)

    head(mm(data.mas.ln)/mm(data.mas),5)

可以看出，在应用线性缩放方法时，默认将第一块芯片作为基准，并对其数值未经过缩放处理的状态进行记录；而对于其余所有芯片，则进行了相应的放大处理。对于同一块芯片而言，在各个探针区域中其放大系数是相同的；此外，在PM与MM区域中所采用的放大策略完全一致

3.2 非线性缩放方法

此方法获得的结果显著优于线性方法，在非线性拟合过程中，并未取整个芯片而是选取其中一列作为基线

    data.mas.nl<-normalize(data.mas,method="invariantset")head(pm(data.mas.nl)/pm(data.mas),5)

可以看到，同一芯片不同探针的信号值的缩放倍数是不一样的。

3.3 分位数（quantile）方法

该方法假定所有单片探针信号的经验分布函数具有相同的形态，并且将任意两张芯片的数据绘制到QQ图上时应当呈现一条斜率为1且截距为0的直线

    data.mas.qt<-normalize(data.mas,method="quantiles")head(pm(data.mas.qt)/pm(data.mas),5)

3.4 其他

如循环局部加权回归法（Cyclic loess）和 Contrasts方法。

    data.rma.lo<-normalize(data.rma,method="loess")

    data.rma.ct<-normalize(data.rma,method="contrasts")

4 汇总（Summarization）：

最后一步汇总时采用合适的统计方法，并通过probeset（包含多个探针）进行杂交信号的整合来计算得到相应的表达量。affy包中提供的函数是computeExprSet()。需要注意的是，在调用computeExprSet()函数时除了需要指定统计方法之外还需要设定PM校正的方式：具体来说就是调用computeExprSet(x, pmcorrect.method, summary.method, ...)的形式进行操作

两个参数可以设定的值可以通过下面函数获得：

    pmcorrect.methods()

    generateExprSet.methods()

常见的汇总方式包括medianpolish、liwong以及mas。liwong的方法仅采用pm作为背景校正（即pmcorrect.method="pmonly"）。例如：

    eset.rma.liwong<-computeExprSet(data.rma.lo,pmcorrect.method="pmonly",summary.method="liwong")

计算过程耗时若干小时，并请静静等待最终完成。最后的结果为 ExpressionSet 类型数据：

    class(eset.rma.liwong)

    class(data.raw)

5 三合一处理和“自动化”函数：

掌握芯片预处理的基本流程后无需额外学习你就可以利用一个R函数来进行三步操作

    # NOT RUN. 函数说明，不可运行。expresso(afbatch,bg.correct=TRUE,bgcorrect.method=NULL,bgcorrect.param=list(),normalize=TRUE,normalize.method=NULL,normalize.param=list(),pmcorrect.method=NULL,pmcorrect.param=list(),summary.method=NULL,summary.param=list(),summary.subset=NULL,verbose=TRUE,widget=FALSE)

例如：

    # NOT RUN:eset.mas<-expresso(data.raw,bgcorrect.method="mas",normalize.method="constant",pmcorrect.method="mas",summary.method="mas")

可调用bgcorrect.methods、pmcorrect.methods以及express.summary.stat.methods函数进行多种分析方法的调用

expresso运行速度较慢，在现有的R软件包中有一些提供了更快捷的预编译三合一功能模块，并且例如affyPLM软件包中的three-in-one function就可以实现这一功能。通过使用affyPLM软件包提供的各种处理方案可以显著提升分析效率与准确性；若感兴趣可自行深入学习其详细实现机制

本文将简述三种自动化的功能模块。这些函数已预先设置了各项处理环节所需的技术手段及参数设置，并且无需自行配置。

5.1 mas5函数

由affy包提供：

    eset.mas5<-mas5(data.raw)

mas5据说是由expresso函数与mas方法打包而成的集成模块。然而，在运行以下代码所得的输出结果似乎与 mas5(data.raw)的行为存在差异（建议自行测试以确认）：

    # NOT RUN:expresso(data.raw,bgcorrect.method="mas",normalize=FALSE,pmcorrect.method="mas",summary.method="mas")

5.2 rma函数

也是由affy包提供的同样采用了rma法作为背景处理；归一化处理采用了分位数法；而汇总采用了medianpolish：

    eset.rma<-rma(data.raw)

它等价于：

    # NOT RUN：eset.rma<-expresso(data.raw,bgcorrect.method="rma",normalize.method="quantiles",pmcorrect.method="pmonly",summary.method="medianpolish")

但rma函数是预编译好的C语言函数，速度非常快！

5.3 gcrma函数

基于gcrma包实现的Affy软件中包含一个名为mas的方法用于进行微芯片分析

    library(gcrma)eset.gcrma<-gcrma(data.raw)

6 芯片处理方法的优劣评估

面对如此多的芯片预处理方法，选择最优方案确实让人感到困惑。究竟该如何决定呢？知道得越多反而让人难以抉择。幸好这些已经有人做了，在此领域已经有了深入研究，并且研究者Rafael Irizarry和Leslie Cope开发了一个名为affycomp的R软件包来系统性地评估各种预处理方法。

评估方法的有效性必须基于数据, 且这些数据应包含可追溯的因素. affycomp要求两种类型的数据显示: 一种是通过DNA分子杂交技术获得的稀释数据系列,称为稀释数据;另一种则是利用内参或外标RNA与基因序列结合的技术得到的结果,称为内参/外标数据系列. 这里的Spike-in RNA指的是在目标物种中没有存在但能在该芯片上被检测到的特定RNA分子. 例如, Affymetrix公司的拟南芥芯片上就标记了几个人类或细菌基因的独特探针. 因为所有样本都经过预先设定好的稀释倍数处理, 所以在计算过程中这些参数都可以被精确量化. 因此, 在计算过程中这些参数都可以被精确量化. 为了确保实验结果的高度准确性, 这些严格的数据验证步骤都是不可或缺的. 然而, 大多数研究者出于成本考虑而选择跳过这一环节: 只进行几张芯片上的常规分析工作即可采用已被广泛认可的方法如RMA或MAS

7 Session Info

    sessionInfo()

Author: ZGUANG@LZU

Created: 2013-09-04 三 14:42

Emacs 24.3.1 (Org mode 8.0.5)

Validate XHTML 1.0

转载于:https://www.cnblogs.com/huzs/p/3741997.html