Bioconductor分析基因芯片数据

使读者初步了解使用Bionconductor完成基因芯片预处理的流程

接着详细讲解戏弄i按预处理和数据分析等内容

最后深入了解实际工作中会遇到的芯片处理问题以及如何用学到的只是解决问题

目的：掌握芯片分析的整体框架，自行学习其他厂商或种类（例如SNP芯片或CHIP-chip芯片）的芯片处理方法

5.1快速入门

例5-1 从数据包CLL中载入芯片数据，完成预处理，最后获得基因（探针组）表达矩阵。注意，探针组表达矩阵的行对应的行对应的探针组，而不是基因，基因和探针组的关系见5.2.1.这段程序从载入原始数据（CEL文件）开始，通过预处理得到基因表达矩阵，是芯片数据处理的一个必须步骤

CLL 数据是慢性淋巴白血病（Chronic lymphocytic leukemia,CLL)数据集，采用了Affymetrix公司的HG_U95Av表达噗芯片，测量了24个样本，12625个探针

例5-1采用的实验设计方式：两组之间是对照试验（control test)，每组内都是平行实验（parallel test),对照实验，简单来说就是为了阐明某种单一因素的效应或者影响，在保持其他因素不变的前提下，测试一定数据的实验组样本呢和对照组样本，并对结果进行比较。平行实验，简单来说就是对同样的一组样本取两个以上相同的样品，以完全一致的条件下进行试验，测试结果的稳定性

5.2基因芯片基础知识

5.2.1探针组

一张基因芯片（以affymetrix表达谱芯片为例）可以包含上万个的探针（通常由25个碱基组成），他们整齐有序地印刷在芯片上。一组探针或者探针组（probe set),来自于一个基因，通常由20对或者11对探针组成，每一对探针都由匹配探针（perfect match，PM）和错配探针（Mi是match，MM)组成，成为探针对（probe pair),MM与PM的序列只有正中央的那个碱基不同，其余的都一致。但是，在一些高密度芯片中，例如外显子芯片（Exon array),每个探针组只有4个PM探针，没有MM探针。

探针序列的来源叫做参考序列，通常来自于公开的核酸数据库（例如NCBIGeneBank或RefSeq)对于不同的芯片类型，探针组在参考序列中的分布不同，3’表达谱芯片的探针组排布在参考序列3‘末端附近的一至两外显子上，外显子芯片中，每个长度大于25个碱基的外显子都有针对他的探针组：铺瓦芯片（Tilling array)中，探针组覆盖了几乎所有的外显子和内含子

需要强调的是，芯片数据领域提到的基因表达矩阵往往是以探针组而忽视以基因为单位的，即每行都对应一个探针组的表达量。后面将要降到的差异基因分析也是找打显著性差异的表达的探针组，然后通过ID映射才对应到探针组代表的基因，探针组与基于的关系往往是多个探针组对应一个基因。但是在实际应用中，经常不太注意区分，探针组有时也会被叫做基因

5.2.2主要的芯片文件格式

主要的是CEL文件

affymetrixi芯片原始数据最常用格式为CEL文件，也是芯片预处理和分析的出发点。CEL文件的主要内容就是每个“cell"的灰度信息，"cell"是整个芯片图像划分后得到的小网格，每个小网格中的图像被看作来自一个探针，自CellHeader开始，每行数据对应芯片上的一个”cell"位点，包含5列信息，依次为X坐标,Y坐标，灰度的平均值，灰度的标准差以及用了多少个像素来求这个平均值

CEL文件只提供了每个探针的灰度信息，还需要基因芯片探针排布的信息（即哪个探针来自哪个探针组），才可以得到芯片上每个探针组对于的表达数据，这就需要CDF文件。另一个重要的是probe文件，他提供了探针的序列信息。afftymetrix公司为每种型号的芯片都提供了对应的CDF文件和Probe文件。CDF文件中的对应关系用户可以自行更改，例如为了应多多个探针组 的ID对应到同一基因ID的现象，有些研究季候就把对应到同一个基因的多个探针合并为一个探针组，并提供修改后的CDF和Probe文件

图5-3B是affymerixHG-U133A芯片的Probe文件的部分内容，他只包括了一个探针组（名称是“200688_at")的所有探针，共11条序列，文件中第2和3列是对应探针所在的X和Y 坐标，第4列是序列的第13个奸计（中心）位置对齐到一致性序列的相对位置，第5列是对应探针的序列，最后是样品与探针杂交的方向

5.3基因芯片数据预处理

基因芯片数据预处理的目的是将探针水平的数据（杂交信号）转换成基因表达数据，主要的数据结构有Affybatch类和Expressionset类，前者用于存储探针水平的数据（相当于CEL文件的内容），而后者用于存储表达水平的数据（相当于基因表达矩阵的内容）。预处理通过质量控制，剔除不合格的芯片（数据），只保留合格的进入下一步处理。然后通过标准化，将所有芯片数据中的基因表达只变换到一个可以比较的水平，用于后续分析

5.3.1数据输入

例5-1中，芯片数据的输入是从数据包中得到的，但是在实际应用中，更常见的情况是从CEL中获得数据，无论是数据包还是文件输入，读入的数据会存入一个“affybatch"类型的对象中，可以通过执行help（affybatch)获得更详细的介绍

头文件：用于描述实验样本、平台等相关信息，其中包括phenoData,featureData,protocolData以及annotation等几个类

assayData:这是affybatch类必不可少的，他的第一个元素是矩阵类型，用于保存基因表达矩阵。该矩阵的行对应不同的探针组（probe sets)，用一个无重复的索引值表示，列对应不同的样品。当使用exprs方法时，调取的就是这个基因表达矩阵

experimentData:一个MIAME类型的数据，设计这个MIAM类的目的就是用于保存MIAME原则建议的注释信息.MIAME原则是一组指导方针，他建议了一组标准来记录与基因芯片实验设计相关的资料

5.3.2质量控制

质量控制对于后续的分析至关重要，原始图像（DAT文件）级别的质量控制一般用个芯片公司自带的软件（如affymetrix公司的GCOS)完成。本节中，质量控制主要集中在CEL文件级别的处理，从最简单的直观观察，到平均值方法，再到比较高级的数据拟合方法。这三个层次的质量控制分别功能分别用image函数simpleaffy包和affyPLM包实现

直观的查看一下芯片上所有位点的灰度图像

image函数表示选取的CLLbatch中的第一个基因芯片（即“CLL10.CEL"）的数据，然后调用image函数根据CEL文件中的灰度信息画图，affymetrix芯片在印刷时会在四个角印刷特俗的花纹，并且在左上角印刷芯片的名称，花纹与芯片名称可以帮助我们借助这个图像分辨率来了解芯片数据是否可靠。如果无法分辨四角花纹或芯片名称，很可能数据有问题

根据image函数的图像信息，可以对芯片的信号强度产生一个总体认识：如果图像特别黑，说明信号强度低；如果图像特别亮，说明信号强度很可能过饱和

尺度因子affymetrix公司规定，用于比较的芯片之间的尺度因子的比例必须小于三

检测值（detection call)和检出率（percent present）:一组探针能否被检测到，用检测值有（present，简称R)、无（Absent,简称A)和不确定（Marginal presen,简称M)来表示检测范围的上下边界（a1及a2)选用了默认值0.04和0.06.检出率，是用所有检测值为p的探针数量除以芯片所有探针组数控得出的百分比。如果检出率过低，表示大部分的基因都未被检测到，很难说明是该芯片实验有问题，还是这个样品的大多数基因本身就很难检测到，有原因是表达量基地或是其他。因此，需要看多个样品之间的相对差别，如果有的样品的检出率与其他的有比较大的差别，那很可能该样品出现了问题

平均背景噪声（average background）：对于每一块芯片，根据所有的MM值作出统计，可以得到背景噪声的平均值、最小值和最大值。往往较高的背景噪声都伴随着最低的检出率，因此这两个指标可以结合使用

标准内参（internal control genes）:mRNA是按照5‘端到3’端的顺序来降解的，芯片探针组也是根据这个顺序来设计的，因此探针组的测量结果可以体现这一趋势。因为大部分的细胞都有β-action和GAPDH基因，所有affymetrix在大部分的芯片里都将他们设置为一组观察RNA降解成都的内参基因。根据这两个基于设计的探针组很好的涵盖了他们3‘端和5’端的每一个区段。通过比较他们3‘端相对于中间或者5’端的信号强度，可以很好地指示出实验质量。affymetrix建议这个比值对于β-action不大于3，对于GAPDH不大于1.25，即可以说明这个芯片的质量可以接受。如果这个比值很高，表明不完整的β-action或者GAPDH的存在，可能源于体外转录不好或者降解非常严重。如果使用的是affymetrix的小样本实验流程（small sample protocol）而不是常用的标准流程（standard protocol），建议使用3’端相对于中间的比值。原因是小样本流程有更扩增次数，有可能产生更多较短的转录序列，不可避免的带来3‘端的偏倚。为了验证杂交的质量

根据上述标准，可以使用Bionconductor的simpleaffy包对affymetrix芯片数据进行质量评估，最后得到质量控制总览图（图5-8）

qc图的看法，图5-8是CLL数据集中全部24个芯片数据的质量控制总览图。图5-8中从左至右，第一列是所有样品的名称；第2列是两个数字（对应每个样品），上面是以百分比形式出现的检出率，下面的数字表明平均背景噪音；第3列（"QQ stats")最下面的横轴是尺度因子等指标对应的坐标，取值范围从-3到3，用浅蓝色虚线作为边界。第3列用到了三项指标：尺度因子、GAPDH3'/5'比值和action3比值（记做graph3/graph5和action3/action5),分别用实心圆、空心圆和三角标志表示出来。另外，如果第三列中出现红色的”bioB"字样，说明该样品中未能检测到BioB

简单地讲，所有指标出现蓝色表示正常，红色表示可能存在质量问题。但是根据实际情况不同，还要进一步分析。一般来讲，如果有一个芯片各项指标都不太正常，尤其是BioB无法检测到，建议判定为该芯片实验失败。如图5-3中的样品”CLL15.CEL",这个数据的检出率（38.89%）明显低于其他样品，action3/action5远大于3，而且没有检测到BioB,因此可以判定此数据无效。如果多个芯片都出现了相同的问题，原因则可能是多方面的；如左侧第2列24个芯片的检出率和背景噪声都很高，原因是阈值设定过高，如果降低阈值，大部分就会变蓝；再如，全部芯片都不能检测到BioB,有可能是嵌入探针所针对的DNA溶液加入比例不对

基于平均值家建设的评价指标都有一个,默认的假设，那就是对于每一块新片，质量是均匀的，不会随着位置变化发生较大的变化。但如果关注芯片的每个小格（Grid),就会发现格与格之间的质量也是有差异的，这可能由于芯片印刷的问题，也可能是杂交过程中出现的问题。那么如何才能得到比较可靠的质量评估，这需要设计多种能反映芯片数据全貌的指标综合分析从而得出最终的结论。这些指标要在对原始数据拟合（回归）的基础上计算得到，然后以图的形式显示，包括：权重（weights)&（residuals)图、相对对数表达（relative log expression，RLE)箱线图、相对标准差（normalized unscaled standard errors,NUSE)箱线图、RNA降解曲线、聚类分析（cluster analysis）图、主成分分析（principal component analysis，PCA)图、信号强度分布图及MA图等，以上功能由Bionconductor中的affyPLM包实现

一般情况下，在权重图中，绿色代表较低的权重（接近0），白色、灰色代表较高的权重（接近1）；在残差图中，红色代表正的高残差，蓝色代表负残差；在残差符号中，红色代表正的残差，蓝色代表负的残差。如果权重和残差都是随机分布的，应该看到绿色均匀分布的权重图和红蓝均匀分布的残差图。图5-9中，左上为原始图像，右上为权重图，左下为残差图，右下为残差符号图。另外，还可以看到，图中左上部出现了一些白色的条块，这是正常的现象，因为有些时候，探针会按照GC比率（GC ratio）排布从而导致白斑的，那什么样的权重和残差图是不可接受的呢

在对比实验中，即使是相互比较的对照组与实验组之间，大部分基因的表达量还是应该保持一致的，平行实验之间一致性更强。相对对数表达（RLE)箱线图可以反映上述趋势，它定义为一个探针组在某个样品的表达值除以该探针组在所有样品中表达值的中位数后取对数。一个样品的所有探针组的RLE的分布可以用一个统计学中常用的箱型图形表示。如果使用RLE箱线图来控制CLL数据集的实验质量，每个样品的中心应该非常接近纵坐标0的位置（图5-11）。如果个别样品的表现与其他样品的表现与其他大多数明显不同，那说明这样品有问题

NUSE是一种比RLE更为敏感的质量检测手段。如果根RLE箱线图对某个芯片的质量产生怀疑，那么再结合NUSE图，这种怀疑就可以确定下来。NUSE定义为一个探针组在某个样品的PM值的标准差除以该探针组在各样品中PM值标准差的中位数，如果所有芯片的质量就是非常可靠的话，那么他们的标准差会十分接近，因此他们的NUSE值会都在1附近。然而，如果有某些芯片质量有问题的话，就会严重地偏离1，进而导致 其他芯片的NUSE值偏向相反的方向。当然，还有一中非常极端的情况，那就是大部分芯片有质量问题，但是他们的标准差却比较接近，反而会显得没有质量问题的NUSE值明显偏离1，所以必须结合RLE及NUSE两个图才能作出更可靠的判断。例如结合图5-11和6-12，可以看出CLL1和CLL10的质量明显有其他yan