CN103745137A

CN103745137A - 一种跨芯片平台的基因表达数据整合方法

Info

Publication number: CN103745137A
Application number: CN201410044121.6A
Authority: CN
Inventors: 杭兴宜; 陈�胜
Original assignee: Faith Bo Auk Biological Information Technology (beijing) Co Ltd
Current assignee: Aiji Taikang (Jiaxing) Biotechnology Co., Ltd.
Priority date: 2014-01-30
Filing date: 2014-01-30
Publication date: 2014-04-23
Anticipated expiration: 2034-01-30
Also published as: CN103745137B

Abstract

本发明属于生物信息技术领域。本发明提供了一种跨多个不同芯片平台整合基因表达数据的方法，包括如下步骤：多个芯片平台的基因表达谱标准化预处理；合并不同芯片平台中共同的基因表达数据；根据基因间在多个芯片平台的表达相似性将基因分为k个子集；最小二乘法计算每个基因子集内不同芯片平台的表达线性关系

；利用公式

将所述不同芯片平台的基因表达值标准化到同一变化范围，得到标准化的基因表达矩阵，其中符号含义如说明书中所定义。

Description

一种跨芯片平台的基因表达数据整合方法

技术领域

本发明属于生物信息技术领域，特别是基因表达数据分析领域。

背景技术

目前，微阵列芯片已经发展成为系统研究生物学问题的常用高通量实验技术，并且有不同类型的芯片平台及其制造厂商。多年来，已经积累了大量的芯片数据集如美国国立生物技术信息中心NCBI的GEO芯片数据库和欧洲生物信息研究所EBI的ArrayExpress芯片数据库。其中，NCBI GEO芯片数据已收集了大约1,008,760个样本共12,090个实验的数据，EBI ArrayExpress芯片数据已收集了43,124个实验共计1,223,250个微阵列芯片数据。由于芯片实验的价格相对昂贵，样本收集耗费工作量等客观原因，大多数科研项目仅会选择相对少量的样本进行实验甚至仅进行一个样本的检测，因此未能严格满足科学实验设计的要求，使得结果的可信度不足。所以，怎样合理充分的利用现有的数据资源，筛选与自己研究目相同相近的芯片数据进行整合分析研究成为一个亟需解决的问题。对于来自同一芯片平台的数据可以较直接地进行数据整合，但对于不同类型芯片间的数据整合目前还没有一个行之有效的方法。

已有基因表达数据整合的研究，其算法的步骤如下：1) 芯片数据预处理进行log₂转换；2) 根据不同芯片平台间信号的最大值，确定两款芯片间的基因表达转化的比值；3) 整合两个芯片间共同存在的基因；4) 消除实验批次影响。上述算法中存在以下缺点：1) 芯片预处理方法对后续的分析至关重要，仅仅log2转化而忽略了采用了何种背景校正方法；2) 评估基因在两个平台间表达的线性关系时，仅仅考虑了芯片的最大表达值，而未考虑基因在不同区段间的线性比例不同；3) 没有对数据进行样本间中值和方差标准化，不利于数据间的比较。

鉴于上述原因，本领域中仍然需要克服上述一项或多项缺点的跨芯片平台的基因表达数据整合方法。

发明内容

在第一方面中，本发明提供一种跨芯片平台的基因表达数据整合方法。包括以下步骤：

步骤1：获得要整合的两个或多个基因芯片平台表达数据，对于每个基因芯片平台而言，表达数据为一个矩阵，矩阵的每一行表示一个基因，矩阵的每一列表示一个被检测的样本。

步骤2：对两个或多个基因芯片平台中共同存在基因，生成合并的基因表达矩阵，矩阵的每一行表示一个共同存在基因，矩阵的列由所述不同基因芯片平台检测的样本构成，每一列表示一个被检测的样本。

步骤3：基于基因在所述合并的基因表达矩阵中的表达模式相似性将基因分割为多个基因子集。

步骤4：对于步骤3得到基因子集，利用基因子集内基因在所述每个不同芯片平台的表达值，得出不同芯片平台的表达线性关系，利用所述关系将所述不同芯片平台的基因表达值标准化到同一变化范围，得到标准化后的基因表达矩阵。

在一个具体的实施方案中，本发明是这样实现的，主要包括如下流程：

步骤2：对两个或多个基因芯片平台中共同存在的基因，生成合并的基因表达矩阵，矩阵的每一行表示一个共同存在基因，矩阵的列由所述不同基因芯片平台检测的样本构成，每一列表示一个被检测的样本。

步骤3：基于基因在所述合并的基因表达矩阵中的表达模式相似性将基因分割为m个基因子集。

步骤4：对于步骤3得到每个基因子集，计算基因子集内基因在所述每个不同芯片平台的表达值均值，利用下式得出不同芯片平台在该基因子集内的表达线性关系

，其中s表示m个基因子集中的一个子集，exp_s1是由基因子集s内基因在芯片平台1中的表达值均值构成的数组，exp_s2指基因子集s内基因在芯片平台2中的表达值均值构成的数组，采用最小二乘法估计

，

，其中h为基因子集s包含的基因数；利用公式

将所述不同芯片平台的基因表达值标准化到同一变化范围，得到标准化的基因表达矩阵。

其中X为一个n行m列由0和1构成的矩阵，且如果基因n属于基因子集m，则X_nm=1；基因n不属于基因子集m，则X_nm=0，n为两个芯片平台共同存在基因的数目，m为基因子集的数目，A和B为m行k列的矩阵，A和B矩阵的行元素分别为对应基因子集中计算的a_s和b_s，“

”表示两个矩阵相乘，“

”表示两个矩阵的点乘(两个矩阵相同位置的元素相乘) 。

A和B矩阵表示如下：

在一个更具体的实施方案中，本发明是这样实现的，主要包括如下流程：

步骤1：获得要整合的两个或多个基因芯片平台表达数据，所述数据为基因芯片的探针荧光值，代表各个基因表达，优选所述基因芯片的探针荧光值数据已经经过背景校正和标准化，所述背景校正和标准化优选采用RMA方法进行，将每个芯片基因平台的探针荧光值表示对应基因的表达值（例如，具体方法如下：如果一个基因仅对应一个探针，则该基因的表达值为其对应探针的荧光值；如果一个基因对应多个探针，则此基因的表达值采用其对应的多个探针荧光值的中值或均值表示），即基因转录本多少的定量值，全部基因(Symbol 1-n)在全部样本(sample 1-k)中的表达值(exp 1-nk)构成即基因表达矩阵(exp)，矩阵的每一行表示一个基因，矩阵的每一列表示一个检测样本，基因可采用任何标识形式表示，例如以“Official Symbol”或“Entrez Gene ID”表示，表达矩阵如下。

步骤3：基于基因在所述合并的基因表达矩阵中的表达模式相似性将基因分割为m个基因子集（常用方法有层次聚类、k均值聚类和SOM自组织映射聚类，优选m>=10，同时每个基因子集包含的基因数目不少于30）。

步骤4：对于步骤3得到每个基因子集，计算基因子集内基因在所述每个不同芯片平台的表达值均值，利用下式得出不同芯片平台在该基因子集内的表达线性关系，其中s表示m个基因子集中的一个子集，exp_s1是由基因子集s内基因在芯片平台1中的表达值均值构成的数组，exp_s2指基因子集s内基因在芯片平台2中的表达值均值构成的数组，采用最小二乘法估计

，

，其中h为基因子集s包含的基因数；利用公式

将所述不同芯片平台的基因表达值标准化到同一变化范围，得到标准化的基因表达矩阵，其中X为一个n行m列由0和1构成的矩阵，且如果基因n属于基因子集m，则X_nm=1；基因n不属于基因子集m，则X_nm=0，n为两个芯片平台共同存在基因的数目，m为基因子集的数目，A和B为m行k列的矩阵，A和B矩阵的行元素分别为对应基因子集中计算的a_s和b_s，“

”表示两个矩阵相乘，“

”表示两个矩阵的点乘(两个矩阵相同位置的元素相乘)。

A和B矩阵表示如下：

在优选的实施方案中，本发明的跨基因芯片平台的基因表达数据整合方法还包括：

步骤5：对基因表达矩阵进行全局中值标准化和样本间方差标准化，其中所述全局中值标准化通过

进行，m_j指基因表达谱矩阵第j列的中值，m为整个矩阵的中值；所述样本间方差标准化通过

进行，

指基因表达谱矩阵第j列的均值，

值矩阵第j列数据的标准差。

步骤6：基于非参数经验贝叶斯模型（例如利用Combat软件）对基因表达谱进行实验批次影响校正（W. Evan Johnson and Cheng Li 2006 Biostatistics），得到不同芯片平台间校正后的基因表达谱数据。

在第二方面中，本发明提供了一种利用两个或多个基因芯片平台获得一套基因表达数据的方法，所述方法包括如下步骤：

步骤1：利用所述两个或多个基因芯片平台分别对样品进行实验，获得各自的基因表达数据，所述数据为基因芯片的探针荧光值，代表各个基因表达，优选所述基因芯片的探针荧光值数据已经经过背景校正和标准化，所述背景校正和标准化优选采用RMA方法进行，将每个芯片基因平台的探针荧光值表示对应基因的表达值（例如，具体方法如下：如果一个基因仅对应一个探针，则该基因的表达值为其对应探针的荧光值；如果一个基因对应多个探针，则此基因的表达值采用其对应的多个探针荧光值的中值或均值表示），即基因转录本多少的定量值，全部基因（Symbol 1-n）在全部样本（sample 1-k）中的表达值（exp1 1-nk）构成即基因表达矩阵（exp），矩阵的每一行表示一个基因，矩阵的每一列表示一个检测样本，基因可采用任何标识形式表示，例如以“Official Symbol”或“Enrez Gene ID”表示，表达矩阵如下：

步骤2：对两个或多个基因芯片平台中共同存在基因，生成合并的基因表达矩阵，矩阵的每一行表示一个共同存在基因，矩阵的列由所述不同基因芯片平台检测的样本构成，每一列表示一个被检测的样本；

步骤3：基于基因在所述合并的基因表达矩阵中的表达模式相似性将基因分割为m个基因子集（常用方法有层次聚类、k均值聚类和SOM自组织映射聚类，优选m>=10，同时每个基因子集包含的基因数目不少于30）；

，

，其中h为基因子集s包含的基因数；利用公式

”表示两个矩阵相乘，“

”表示两个矩阵的点乘(两个矩阵相同位置的元素相乘)。

A和B矩阵表示如下：

在优选的实施方案中，本发明的利用两个或多个基因芯片平台获得一套基因表达数据的方法还包括：

进行，指基因表达谱矩阵第j列的均值，

值矩阵第j列数据的标准差；和/或

在本发明第一、第二方面的方法中，所述两个或多个基因芯片平台选自Affymetrix基因芯片平台和Agilent基因芯片平台。

在本发明的方法中，所述Affymetrix基因芯片平台为Affymetrix HG U133A，所述Agilent基因芯片平台为Agilent G4112F

本发明的优点在于：

1）不同的芯片平台均采用RMA (Robust Multi-Array Average expression measure)标准化方法；

2）根据基因在多个芯片平台的表达模式聚类结果，将基因分割成不同的子集，评估每个基因子集中两个平台间的表达线性关系；

3）对于合并后两个平台的基因表达数据，我们的方法不仅考虑全局样本的表达中值，同时还对基因表达方差变化进行了标准化。一般情况下同研究组内的样本基因表达具有一定相似性，组内样本基因表达会服从相同的分布，因此进行中值和方差标准化会使整合后的基因表达数据更合理，更具可比性；

4）我们利用非参数经验贝叶斯模型算法Combat消除基因表达数据中的实验批次影响，Combat相比其他的实验批次校正方法具有更好的表现。

本发明的创新点在于：对于不同芯片平台间共同存在的基因集，依据基因在多个芯片平台的表达模式分割为多个子集，评估每个基因子集中两个芯片平台的线性关系，最后依据基因所属基因子集信息，选用相应的a_k和b_k进行平台间数据标准化。

附图说明

图1示出本发明所述跨芯片平台基因表达数据整合算法的流程图；

图2示出HG U133A和Agilent G4112F某一基因子集的表达数据散点图；

图3示出线性校正后HG U133A和Agilent G4112F芯片基因表达箱式图；

图4示出中值标准化后HG U133A和Agilent G4112F芯片基因表达箱式图；

图5示出方差标转化后HG U133A和Agilent G4112F芯片基因表达箱式图；

图6示出HG U133A和Agilent G4112F芯片跨平台数据整合最终结果；

图7示出HG U133A和Agilent G4112F芯片跨平台数据整合前后基因表达总体相关性比较图；

图8示出GSE14520和GSE46480跨芯片平台数据整合最终结果。

具体实施方式

RMA：Robust Multi-Array Average expression measure，是一种较为鲁棒的芯片数据预处理方法。

Official Symbol：由基因命名委员会提供的一种通用的基因名称标识。

Enrez Gene ID：由美国国家生物技术信息中心提供的一种基因编号。

Combat：W.Evan Johnson和Cheng Li 2006发表的一个消除实验批次影响的软件。

实施例一、为了评估不同芯片平台间基因表达的差异，实现基因表达数据的跨平台整合，分别以Affymetrix HG U133A和Agilent G4112F芯片平台各检测了10例乳腺癌组织样本和10例正常癌旁组织样本的基因表达。

（一）、Affymetrix HG U133A芯片实验过程如下：

第1步：RNA的抽提。

按照制造商提供的说明书，使用QIAGEN’s RNeasy Total RNA Isolation kit抽提人类乳腺癌组织和癌旁组织的总RNA；使用QIAGEN’s Oligotex Direct mRNA kit从总RNA中抽提mRNA。

第2步：RNA沉淀。

用QIAGEN’s RNeasy Total RNA Isolation kit分离或洗涤后没有必要沉淀总RNA。调整洗脱体积以制备cDNA合成接近希望的RNA浓度。大多数Poly(A)⁺mRNA分离过程都会导致得到较稀的RNA，所以需要在cDNA合成前浓缩mRNA。TRIzol分离和热酚提取后需要乙醇沉淀：1)加1/10体积3M NaOAc, pH5.2, 和2.5倍体积乙醇；2)混匀，-20℃放置最少1小时；3) 4℃，≥12000 × g离心20分钟；4) 80%乙醇洗涤沉淀2次；5)空气干燥沉淀；6) DEPC（焦磷酸二异质）处理水重新溶解沉淀。最合适的溶解体积由cDNA合成中需要的RNA的浓度和量来决定。RNA测定，用分光光度计分析RNA浓度，在260 nm下以1单位吸光度等于40 μg/ml RNA，需要在260 nm和280 nm测定吸光度来确定样品的浓度和纯度，A₂₆₀/A₂₈₀应接近2.0为较纯的RNA (即比值在1.9-2.1)。

第3步：由纯化的总RNA合成双链cDNA。

HPLC纯化T7-(d7)₂₄引物；

1) 第一链cDNA合成，cDNA合成前，DEPC处理水和逆转录的正确体积必须确定。它由加到反应中的RNA浓度和总体积决定。

RNA和SuperScript ⅡRT体积不要超过12 μl。

合成反应可在1.5 ml离心管中进行(RNase-free)，按下列组分合成cDNA。

2) 第二链cDNA合成。第一链反应放置冰上，稍微离心甩下管壁试剂，在第一链合成的管中加入下列第二链反应试剂，混匀。

大于等于12000 g离心10分钟，16℃放置2小时加2 μl 10 U T4 DNA聚合物16℃放置5分钟，加10 μl 0.5 M EDTA继续纯化cDNA步骤或-20℃储存。

3) 纯化双链cDNA。Phase Lock Gels (PLG)-酚/氯仿提取，≥12000 g离心PLG管20-30秒，离下管壁PLG；加162 μl（等体积）的（25：24：1）酚：氯仿：异戊醇（10 mM Tris-HCL pH8.0, 1 mM EDTA饱和）到cDNA最后合成产物中（162 μl），最后体积到324 μl，混匀，≥12000 g离心10分钟；转移上清液至PLG管，不要混合，PLG会混入溶液中；≥12000 g离心2分钟；转移上层水相到一个新的1.5 ml离心管中；加0.5倍体积7.5 M NH₄OA_C和2.5倍体积乙醇（-20℃储存）到样品中，混匀；立即在室温下≥12000 g离心20分钟；去上清，0.5 ml 80%乙醇（-20℃储存）洗涤沉淀，在室温下≥12000 g离心5分钟；小心去掉80%乙醇，80%乙醇再洗涤一次；空气干燥沉淀，检查是否干燥，Rnase-free水重新溶解沉淀。

第4步：生物素标记cDNA合成。

利用BioArray High Yield RNA Transcript Labeling kit进行cDNA标记，cDNA体外转录产物(总RNA)如下表。

cDNA体外转录成分如下，37℃，4.5小时，600 rpm振荡10秒/35分种。

第5步：纯化和质控体外转录(IVT)产物。

QIAGEN RNeasy Columns纯化体外转录产物，洗涤和洗脱之前将样品过柱两次；洗脱RNA时加水到柱子后，静置一分钟，再离心。cDNA质控，用分光光度计分析RNA浓度，A₂₆₀/A₂₈₀应接近2.0为较纯的RNA (即比值在1.9-2.1)。根据下面的计算公式确定调整cRNA的含量：

cRNA含量=RNA_m-(总RNA_i) (y)

RNA_m=体外转录后测得cRNA量 (μg)

总RNA_i=开始总RNA的量 (μg)

y=在IVT过程中使用的cDNA的倍数

凝胶电泳检测样品，同时进行纯化和没有纯化的体外转录产物的凝胶电泳有助于检测纯化过程丢失的范围，0.1%琼酯糖凝胶电泳分析0.1%的样品，RNA和溴化乙锭混合，加热到65℃，5分钟。

第6步：片段化cDNA。

在新的1.5 ml RNase-free离心管中按下表加入样品，

94℃，35分钟。然后放置冰上。变性凝胶电泳，至少需要1 μg cRNA。-20℃储存样品。

第7步：杂交。

在新的1.5mL RNase-free离心管中按下表加入样品。

20×真核生物杂交控制试剂冻存，在使用前在65℃下5分钟。使用前室温平衡探针；在99℃下5分钟；通过加样孔加入适量体积1×杂交缓冲液湿润芯片；在45℃下60 rpm预杂交芯片10分钟；处理过的样品45℃，5分钟；最大速离心5分钟；从芯片中取出缓冲液，加等体积处理好的杂交液45℃，60 rpm杂交芯片16小时。

第8步：洗脱和染色。

杂交16小时后，从芯片中取出杂交液装入一个新的离心管，放置冰上或-80℃长时间保存；洗脱缓冲液A充满芯片；配制下列溶液：

SAPE液（使用前配制，4℃储存）

成分	体积	终浓度
			2×MES染色剂缓冲液	600.0 μl	1×
50 mg/ml乙酰化BSA	48.0 μl	2 mg/ml
			1 mg/ml Streptavidin-Phycoerythrin(SAPE)	12.0 μl	10 μg/ml
去离子水	540.0 μl
			总体积	1200 μl

抗体溶液

洗脱工作站按下表工作。

第9步：扫描。

采用GeneChip Scanner 3000 7G扫描芯片，并进行芯片的图像分析处理。

（二）、Agilent-014850 Whole Human Genome Microarray 4x44K G4112F芯片实验过程如下：

第1步：总RNA的抽提，每2×10⁷细胞加入1 ml Trizol，在旋涡震荡器上混匀；用液氮研磨或采用电动匀浆器充分打碎组织块。加入约1/5体积的氯仿，上下颠倒充分混匀1分钟左右，室温下静置5分钟。4℃，12,000 rpm离心15分钟后小心取出上清液，将上清夜转入新的1.5 ml离心管，加入等体积的异丙醇，轻轻颠倒混匀，室温静置5 分钟。（15 ml离心管用7,500 rpm离心20 min）4℃，12000 rpm离心10分钟后，去上清，向沉淀中加入2/5体积的70%乙醇，4℃，12000 rpm离心洗涤沉淀15分钟。（15ml离心管用7,500 rpm离心20 min）去上清，沉淀室温晾干后加入适量无RNA酶的水充分溶解沉淀，测定OD₂₆₀和OD₂₈₀值。

第2步：总RNA质量检测。

1) 琼脂糖凝胶电泳，配置用DEPC处理的电泳缓冲液50×TAE，高压灭菌后待用。使用电泳槽前用3％H₂O₂浸泡15 min，然后用DEPC处理的水冲洗，倒入适量1×TAE电泳缓冲液。称取适量琼脂糖，加入1×TAE电泳缓冲液，制备1％的胶（注意使用专用的溶液和相关设备，避免引入外源RNA酶）。用专用6×Loading缓冲液做指示剂，取10 μl上样电泳（电压100伏）15 min。关闭电源，取出电泳胶，在凝胶成像仪上观测、拍照，保存图像。评价总RNA或mRNA质量，通过测28S和18S的亮度比例可以初步评价总RNA的质量。一般28S:18S≥2可以初步判定总RNA质量较好。

2) lab-on-chip，胶制备取出400 μl RNA凝胶基质，加入Spinfiter柱子，离心过滤凝胶（1500 × g，10 min）。取出130 μl过滤胶到1.5 ml离心管中，再加入2 μl RNA Dye Concentrate，在涡旋仪上振荡混匀。用RNA ZAP清洁操作区，同时在Electrode cleaner中加入350 μl ZAP，放入Lab on chip正确位置，合盖清洁探头1 min。再用另外一Electrode cleaner中加入350 μl DEPC处理的H₂O，重复。移开Electrode cleaner,让探头自然干燥。取出一新的RNA Chip，吸取9 μl步骤2制备的凝胶加入G孔中。将chip放置于带有活塞的水平台上（chip priming station），将活塞拉杆向上拉倒1 ml刻度出，再将chip priming station的盖子合上，压紧chip。同时向下推动拉杆至底部并维持30秒左右，松开让拉杆自动弹开。从chip priming station上取出chip，用放大镜检查是否有气泡存在，如果在微通道中有气泡，需要重复步骤。再在标有G的两个孔中各加入9 μl步骤2制备的凝胶，加入5 μL rna 6000 Nano Marker。同时在12个加sample的孔中滴加5 μl RNA 6000 Nano Marker（不能空一个孔）。滴加1 μl RNA 6000 ladder，再各取1 μl样品（RNA）加入12个样品孔中，并将chip放入IKA涡旋仪上，在set-point振幅处振动1 min。振荡好后的chip在5 min之内必须放入Agilent 2100分析仪上，按照软件提示进行RNA电泳操作。

第3步：总RNA的纯化。

使用QIAGEN RNeasyKit纯化总RNA，取总RNA ≤ 100 μg溶解于100 μl RNase free水中，加入350 μl Buffer RLT并充分混匀。加入250 μl无水乙醇，Tip头充分混匀。将共计700 μl含总RNA的溶液转入套在2 ml离心管内的RNeasy柱子内，≥8000 g离心15-30 sec，弃去滤过液。吸取500 μl Buffer RPE到RNeasy mini 柱子内，≥ 8000 g离心洗涤15-30 sec，弃去滤过液，再用500 μl Buffer RPE在≥8000 g离心洗涤2 min，弃去滤过液和2 ml的套管，将RNeasy mini柱子转入一新的1.5 ml Eppendorf管中。吸取40 μl RNase free的水，≥8000 g离心洗脱1 min。重复上一步骤。

第4步：cDNA双链合成。

取2 ug总RNA于1.5 ml离心管中，如下配置反应溶液：

总RNA 2 ug

T7启动子引物5 ul

无核糖核酸酶的水X ul

65℃保温10分钟，冰浴5分钟。配置如下cDNA合成体系：

5×第一链缓冲液4 ul

0.1M DTT2 ul

10mM dNTP mix1 ul

MMLV-RT1 ul

核糖核酸酶 OUT 0.5 u

	2.2r(ul)	3.3r(ul)	4.3r(ul)	5.3r(ul)	6.3r(ul)	7.3r(ul)	8.3r(ul)	9.3(ul)	10.3(ul)
										5×第一链缓冲液	8.8	13.2	17.2	21.2	25.2	29.2	33.2	37.2	41.2
0.1M DTT	4.4	6.6	8.6	10.6	12.6	14.6	16.6	18.6	20.6
										10mM dNTP mix	2.2	3.3	4.3	5.3	6.3	7.3	8.3	9.3	10.3
MMLV-RT	2.2	3.3	4.3	5.3	6.3	7.3	8.3	9.3	10.3
										核糖核酸酶-OUT	1.1	1.65	2.15	2.65	3.15	3.65	4.15	4.65	5.15

将上述混合液加入变性后冰浴的RNA中混匀，40℃保温2 h。

第5步：荧光标记cDNA合成。

cDNA合成结束后将离心管放65℃保温15分钟，冰浴5分钟。

如下配置转录本混合物：

无核糖核酸酶的水5.7 ul

4×转录本缓冲液20 ul

0.1 M DTT6 ul

NTP Mix(ATP、GTP、CTP、UTP) 16 ul

50% PEG6.4 ul

Rnase OUT0.5 ul

无机焦磷酸酶0.6 ul

aa-UTP(25 mM) 4 ul

T7 RNA 聚合酶0.8 ul

加入60 ul Transcription mix并混匀，40℃保温2 h，冰上放置。

第6步：cDNA纯化。

采用QIAGEN RNeasy Mini Kit进行cDNA纯化。加入20 ul Rnase-free Water, 加入350 μl Buffer RLT并充分混匀。加入250 μl无水乙醇，Tip头充分混匀。将共计700 μl含总RNA的溶液转入套在2 ml离心管内的RNeasy柱子内，≥8000 g离心15-30 sec，弃去滤过液。吸取500 μl Buffer RPE到RNeasy mini柱子内，≥8000 g离心洗涤15-30 sec，弃去滤过液，再用500 μl Buffer RPE在≥8000 g离心洗涤2 min，弃去滤过液和2 ml的套管，将RNeasy mini柱子转入一新的1.5 ml Eppendorf管中。吸取30 μl RNase free的水，静置1 min，≥8000 g离心洗脱1 min。吸取30 μl RNase free的水，静置1 min，≥8000 g离心洗脱1 min。

第7步：cDNA探针后标记。

取4ug上述cRNA并浓缩至6.66ul，加10ul DMSO混匀，加3.4ul 0.3M 碳酸氢钠缓冲液, pH 9.0并混匀，将上述20ul cRNA混合物加入到荧光染料中并混匀，室温25度保温1小时，加9ul 4M羟胺混匀后室温即25度保温15分钟。

第8步：cDNA探针纯化。

加入70ul Rnase-free Water, 加入350 μl Buffer RLT并充分混匀。加入250 μl无水乙醇，Tip头充分混匀。将共计700 μl含总RNA的溶液转入套在2 ml离心管内的RNeasy 柱子内，≥8000 g离心15-30sec，弃去滤过液。吸取500 μl Buffer RPE到RNeasy mini柱子内，≥8000 g离心洗涤15-30sec，弃去滤过液，再用500 μl Buffer RPE在≥8000 g离心洗涤2 min，弃去滤过液和2 ml的套管，将RNeasy mini柱子转入一新的1.5 ml Eppendorf管中。吸取 30 μl RNase free的水，静置1 min，≥8000 g离心洗脱1 min。吸取30 μl RNase free的水，静置1 min，≥8000 g离心洗脱1 min。

第9步：芯片杂交。

cRNA的片段化：

Cy3 cRNA X ul（100 pmol）

Cy5 cRNA X ul（50 pmol）

10× 控制靶标50 ul

无核糖核酸酶的水X ul

总体积 240 ul

加入片段化缓冲液10 ul，轻轻混匀后60℃保温30分钟（不超过30分钟）。加入250 ul杂交缓冲液轻轻混匀离心。取490 ul杂交溶滴加在盖玻片上，盖上芯片并密封芯片于杂交盒中，60℃滚动杂交16小时。

第10步：芯片洗涤。

洗液配置：

洗液1（1升）：焦磷酸二异酯处理过的水700 ml

20×SSPE300 ml

20%十二烷基肌氨酸钠0.25 ml

洗液2（1升）：焦磷酸二异酯处理过的水997 ml

20×SSPE 3.0 ml

20%十二烷基肌氨酸钠0.25 ml

洗液3：稳定液和干燥剂

取出芯片于洗液1中洗涤1分钟；再将芯片放入洗液2中洗涤1分钟；最后将芯片放入洗液3中洗涤30秒。

第11步：芯片扫描。

采用Agilent SureScan Microarray Scanner扫描芯片，并进行芯片的图像分析处理。

（三）、跨平台数据整合

Affymetrix HG U133A芯片平台共包含22283个探针，对应13434个不同的Entrez GeneID，平均1个Entrez Gene ID对应1.66个探针；而Agilent 014850 Whole Human Genome Microarray 4x44K G4112F包含43376个探针，对应19596个不同的Entrez GeneID，平均一个基因对应2.214个探针；两款芯片平台共有的基因（GeneID或Symbol）有12,810个，这意味着95.35%的HG U133A检测到的基因都可以在Agilent G4112F中找到与之相对应的表达值。分别扫描20例人类乳腺癌细胞系的Affymetrix芯片和Agilent芯片得到包含有20个样本的HG U133A的基因表达谱和G4112F基因表达谱数据。

对于上述两个芯片平台间的表达谱数据按照下面方法步骤进行跨平台数据整合：

第1步：对基因芯片的探针荧光值数据预处理、背景矫正和中位数标准化。Affymetrix和Agilent的芯片均采用RMA方法进行数据预处理、背景矫正，并进行中位数标准化。Affymetrix芯片数据处理采用R统计软件中的affy包rma方法，Agilent芯片数据处理采用R统计软件中的preprocessCore包rma.background.correct方法，芯片间中值标准化采用的是limma包normalizeBetweenArrays方法。从Affymetrix和Agilent提供的芯片探针注释文件中，得到芯片探针和基因的对应关系，HG U133A平均每个基因具有1.66个探针重复，Agilent G4112F平均每个基因具有2.214个探针重复，为了计算各个芯片平台下基因的表达信息，当一个基因在芯片中对应2个以上探针时，以探针荧光值的均值表示该基因的表达。基于上面的规则将得到的探针荧光值统一转化为基因表达值，Affymetrix HG U133A检测的20例样本就生成了1个13,434

20的基因表达谱矩阵，相应地Agilent G4112F芯片检测的20例样本就生成了1个19,596

20的基因表达谱矩阵，其中表达谱矩阵的行表示每一个基因，列表示20个被检测的乳腺癌细胞系样本。

第2步：合并两款芯片间共同存在基因数据，生成了一个12,434

40基因表达谱矩阵，前20列为Agilent检测平台数据，后20列为Affmetrix检测数据。采用k均值聚类将12,434基因聚为100（m=100）类，如图2所示，两个芯片的基因表达值在某一基因子集范围内，呈现显著正相关关系。

第3步：利用最小二乘法线性拟合

（需满足exp_s1>0同时exp_s2>0），其中exp_s1是由基因子集s在芯片平台1中的表达值均值构成的数组，exp_s2指基因子集s在芯片平台2中相应的的表达值均值构成的数组，h为基因子集s包含的基因数目，s属于1--m(此处m为100)，代入最小二乘法线性拟合公式

，，计算100个基因子集的a_s和b_s。

第4步：将第3步得到的a_s和b_s代入公式

，Affymetrix和Agilent的基因表达值均标准化到0-12范围内波动，如图3所示。但也可以看出两款芯片间样本基因表达的中值及四分位数具有明显的区别。

第5步：全局中值标准化。将第4步得到的基因表达谱矩阵，利用公式

进行基因表达的全局中值标准化，其中m_j指表达谱矩阵第j列的中值，m为整个矩阵的中值。标准化后的结果如图4所示。

第6步：样本间方差标准化。第5步得到的基因表达谱矩阵如图4所示，Affymetrix和Agilent的基因表达值中值处于同一水平，但基因在各个芯片的方差变化存在显著地差异，同一样本或同一组样本的虽然采用了不同的芯片检测平台，但他们的基因表达数据表达应服从同一分步。利用公式

，其中

指表达谱矩阵第j列的均值，值矩阵第j列数据的标准差，芯片内标准化样本的基因表达方差，结果如图5所示。

第7步：基于非参数经验贝叶斯模型，利用R统计软件中的ComBat方法，校正Affymetrix和Agilent芯片间的实验批次影响，最终跨平台数据整合结果如图6所示。

跨基因芯片平台的数据整合目标是挖掘公开的数据源中与自身研究目的相似的数据集，并对这部分数据进行整合分析，其中最大的难题是基因在多个平台的表达一致性问题。同组样本或相似的样本即使采用不同的芯片平台测量基因表达，相同基因在两个或多个芯片平台的表达值也应具有较高的相关性，但由于实验中各个环节及平台探针设计的不同等因素，共同存在基因的表达在两个平台间的总体相关性并不高。如图7所示，经过本发明方法对两个芯片平台的数据整合后，同一样本相同基因在两个芯片平台间表达值的相关性有了显著地提高（整合前20例样本在两个平台得到的基因表达值的总体相关性R²为0.45，整合后得到的总体基因表达相关性R²为0.80），可以将整合后的数据作为一个整体进行后续的乳腺癌肿瘤分子研究。

Affymetrix HG U133A检测的10例乳腺癌组织样本和10例癌旁组织样本经第1步和第2步处理得到的基因表达谱，t-test检验统计13,434个基因在癌和癌旁组织中的差异表达分析，并进一步采用BH多重比较校正方法对t检验结果p进行校正得到q，以q<=0.05为阈值，Affymetrix HG U133A共检出794个基因在乳腺癌组织和癌旁组织中显著差异表达。Agilent G4112F检测的10例乳腺癌组织样本和10例癌旁组织样本经第1步和第2步二处理得到的基因表达谱，t-test检验统计19,596个基因在癌和癌旁组织中的差异表达分析，并进一步采用BH多重比较校正方法对t检验结果p进行校正得到q，以q<=0.05为阈值，Agilent G4112F共检出1,376个基因在乳腺癌组织和癌旁组织中显著差异表达。比较Affymetrix检测到的794个差异表达基因和Agilent检测到的1,376个差异表达基因，得到719个共同的差异表达基因。t检验统计第7步得到的合并后基因表达谱数据，进一步采用BH多重比较校正方法对t检验结果p进行校正得到q，以q<=0.05为阈值，得到812个显著差异表达基因。合并后基因表达谱数据显然得到了比分别显著表达的基因的交集更多的显著差异表达基因。812个差异表达基因和719个差异表达基因的交集为643，跨芯片平台数据整合分析方法挖掘了到89.4%的差异表达基因，此外还有163个差异表达基因是跨平台整合后发现的，此163个基因在Affymetrix HG U133A和Agilent G4112F实验中以q<=0.05均为检出差异表达，并且有文献表明着163个基因中的部分基因参与的血管新生、调控细胞分析和调控细胞程序性凋亡等功能与乳腺癌发病机制相关，见下表：

GO term	Count	GeneIDs	GeneNames
				GO:0001525血管新生	5	302; 8742; 83478; 14; 80155	门冬氨酸受体调节蛋白1； Rho GTP酶激活蛋白24；肿瘤坏死因子（配体）超家族成员13；血管相关的迁移细胞蛋白；膜联蛋白A2假基因1
GO:0045595调控细胞分化	10	23318; 19; 80070, 100133941; 2103, 7068; 9607; 4897; 8359; 5458	锌指CCHC结构域包含11；ATP结合盒子家族成员1；ADAM金属肽酶血小板粘合素类型1 motif 20；CD24分子；雌激素相关受体beta；甲状腺激素受体beta；CART神经肽；神经细胞粘附分子；组蛋白类H4a；POU4同源框2
				GO:0010941调控细胞的程序性死亡	13	2903; 7015; 2065; 2923; 80070; 100133941; 8742, 308; 350; 3070; 4891; 2895; 29949	谷氨酸受体2A；端粒酶反转录酶；禽白血病病毒癌基因同源物3；蛋白质二硫键异构酶家族A成员3；ADAM金属肽酶血小板粘合素类型1 motif 20；CD24分子；肿瘤坏死因子（配体）超家族成员12；模联蛋白A5；脂蛋白（β2糖蛋白I）；淋巴细胞特异性解旋酶；溶质载体家族11成员2；谷氨酸受体D2；白细胞介素19

上面的跨平台整合结果说明本方法不仅可以满足疾病分子机制研究的需求，同时还可以满足增大研究统计样本量提高分析精度的需求。

实施例二、以公开数据库Gene Expression Ominbus中获取的人类肝癌和正常组织基因表达数据GSE14520-GPL571和GSE46408为例阐述本发明：

GSE14520-GPL571包含19个人类的肝癌组织样本和19个对应非肿瘤组织样本的基因表达数据，基因表达数据集的ID为GSM362950、GSM362951、GSM362952、GSM362953、GSM362954、GSM362955、GSM362956、GSM362957、GSM363420、GSM363421、GSM363422、GSM363423、GSM363424、GSM363425、GSM363426、GSM363427、GSM363428、GSM363429、GSM363430、GSM363431、GSM363432、GSM363433、GSM363434、GSM363435、GSM363436、GSM363437、GSM363438、GSM363439、GSM363440、GSM363441、GSM363442、GSM363443、GSM363444、GSM363445、GSM363446、GSM363447、GSM363448和GSM363449，基因表达检测所用的芯片平台为Affmetrix HG U133A 2.0。GSE46408包含6个人类的肝癌组织样本和6个对应非肿瘤组织样本的基因表达数据，基因表达数据集的ID为GSM1129834、GSM1129835、GSM1129836、GSM1129837、GSM1129838、GSM1129839、GSM1129840、GSM1129841、GSM1129842、GSM1129843、GSM1129844和GSM1129845，基因表达检测所用的芯片平台为Agilent G4112F。

第1步：由于GSE14520-GPL571的38个芯片数据已经过RMA方差标准化，公共数据库下载的探针表达谱矩阵数据可以直接用于后续的分析；获取GSE46408的12个芯片原始数据，采用R统计软件中的preprocessCore包rma.background.correct方法，芯片间中值标准化采用的是limma包normalizeBetweenArrays方法，得到标准化后的探针信号值表达谱矩阵。从Affymetrix和Agilent提供的芯片探针注释文件中，得到芯片探针和基因的对应关系，当一个基因在芯片中对应2个以上探针时，以探针荧光值的均值表示该基因的表达。基于GSE14520包含的38例样本就生成了1个13,434

38的基因表达谱矩阵，相应地GSE46408包含的12例样本就生成了1个19,752

12的基因表达谱矩阵，其中表达谱矩阵的行表示每一个基因，列表示被检测人类肝癌或正常肝脏组织样本。

第2步：合并两款芯片间共同存在基因数据，生成了一个12,809

50基因表达谱矩阵(n=12809，k=50)。采用层次聚类方法，将12,809基因分到了102个聚类树中，得到102个基因子集（m=102）。

第3步：利用最小二乘法线性拟合

（需满足exp_s1>0同时exp_s2>0），其中exp_s1是由基因子集s在芯片平台1中的表达值均值构成的数组，exp_s2指基因子集s在芯片平台2中相应的的表达值均值构成的数组，h为基因子集s包含的基因数目，s属于1--m，代入公式，

计算102个基因子集的a_s和b_s。

第4步：将第3步得到的a_s和b_s代入公式

，GSE14520和GSE46408基因表达值标准化为在0-14范围内变化。

进行基因表达的全局中值标准化，其中m_j指表达谱矩阵第j列的中值，m为整个矩阵的中值。标准化后得到50个样本的全部12,809基因的表达值中值均为4.886。

第6步：样本间方差标准化。第5步得到的基因表达谱矩阵，利用公式

，其中

指表达谱矩阵第j列的均值，值矩阵第j列数据的标准差，芯片内标准化样本的基因表达方差。得到50个样本的全部基因表达值服从均值为0方差为1的正态分布。

第7步：基于非参数经验贝叶斯模型，利用R统计软件中的ComBat方法，校正GSE14520和GSE46408芯片间的实验批次影响，最终跨平台数据整合结果如图8所示，50个样本的的基因表达总体上处于同一水平，波动范围也类似，使得后续的基因表达分析具有了可比性。

GSE14520数据集中检测的19例肝癌组织样本和19例配对的正常肝脏组织样本经第1步和第2步处理得到的基因表达谱，采用配对t-test检验统计13,434个基因的差异表达，并进一步采用BH多重比较校正方法对t检验结果p进行校正得到q，以q<=0.05为阈值，GSE14520共检出1,770个基因在肝癌组织和正常肝脏组织中显著差异表达。Agilent G4112F检测的GSE46408的6例肝癌组织样本和6例配对正常肝脏组织样本芯片原始数据经第1步和第2步处理得到的基因表达谱，配对t-test检验统计19,752个基因的差异表达分析，并进一步采用BH多重比较校正方法对t检验结果p进行校正得到q，以q<=0.05为阈值，GSE46408共检出2,580个基因在肝癌癌组织和正常组织中显著差异表达。比较Affymetrix检测到的1,770个差异表达基因和Agilent检测到的2,580个差异表达基因，得到1,178个共同的差异表达基因。配对t检验统计第7步得到的合并后基因表达谱数据，进一步采用BH多重比较校正方法对t检验结果p进行校正得到q，以q<=0.05为阈值，得到1,690个显著差异表达基因。合并后基因表达谱数据显然得到了比分别显著表达的基因的交集更多的显著差异表达基因。1,690个差异表达基因和1,178个差异表达基因的交集为1162，即98.6%的差异表达基因都可以在跨平台整合的基因表达数据中被挖掘，可见本发明方法对于公开的单通道芯片数据依然有较好的效果。医学统计学上以a=0.05为显著性阈值进行t检验统计分析时要求样本应不小于30，对于GSE46408这套基因表达数据仅包含了12个样本，没有达到统计学上对最小样本数的要求；但经过本方法芯片数据整合后，该研究项目的样本数达到50，满足了统计学对样本量的要求。

以上结果说明利用本发明的方法成功实现了不同平台的基因表达谱数据的整合。

以上是对本发明的描述，而非限定，基于本发明的思想的其他实施方式，均在本发明的保护之中。

Claims

1.一种跨芯片平台的基因表达数据整合方法，所述方法包括如下步骤：

步骤1：获得要整合的两个或多个基因芯片平台表达数据，对于每个基因芯片平台而言，表达数据为一个矩阵，矩阵的每一行表示一个基因，矩阵的每一列表示一个被检测的样本；

步骤3：基于基因在所述合并的基因表达矩阵中的表达模式相似性将基因分割为多个基因子集；

步骤4：对于步骤3得到基因子集，利用基因子集内基因在所述每个不同芯片平台的表达值，得出不同芯片平台的表达线性关系，利用所述关系将所述不同芯片平台的基因表达值标准化到同一变化范围，得到标准化的基因表达矩阵。

2.权利要求1的方法，所述方法包括如下步骤：

步骤3：基于基因在所述合并的基因表达矩阵中的表达模式相似性将基因分割为m个基因子集；

，

，其中h为基因子集s包含的基因数；利用公式

将所述不同芯片平台的基因表达值标准化到同一变化范围，得到标准化的基因表达矩阵，

”表示两个矩阵相乘，“

”表示两个矩阵的点乘。

3.权利要求2的方法，所述方法包括如下步骤：

步骤1：获得要整合的两个或多个基因芯片平台表达数据，所述数据为基因芯片的探针荧光值，代表各个基因表达，将每个芯片基因平台的探针荧光值表示对应基因的表达值，即基因转录本多少的定量值，全部基因(Symbol 1-n)在全部样本(sample 1-k)中的表达值(exp11-nk)构成即基因表达矩阵(exp)，矩阵的每一行表示一个基因，矩阵的每一列表示一个检测样本，基因可采用任何标识形式表示，表达矩阵如下：