本发明涉及一种复杂数据预测模型的构建方法。
背景技术：
新一代测序技术的出现，使研究人员能够处理收集的大数据(例如，使临床研究人员能够处理收集自患者的上百个生物样本)，并进行如全基因组表达水平、甲基化水平或体细胞突变的分析，这里称为高维组学数据(hdod，highdimensionomicsdata)。虽然可获得的临床样品量通常有限，但由于每个样本被观测的变量的数目可以达到数千或数百万，因此临床研究的瓶颈，已经从样品采集转移到了数据管理和数据分析上。利用hdod连同其它临床变量建立特定临床结果的预测模型，已经是生物医学信息学的研究人员的众多分析目标之一。
建立预测模型已经成为一些学科的定量研究员共享的研究点。研究员一直在积极利用来自数据库的大数据集进行预测模型的开发，采用的方法包括机器学习算法、支持向量机和遗传算法。此外，基于对数据库技术和可视化工具的熟练掌握，研究员可以有效地构建hdod，通过缩放p计算分析hdod，并使得hdod衍生的结果可视化，从而使生物医学研究人员可以对hdod进行处理，并可以直观地观测结果。
构建预测模型已经是现有技术，通常是根据已知预测多变量的结果，构建基于回归的预测模型，且大多采用广义线性模型(glm)。hastie和tibshirani放宽了参数假设，描述了广义相加模型(gam)，用非参数回归方法结合几十年的研究。近年来，统计学家一直在研究使用惩罚似然技术(包括lasso、gbm和弹性网络技术)来使hdod自动的选择协变量。这些方法是转化研究中用于处理维度灾难的主要工具。
虽然计算机科学与统计学之间存在交叉，但之间的基本区别在于，计算机科学通常从系统的角度探索带有多变量的图谱，而统计学倾向于遵循节省原则确定几个协变量预测模型。统计学面临的一个主要挑战是如何控制根据hdod选择预测器的假阳性错误率的过度膨胀，其将导致“过度拟合”预测模型。与此相反，计算机科学或生物信息学中，则主要对hdod图谱感兴趣，常常想要量化直观的图谱，重复生成图谱独立的数据集。
本发明保留这两种分析方法的特点，提出一种混合算法，包括两个步骤：在第一步骤中，要确定一组代表对象hdod图谱的“范例”，所述“范例”一般通过无监督学习的聚类分析法获得。为了代表集群图谱，选择单一集群的质心作为范例。每个范例通过p个元素的hdod向量分类。范例的数目(q)通常小于等于样本量(n)。参照各范例，可以计算每个对象的相似性度量，生成具有维度(n×q)的相似性度量的矩阵，通常情况下，p>>n≥q。本步骤可以有效地将高维稀疏矩阵hdod(n×p)转换成“稠密数据矩阵”(n×q)。在第二步骤中，使用惩罚似然方法来选择出那些符合预测结果的范例。由于维数从p大幅减小到q，惩罚似然方法可以很好地选择出包含信息的范例，大大减少了惩罚计算的步骤。本过程首先基于“无监督学习”的范例，然后通过“有监督学习”选择与结果关联的包含信息的范例。由于结果回归范例特异性的相似性，这种方法被称为“面向对象的回归”，或简称为oor。
随着新一代测序技术，一些生物技术学家/生物技术公司已经将其创新研究转向于生产人类基因组的完全分相二倍体，即，一对带有多个单核苷酸多态性(snps)的分相单倍型。在功能基因内，多个分相snp等位基因，连同所有单型核苷酸，代表可用于破译官能转录物或蛋白序列的完全分相序列。实际上，目前这种双等位基因多态性可以构建成多等位基因多态性，能对遗传分析提供更多的信息。最好的范例基因包括主要用于组织相容性复合体(mhc)的人类白细胞抗原(hla)基因，位于染色体6上的6p22.1和6p21.3之间。例如，下面将要详述的hla*drb1基因，由一对等位基因组成，各等位基因对应一分相序列。根据最近的计数统计(http://www.ebi.ac.uk/ipd/imgt/hla/)，hla*drb1拥有超过1868个等位基因，编码1364个蛋白质。虽然对它们的功能已经进行了几十年的研究，但其特殊的多态性对如i型糖尿病(t1d)等疾病的遗传关联性的研究，则是个新的课题。此外，因为与许多较不常见的等位基因关联的样本数量有限，以及多个等位基因测试的多重性的原因，阻碍了多态性从基础研究到临床应用的转化。
为了克服上述问题，需要一个新的分析框架。在大多数科学事业，如遗传学中，通常采用简化论方法进行分析，即专注于与单一染色体、基因、等位基因或核苷酸的疾病的关联性。但这种简化论方法在同时处理太多元素时受到挑战，不适于用于同时处理太多的元素。近年来组学研究中，越来越多的科学团体开始关注多基因及其与表型联合关联的“系统生物学”，即“整体”的方法。从“整体”方法的角度看，当两个对象共享相同的疾病表型时，可能是因为两者有相似的基因分布(所述基因分布是基于多基因的基因型进行分类的)，这促进了上述面向对象的回归(oor)的应用和发展。
oor的核心思想如下简述：基于一组选定的基因/snp，构造一组以多基因/snp的基因型分布为特征的范例。通过比较个体的基因型和范例，有效地将基因型的域转化为相似性值的域。通过这些相似性度量，oor评估疾病表型是否与每个范例的相似性度量相关联。如果发现范例的相似性与表型显著关联，意味着该范例的基因型代表一种风险/保护基因型的类别。
oor实际上与一些统计文献以及计算机科学中的机器学习文献中的方法存在关联。从根本上讲，如果所有范例的基因型分布是由内部衍生，并包括所有观测到的基因型，可以认为上述oor是基于kimeldorf和wahba(1971)所描述的表现定理中的核表示进行的改进。基于同样的类比，oor与核逻辑回归和支持向量机也存在密切联系。联系到计算机科学的文献，由于密切相关的“邻居”倾向于拥有相同的表型，oor与近邻法拥有相同的动机。在处理复杂、稀疏、高维的数据时，通过“相似性度量”搜索数据库，对数据挖掘是至关重要的。近年来，统计和数据挖掘的融合促进了核机器学习技术在文本挖掘、蛋白质序列分析以及全基因组关联分析方面的应用。
尽管oor与核机器方法紧密联系，但仍有区别。首先，比核机器学习方法先进的是，oor的范例可以从外部获得或从内部数据衍生。其次，通过把所有计算得到的相似性度量作为协变量，oor通过惩罚似然法使用“变量选择算法”，如lasso、脊回归或弹性网络，来侧重于不同于零的有意义的项。第三，以“整体观”对待复杂的变量，oor提供了一种天然量化工具来发现和验证复杂的变量之间的相互作用，所述复杂的变量之间的相互作用已成为在生物医学研究和系统生物学的一个长期的研究课题。最后，从oor分析得到的预测模型很适合于将基于相似性的搜索应用到大型数据库。
在下文中，本发明第一部分示出了oor的统计学动机，勾画出了oor框架，确定了选择范例的方法，并构建出预测模型。此外，本发明还介绍了从协变量转换到相似性度量，然后建立预测模型的部分流程。除了详细介绍了对范例以及预测器的选择，还介绍了如何评估选择惩罚参数的稳定性以及如何通过自助法评估所含信息范例的一致性。为了说明oor，应用部分介绍了i型糖尿病的研究，并说明了oor在探索疾病与hla基因的关联以及构建预测模型的应用。结果部分介绍了应用于hla-drb1基因以及八个hla基因的所有结果。
技术实现要素：
为了解决现有技术的不足，本发明提供了一种复杂数据预测模型的构建方法。
本发明提供的一种复杂数据预测模型的构建方法，所述方法包括以下步骤：
a、获取一高维组学数据hdod，确定出对于所述hdod具有代表性的一组数据对象作为范例；
b、确定出所述hdod中每个数据对象与每个所述范例的相似性度量，并据此构建数据对象与范例的相似性度量矩阵；
c、由所述数据对象与范例的相似性度量矩阵，通过惩罚似然法从所述范例中选出包含信息的范例；
d、基于所选出的范例构建预测模型。
可选的，步骤a所述具有代表性的一组数据对象包括：可分别代表hodo的不同组数据的各特性的各数据对象。
可选的，步骤a所述具有代表性的一组数据对象包括：具有远离集群、相对独特性的hodo特性的数据对象，即不易由质心或它们的组合来表示其hdod特性的对象。
可选的，所述独特性的对象用下面的回归方法来确定：
假设预先已确定了一组起始的t个质心作为范例，表示为[1],[2],...,and[t]，
首先，通过公式
然后，针对每个xi估算来自上述线性回归的残差平方的总和(srsi)，并且计算由这些信息代表的残差变化的分数；当第i个对象xi满足公式i＝argmax(srsi/srs0),subjectto(srsi/srs0)≥f时，其加入到所述范例中；其中srs0是不包含范例的srs，f是一预先选择的阈值。
可选的，步骤a所述确定出对于所述hdod具有代表性的一组数据对象作为范例的步骤包括：基于面向对象的回归oor方法通过无监督聚类分析来对hdod确定出所述一组数据对象作为范例。
可选的，步骤b所述相似性度量的确定步骤包括：所述相似性度量通过相似性度量函数进行计算时，相似性度量的选择取决于所述hdod的特性和/或对于要回归的目标结果的解释。
可选的，步骤c包括：通过惩罚似然法，由所述相似性度量矩阵的各个范例对应的各相似性度量，计算出的各个范例对应的回归系数值，选择该回归系数值与设定阈值区间匹配的回归系数值所对应的范例为包含信息的范例。
可选的，步骤c所述惩罚似然法为lasso算法。
可选的，步骤c所述各个范例对应的回归系数值的计算采用下述公式计算：
其中，f表示密度函数，yi是对应第i个范例的要回归的结果，si是对应第i个范例的相似性度量值，α是截距，β是回归系数，n是当前范例对应的样本的个数，q是所述范例的个数，λ是调谐参数。
本发明方法的范例可以从外部获得或从内部数据衍生。其次，通过把所有计算的相似性度量作为协变量，本发明方法可通过惩罚似然法使用“变量选择算法”，如lasso、脊回归或弹性网络，来侧重于不同于零的有意义的项。第三，以“整体观”对待复杂的变量，本发明方法提供了一种天然量化工具来发现和验证复杂的变量之间的相互作用，所述复杂的变量之间的相互作用已成为在生物医学研究和系统生物学的一个长期的挑战。最后，根据本发明方法得到的预测模型很适合于通过基于相似性的搜索应用到大型数据库。
附图说明
图1示出了面向目标的回归的流程图，其中a)协变量矩阵的高维组学数据(hdod)，b)通过无监督学习方法组织hdod，c)通过双向聚类分析聚类的hdod,以进行范例的确定，d)计算每个范例的相似性度量，将它们当作协变量，e)相似性度量的稠密协变量矩阵，可用于构建预测模型，f)在广泛线性模型下，使用惩罚似然来选择包含信息的范例，g)在训练集和验证集上进行roc分析，以检查预测模型的有效性；
图2示出了训练集和验证集的所有患者的诊断年龄分布；
图3示出了通过对数秩检验法，使用kaplan-meier曲线探索年龄、性别、肿瘤类型和分期与5年内预后存活的边际关联，其中左栏为训练集，右栏为验证集；
图4示出了度量基因表达与i期指示器的关联的对数p值的估算分布；
图5示出了训练集中296个对象中被观测到的789个基因的高维组学数据(hdod)；
图6示出了训练集中2506个基因的高维组学数据(hdod)的选择的22个范例；
图7示出了训练集中789个基因的高维组学数据(hdod)的选择的11个范例；
图8示出了1000次蒙特卡洛模拟分析的估算惩罚参数值的经验分布，其中所选择的惩罚值(λ＝0.021)已标出；
图9示出了20个固定惩罚值的多次lasso的变量选择的估算(对数值示于对角线上)，其中每一个xy图展示了2个不同lasso的估算评价系数；
图10示出了oor分析的范例特异性预测模式下以及传统回归分析的协变量特异性预测的估算对数相对风险，其中左栏为训练集，右栏为验证集；
图11示出了预测模型下训练集的预测风险评分的分布(左栏)，风险评分1(参考)、2、3和4的估算存活率；
图12示出了训练集中705个对象的计算的相似性矩阵(705×705)，各元素取值分别为0(绿)、0.5(黑)和1(红)，以说明成对对象共享0个等位基因、1个等位基因和两个等位基因；
图13示出了通过lasso、脊回归、弹性曲线和逐步回归选择的范例的预测模型在训练集(实线)和验证集(虚线)估算的敏感度、1-特异性和曲线下面积(auc)，其中图中右侧彩色条为各个模型下相应的风险评分值(
图14示出了估算的各个元素的相似性矩阵，度量hla-drb1,drb345,dqa1,dqb1,dpa1anddpb1的未加权状态同源，其中颜色范围从绿到黑再到红，对应于备注中的低、中、高相似性；
图15示出了训练集(顶部图)和中验证集(底部图)的ii类hla基因(hla-drb1，-drb345，-dqa1，-dqb1，-dpa1和-dpb1的t1d预测模型的评价，箱图显示了训练集和验证集的风险评分分布，roc曲线示于左手区域；
图16示出了通过不同的惩罚参数下图谱偏差函数的反复交叉验证估算(上部区域的图)得到的估算惩罚参数
图17示出了1000个自助样本的平均系数估算的成对xy坐标图，其中x轴为一个惩罚值，y轴为另一惩罚值(y轴)，惩罚参数对数值示于对角线；
图18示出了当惩罚参数固定为15个对数独特系数之一时，通过lasso选择的范例的所有预测模型的roc分析与选择的范本由lasso，计算在训练集(彩色曲线)以及在验证集(黑色虚线)的auc值；
图19示出了惩罚参数固定为
具体实施方式
下面结合具体实施例进一步描述本发明，以更清楚的阐述本发明的优点和特点。下述实施例仅为具体的范例，并不对本发明的保护范围构成任何限制。本领域技术人员应该理解的是，在不偏离本发明的精神和范围下可以对本发明技术方案的细节和形式进行修改或替换，均落入本发明的保护范围内。
第一实施例：下面，以来自临床转化研究的高维组学数据的预测模型的构建过程为例，详细介绍本发明方法。
1、方法：
1.1、动机
1.1.1、问题陈述：以数据库中的n个对象(i＝1,2,…,n)为样品。在每个第i个对象(xi)上，将观测到的一组高维(该例中为p维度)稀疏的协变量，表示为xi＝(xi1,xi2,…,xip)，基于hdod的典型特征，其中协变量的数目通常比样品量大很多。在每个第i个对象上还观测到对应的目标yi的结果变量，它可以是二元的、分类的、连续的或截尾的(即，部分被观测到的)。所有观测到的数据的似然可写成
其中上述求和函数中是对n个对象求和(即i＝1到n)，f(yi|xi)是已知协变量xi的yi的条件密度，并且f(xi)是协变量的多变量分布函数。本公式目的是为了通过对f(yi|xi)的边际均值建模，建立结果(yi)与协变量(xi)的关联，其可以被写成
g[e(yi|xi)]＝h(xi,θ)[2]
其中g(.)通常称为关联函数，由每个结果(yi)限定，e(yi|xi)是条件均值，并且h(xi,θ)是由未知参数θ索引的协变量函数。
1.1.2、表现定理：kimeldorf和wahba(1971)已经表明，当协变量函数是未知的并且未被限定，且已知观测的样品为(x1,x2,…,xn)时，上述函数通常可表示为
其中θk是样本特异性的未知参数，k(x,xk)作为核函数，且需要是半正定的。其中，包括一类用于相似性度量的核函数，用于量化x与xk的相似性：当观测到x与xk相同时，所对应的项是θkk(x,xk)＝θk；当x与xk完全不一样，θkk(x,xk)＝0；当x与xk是相同或几乎相同，对应项可以合并为θkk(x,xk)+θk'k(x,xk')≈(θk+θk')k(x,xk)＝αkk(x,xk)。最后，可期望的是，如果第k个个体的协变量特性不与对应的结果相关联，系数θk很可能等于零，这里的系数θk是用于量化结果与第k个个体的相似性度量k(x,xk)的关联。zhu和hastie使用上述观测中的一些情况，通过对一些k(x,xk)项的分组来描述一向量机的输入方法。现今的理论基础和相关研究提出了oor方法，其可表示为
其中sk(xi)＝k(xi,zk)是第i个对象xi与第k个范例zk的相似性度量，q是范例的数量(将在后文描述)，并且(α,βk)是待被估算的未知回归系数。当回归系数βk不等于零时，意味着当所述xi的hdod的特性与zk相似时，xi通过上述oor与结果关联。
oor将结果回归到对象x与范例的相似性，而不是作为协变量回归到hdod。正如预测的那样，本例中的回归系数是针对于与范例的相似性的，此类情况类似于计算机科学家经常使用的数据查询。正如预测的那样，oor是对范例特异性关联的“整体解释”，而不是对协变量特异性关联的“整体解释”。
1.1.3、临床直觉：oor动机来源于临床医生的直觉。临床医生通常收集来自医疗记录、体检以及诊断实验室测试的多方面的信息，这种信息即为一种hdod数据，然后基于这一信息加上他们掌握的过去的案例经验进行临床判断。一个有经验的临床医生会将新患者与先前治疗的患者或教科书或文献中的典型案例作比较，并且通过样本量为1来减少比较的次数，作出合理的临床判断。可见，在本质上，临床医生的过程也是oor过程。
1.2、oor框架
图1提供了oor过程的示意图。作为输入数据的hdod是一个关于多个单一、连续的元素的大型协变量矩阵(图1a)。作为对于任何有意义的聚类分析的常规要求，过滤掉那些是噪声信息或不可能包含信息的协变量是很重要的。当没有结果数据时，oor首先通过无监督聚类分析来对hdod确定范例zk(图1b和c)。无监督学习的结果形成包括q个范例(z1,z2,...,zq)的阵列。基于选定的相似性度量k(xi,zk)(见如下讨论)，可以计算每个第i个对象xi与每个第k个范例zk的相似性度量(图1d)。通过把相似的度量作为协变量，可得到稠密协变量矩阵(图1e)。在广义线性模型下通过适当地选择关联函数，可以再选择包含信息的范例，来形成预测模型(图1f)。在下文中，通过训练集的roc分析，对oor预测模型的敏感度和特异性进行初步估算，然后对验证集进行roc分析。下面的章节集中描述了oor框架的重要组成部分。
1.3、无监督学习
无监督学习方法不参照结果数据，其目的在于探索跨基因和跨对象的hdod协变量的相关结构。单纯从统计框架来说，无监督学习方法可以不参照结果数据，利用对数似然函数的第二部分，基于hdod的相关结构来构建hdod。以往，聚类分析通过相关性组织基因和/或样品，得到的样品集群可以实现识别目标的质心。因为聚类分析，故，质心与该集群内的样本有相对高的相关性(或相似性)，并作为示例来表示多个样本。
当处理hdod时，通常会出现很多远离集群的含有相对独特的hdod特性的对象，此类“独特的对象”可定义为不易由质心或它们的组合来表示其hdod特性的对象。这种独特的对象被作为范例时，可用下面的回归方法来确定这些对象。假设预先已确定了一组起始的t个质心作为范例，表示为[1],[2],...,and[t]，首先，通过下式将所有对象的hdod回归到质心的协变量，而不是回归到那些由hdod代表的集群：
其中xi表示第i个对象，
针对每个个体估算来自上述线性回归的残差平方的总和(srs)，并且计算由这些信息代表的残差变化的分数。当第i个个体满足下列公式时，其可以加入到该组的范例中：
i＝argmax(srsi/srs0),subjectto(srsi/srs0)≥f,[5]
其中srs0是不包含范例的srs，f是一个预先选择的阈值(例如，0.5)。需注意的是，由于该分析未参照结果选择范例，故并不影响任何下游的监督学习(参见下文)。并且，除了从内部推导范例，还可以包括来自外部资源的范例。
1.4、监督学习
在确定范例之后的下一个步骤是，估算这些范例的相似性是否与目标结果以任何方式相关联。这样的关联分析即称为监督学习(监督学习的来由：计算机科学家创造该词是出于对大众的吸引力)。根据不同的结果函数，如二元函数、分类函数、连续函数或断尾函数，监督学习可以使用广义线性模型(glm)来估算与范例相似性结果的关联。在选择变量时，我们建议使用惩罚似然方法控制过拟合问题，特别是使用最不绝对收缩和某选择算法(例如，lasso)来选择包含信息的范例。
下面对单一结果(yi)与q个选择的范例zi＝(zi1,zi2,…,ziq)的相似性度量进行分析，为了估算结果yi与范例zi的关联，已知相似性度量si，可以通过以下分布函数表示yi的分布
其中ωi称为标准参数，主要由以下的边际均值限定，c(yi)是预先限定的函数，δi是用于确保上述函数整合到相同单位的归一化常数。上述功能可以通过下面边际均值来充分限定
通过上述回归方程，便限定了与相似性度量关联的边际均值。一旦嵌入glm，可以援引似然理论的整体含义来支持参数的估算和推断。
当应用glm来选择包含信息的范例时，可预测的是，许多范例可能不与结果相关联。可以使用lasso来选择那些包含信息的范例，lasso可以理解为是惩罚似然估算的一种版本，并在oor模型中采用估算回归系数
其中，对所有随机样本通过常规对数似然函数进行第一求和，对q个回归系数的所有绝对值进行第二求和，λ用于确定那些非零回归系数的惩罚幅度的调谐参数，且，估算调谐参数λ可得到基于交叉验证的最小预测误差。其中，f表示密度函数，yi是对应第i个范例的要回归的结果，si是对应第i个范例的相似性度量值，á是截距，
1.5、相似性度量
对于oor，选择合适的度量以及对象和范例之间相似性的程度来测量相似性是至关重要的，因为它决定了如何计算相似性、如何确定集群、如何确定范例。通常，所述度量的选择取决于hdod的特性和对于结果的解释。本例中，下面介绍了几种常见的相似性度量。按照惯例，该相似性度量是距离的倒数，即1和0的相似性分别等于零距离和无穷大距离。
1.5.1、欧氏距离函数：考虑hdod协变量向量(xi,xi')的两个对象，其中所有参数均为数字型参数。其欧几里得距离可表示为
dii'＝||xi-xi'||,[8]
其中，||·||代表平方差的和的平方根。由平均值和标准差对协变量归一化处理时，欧几里得距离具有相关系数的单调关系，该单调关系通常用来测量相似性。最近，frey和dueck使用了负欧几里得距离的平方，即-||xi-xi'||2，作为相似性度量。
1.5.2、径向基核函数：该核函数是将上述距离转换为相似性度量的另一种常见的函数，如下：
其中
1.5.3、余弦相似函数：在信息理论中，余弦相似性是两个向量(xi,xi')之间的常用度量，其可被写为
其中
除了上述常用距离或相似性的度量，还有其它的域特异性的度量。在遗传学的背景下，遗传学家曾用“血缘同源”，“状态同源”或亲缘关系系数作为遗传相似性的度量。当处理文字时，也存在很多用于语义相似性的度量。可灵活选择最适用于给定的oor中应用的相似性度量。
1.6、oor与协变量特异性回归的比较
如上所述，oor与协变量特异性回归(csr)的方法存在固有关联，但仍存在一些如在分析对象、应用领域、结果解释和分析能力方面的根本差异。在下文中，将着重描述两个回归方法之间的一些关键差异。
csr的常规设计中，可以对结果与一个或多个协变量的关联进行估算。csr的理想特征是，在对回归中的其它协变量进行控制后，可以将协变量特异性的关联分离出来。基于此及其他原因，过去几十年中csr一直是大部分统计应用的“主力”。然而，在“大数据”时代，csr的应用遇到了一些挑战，首先，在处理hdod时，由于p>>n，不能使用csr的一个典型的应用来同时分析所有的协变量。第二，csr假设的前提条件为，协变量对回归模型具有影响。当包括多个相关协变量时，该假设可以使csr的外推结果与很少或甚至没有被观测到的协变量进行关联。当上述假设成立时，则可发挥csr的作用，否则，csr的外推可能会被误导。第三，csr适用于分析数值的hdod协变量，而其在非结构化数据，如基因组序列的应用方面受到限制。
而提出oor的主要目的就是克服上述限制。oor将协变量矩阵(n×p)转换为相似性打分矩阵(n×q)，其中q比样品量n要小得多(图1)。这种转换使得oor能够处理hdod。可见，oor解决了不同于csr的关联问题，其主要目的不是为了确定哪些协变量与结果明显关联，而是确定出哪个患者群体很可能与结果相关联。通过使用相似性度量，oor适用于将结果与任何维度的hdod进行关联。相对于多元“外推”的问题，oor旨在估算涵盖范例的范围内的参数，自然缓解了外推的问题。
2应用：
在tcga肺癌研究方面的应用
与其他癌症相比，男性和女性由肺癌引起的死亡率均最高，占所有癌症死亡率的28％左右。由于被诊断时大都已经为晚期，因此肺癌预后很差。肺癌早期的预后会好些，五年存活率约为60％。即使在i期患者中，一些患者的存活期也都相对较短。预测i期患者的预后存活率成为了研究热点，以便于肿瘤学家为较低存活率的患者可以设计更积极的治疗计划来改善预后。
2.1、数据源
为了解决这个问题，从xena(http://xena.ucsc.edu/)下载了临床表型数据和rna-seq数据。截止2015年6月10日，该网站发布的数据总共包括1299个样本。在对临床表型数据和rna-seq数据关联，并进行基本的质量控制后，对1124个肺癌病例(571个腺癌病例和553例鳞状细胞癌病例)进行了研究，在此说明的是，上述的临床表型数据和基因表达数据都是完整的。将全部数据随机分配到训练集和验证集，以用于下游的分析。图2示出的训练集和验证集的所有患者的诊断年龄分布，表明了两组患者具有类似的年龄分布。对于性别、肿瘤类型和分期进一步的研究发现，训练集和验证集的频率在很大程度上是类似的(见表1)。关于存活率，与四个协变量相关联的估算kaplan-meier曲线在训练集和验证集间也是类似的(见图3)。
表1训练集和验证集中年龄、分期和肿瘤类型的分布(数据集来自tcga)
2.2、预后存活率
在当前组合的数据集中，包括了腺癌和鳞状细胞癌的患者，从数据来看，存活率并不与年龄(p值＝0.143)、或者与性别(p值＝0.605)、或与肿瘤类型(p值＝0.444)显著关联，存活率而与肿瘤分期(p值<0.001)显著关联。首要目标是构建一个预测模型，用于预测i期患者的预后存活率。在训练集中，有296个i期患者。为了保持用于构建预测模型的样本量，并不会按肿瘤类型、性别或年龄将肿瘤样本进行分类，因为这些并不与存活率显著关联。
2.3、rna-seq数据
tcga研究人员使用illumina公司的hiseq2000测序平台，产生mrna的短读取，然后将这些短读取与20531个基因的基因表达水平的估算进行整合(https://support.illumina.com/sequencing/sequencing_instruments/hiseq_2000.html)。在此说明，对基因表达值进行排序后，用相对应的排序代替表达值，虽然这样会失去了rna-seq的定量信息，但是基于排序的转换消除了样品与样品之间的不均一性。
2.4、基因筛选
在进行oor分析之前，先从训练集的20531个基因列表中筛选出基因。为了保持这种经验的特性，将“分期”作为一个关键变量进行分析，这是因为分期与存活率具有明显的关联，从i期到iii期的变化示出了从早期癌症到晚期癌症的进展情况。正如预测的那样，许多基因在癌症进展中被上调或下调。据推测，甚至在早期的癌症，癌症也会出现进展，但它们的形态特征尚不能被观测到。通过将基因表达水平与期(i期vs其它更高的期)相关联，计算每个基因的z分数和相关联的p值(图4)。使用p值＝0.05的阈值(选此阈值，是考虑到达到传统显著水平的所有单个基因)，可选择出831个基因。在去除一些高度关联的基因后，最终可得到789个基因的列表，并作为oor分析的输入数据。
2.5、基因图谱的选择
使用欧几里得距离和r包'gplots'的heatmap.2函数的完全连锁选择方法(https://www.r-project.org/)，对输入数据进行双向聚类(图5)分析。基因(列)树状图表示存在多组共变基因。其中，采用垂直块(白线)表示一组高度共变的基因。由于该步骤目的在于确定范例，故，用样品(行)的等级聚类来表示存在的多个组，其中对7大聚类进行了突出显示，由六条黄色线区分。由于视觉图谱具有较强的定性化的印象，可直观的显示出多组对象具有的不同的基因表达特性。虽然将数据图像化有其优点，但对于使用合成数据来生成可重复的结果，也存在着挑战。首先，视觉图谱的感知因人而异。第二，呈现的视觉图谱取决于所选择的可视化参数，如所选颜色、所选颜色深浅等。第三，在视觉上几乎无法区别出系统图谱和随机图谱。实际上，已经随机选择了1000个基因进行聚类分析(未示出)，并进行了“模拟实验”。在这些实验中，可以偶尔看到一些由双向聚类所产生的图谱。总的来说，该图谱与通常得到的图谱(图5)区别不太大。
2.6、路径分析
除了采用图谱所提供的视觉印象，还可以想到的是，选择的基因包括生物学上有用的元素，该基因的选择是基于i期的关键指示器所选择的。当然，由于p值0.05是随机选择的，某些基因估计是纯粹被随机选中的。如果要分析被选中的这789个基因的生物学意义，可以采用一种关于路径分析的网络工具targetmine，(http://targetmine.mizuguchilab.org/targetmine/begin.do)进行路径分析，所分析出的10个路径包括对应于基因富集p值小于5％的基因(补充表s1)。表2的第一块区域列出了这些路径，包括细胞周期、有丝分裂的细胞周期、m期和减数分裂重组，所有这些都与癌细胞从i期到更高期的加速细胞生长相一致。更需注意是，除了输卵管上皮细胞，所有相关的组织似乎与气道的上皮细胞相关联(表2)。如下补充表示出了各种组织的基因列表(表s1)。
表s1各种路径和组织的基因列表
表2：根据基因与分期的关联分析选择的一组789个基因的确定的路径和组织(i期与其他更高期相比)
2.7、范例的研究
如上确定范例的步骤完成后，可得到代表观测到的集群或单体的273个范例。考虑到大多数范例都不太可能与预后存活率关联，因此进行边际关联处理，仅保留那些有边际关联的范例。通过关联p值为0.05选出22个范例(该22个范例将被lasso进一步进行选择)。表4列出了边际关联分析的估算系数、风险比、标准误差和p值。
表4：对于训练集、验证集和所有数据的风险评分，并进行cox回归分析，得到的估算的系数、风险比、标准误差、z分数和p值
由所选的22个范例，计算每个对象与每个范例的相似矩阵，生成“稠密协变量矩阵”，可参见图1e所示。图6示出了具有296行×22列的相似性矩阵。灰色、黄色和红色分别对应于对象与范例的弱、中等和强烈的相似性。通过聚类分析将296个对象和22个范例构建为不同的子集。将22个范例分成三组，其中“列”中的彩条代表每个范例的边际关联：红色为保护关联，绿色为风险关联。为了进一步深入了解预后存活率，此处创建了一个为期一年的存活率指标，该指标需要取存活的值(0和绿色)、死亡值(1和红色)和断尾值(丢失和黑色)，并将彩条放入“行”中。为了方便观测，图中突出显示了两个高风险区，所述高风险区包括与拥有风险表达特性的范例高度相似的多个对象。与此同时，图中还突出显示了具有更好的一年存活率的对象。
2.8、根据所选范例构建预测模型
根据所选范例，先通过lasso从“稠密协变量矩阵”选择包含信息的范例。结果示于表3的最后一列，其中的11个范例被选为用于预后结果的包含信息的范例(图7所示)。表中已列出估算的回归系数，而未选择的范例其系数则设置为零。可观测到的是，在第8列中的估算的回归系数往往比其对应的第三列的来自边际回归分析的系数小，这可能反映了lasso将边际关联分配给了与多个范例的关联，而惩罚一些例如第一范例那样的不稳定的范例(表示为例1)。
表3：在通过cox回归模型进行的单变量oor分析和cox回归模型下的基于lasso的多变量oor分析中，与每个范例相关联的估算回归系数(估算系数、风险比、标准差、z分数和p值)。
根据估算回归系数，可以通过下式继续计算任一训练集、验证集或所有i期样本中的每个个体的风险评分
其中，
2.9、惩罚参数的蒙特卡洛稳定性分析
当应用oor时，lasso要求必须估算惩罚参数(λ)，这个参数的选择对变量的选择产生深远的影响。在真实值未知的情况下，常见的方法是使用交叉验证法来估算这个惩罚。不幸的是，交叉验证产生随机估算惩罚参数所带来的问题是“估算惩罚参数有多稳定？”。出于这个目的，重复进行了1000次蒙特卡罗模拟实验，在每次实验时，使用r'glmnet'包的'cv.glmnet'函数(https://cran.r-project.org/web/packages/glmnet/index.html)通过10倍交叉验证来估算惩罚参数。图8示出了用对数比例来表示的估算惩罚参数值的经验分布情况。可见，共有20个不同的惩罚值，范围从0.016到0.091。惩罚值越小，表示越多范例被选中。在当前的应用中，0.091的惩罚值对应没有选定范例的空模型，可参见图7的惩罚值(λ＝0.021)项。通过回顾性的分析，(λ＝0.021)与模型(λ＝0.030)相比是一个稍小的惩罚值。
2.10、通过自助分析选择范例的稳定性
考虑到惩罚参数值的范围，预计所选定范例是可变的。为了估算选定范例的稳定性，重复进行了1000次的自助分析。对每一个自助样本，对其观测到的基因表达值和对应的存活结果随机采样，然后放回该样本，以确保与训练集的分析数据集具有相同的样本量。基于20个固定的惩罚值，通过lasso从同一个分析数据集中选择范例。表5列出了计算的kappa值，kappa值用来表示估算的选定范例与不同惩罚参数值的一致性，其中kappa值的范围从0(无一致性)到1(完全一致)。对应于上述1000次的自助分析，表5中上三角形的参数为kappa平均值，表5中下三角形的参数为估算的平均kappa值的标准差。可见，相邻惩罚值的一致性接近1，该一致性随着相应的惩罚值的发散而减小。为了保证定量水平的一致性，根据上述1000次的重复分析结果，计算与所有22个范例相关联的系数的平均估算值。然后，以成对的xy图(图9)进行显示。同样可见，一致性在定性和定量估算之间基本一致。其中，右上角的xy图，除了有两个极端的惩罚值，大多数范例的平均系数保持一致。
表5：通过lasso选择的范例之间的kappa平均值，其中右上方的三角区内为不同的惩罚值，下面的三角区为1000个自助样本的标准偏差。
2.11、与协变量特异性回归分析的比较
如上所述，也可以使用csr方法，并结合lasso构建预测模型。为了进行比较，选择同一组的所筛选出的789个基因，对存活结果使用cox回归模型、并应用lasso选择预测器，可计算得到6个预测器。然后，对训练集和验证集的每一个对象进行预测值的计算，该预测值类似oor的风险评分。将这些预测值与从oor获得的预测值绘制成如图10所示。其中，来自csr和oor的两组预测值在训练集和验证集中均彼此相关联(r2＝0.318)。其主要的区别在于，来自范例特异性的预测模型的预测值比来自协变量特异性的预测模型的那些预测值具有更广泛的范围。换句话说，所述oor预测模型比csr预测模型具有更好的分化特性。
3、分析
对于数据科学家来说，无论他们的学术根基在生物医学信息学、计算机科学或生物统计学，在转化生物医学研究中越来越多地使用组学技术是一个前所未有的挑战。来自转化研究的hdod都有一个共同的特征，即样本量相对较小，但协变量维度却非常高。为了应对这一挑战，引入了面向对象的回归(oor)方法，其类似无监督学习方法和监督学习方法的结合。oor关键点在于对范例的确定，该范例可理解为：由他们的hodo特性代表的多个集群对象，或者具有相对独特的hodo特性的对象。参考这些范例，oor估算每一个对象与范例的相似性，并侧重于确定出包含信息的范例(即与感兴趣的结果相关联的范例)的特异性的相似性。除了探索范例与结果的边际关联，oor也可用于选择包含信息的范例，并形成预测模型。相比传统的协变量特异性预测模型，范例特异性预测模型看起来具有更大范围的预测值(图10)。
为了说明oor，将其以及从tcga获得的基因表达数据一起应用于肺癌研究，建立一预测模型，该模型用于分类已被诊断为i期肺癌但具有明显不同存活时间的患者(无论是腺癌或鳞状细胞癌)。首先确定来自训练集的11个范例，并生成作为相似性测量的加权的和的风险评分，该风险评分可显著地预测验证集的存活率(p值＝0.0145)。根据假设的预测模型，对验证集的所有对象的预测风险评分进行计算，其分布可参见图11左侧区域所示。在右侧区域中示出的与风险评分1、2、3、4相关联的估算的存活曲线表明，随着风险评分的增大，存活率逐步恶化。oor虽有很多优点，也有一潜在的弱点：用于衡量相似性的度量的选择是有点武断的。在关于聚类分析或无监督学习的文献中，使用了多种相似性度量，不同的相似性度量根据不同的应用环境各有优缺点。从这个角度看，oor提供的相似性度量的选择具有一定的灵活度，适合于上述的应用。
oor在概念上与其他分析方法相关联，k近邻方法(knn)即为相关联的一种方法，knn被广泛用于计算机科学文献的数据挖掘，其核心思想是，由某些特性定义的相对“亲密邻居”的对象往往有类似的结果。所述k近邻方法可以不用于做任何建模假设，而用于构建预测模型，因此也被称为非参数预测模型。但是，k近邻方法并没有考虑到的一个事实是：许多邻居具有同等的结果关联(无论是无效假设或备择假设)。在这方面，oor可以被认为是k近邻方法的延伸或回归函数最近邻的估算。
另一种密切相关的方法是成员等级分析法，简称为gom。从概念上讲，gom通过引入一组潜在成员变量，假设该潜在成员变量的分布是合理的，gom对结果的联合分布以及协变量建模，在整合了所有gom潜在成员变量后可计算似然。gom参数可以理解为是与个体相关联的属性，而不是单一协变量的特定边际。尽管gom和oor提取关于个体或对象的属性信息，拥有相同的概念目标，但是oor侧重于观测到的结果和协变量的经验观测，而无需调用任何潜在随机变量。
oor使用的相似性度量的概念也与统计遗传学中开发并使用的多种方法相关。虽然本文并不意图追踪这些联系，但需注意，经典和现代遗传学旨在发现结果相关联的易感基因，往往会导致相关个体中的相似性，所述相关个体中的相似性比无关个体具有更多遗传变异体。在遗传学的早期，隔离和联系方法用于描述和发现家族聚集性基因。在现代遗传学，一些研究小组提出，通过评估遗传标记的相似性并使用相似性回归来发现疾病基因。虽然有着相似的科学目标，但oor使用相似性评分来发现哪些范例有更高的疾病风险，而不是发现哪些snp(单核苷酸多态性)与疾病相关联。
oor与最近流行的被称为序列核关联测试(skat)方法，也存在着内在的联系，这是因为oor和skat都采用表现定理作为理论基础。在很大程度上，skat使用表现定理来表示所有snp的组合和它们的作用，作出关于所有回归系数的合理多变量假设，并根据无效假设测试他们的偏离。最近，pan(2011)表明，skat测试本质上等同于上述oor提到的相似性回归。而oor比skat更进一步，其将结果回归到相似性评分，而不是假设它们为随机变量。
上述已经介绍了用于分析hdod的新的分析框架。介绍了上述技术推导，以及与现有方法的各种关联，oor给我们介绍了探索hdod的“整体关系”与临床结果的分析框架。协变量特异性研究已经应用于“简化论观点”几十年，上述方法是对协变量特异性研究的补充。在大数据和系统生物学的时代背景下，该整体的框架不仅会促进hdod的系统研究，也会生成组学数据的“可重复结果”。
第二实施例：下面，以构建多态性多等位基因hla基因的疾病预测模型为例，进一步详细介绍本发明的方法。
1、方法
1.1、动机
对从高维多态性基因研究产生的协变量数据进行分析。具体而言，包括将t1d和八个ii类hla基因(hla*drb1,*drb3,*drb4,*drb5,*dqa1,*dqb1,*dpa1,*dpb1)(手稿:zhaoetal2015,待提交)的病例进行对照研究。由于它们结构的多态性，在任何单条染色体中只会出现hla*drb3,*drb4和*drb5等位基因其中的一个，因此，以下用hla*drb345表示所有这三个基因的基因型。其中，每个基因包含两个等位基因，各等位基因代表一个完全分相核苷酸序列。当第j个基因具有mj个可能的序列变异时，如果一对等位基因处于哈迪—温伯格平衡(hwe，即统计上独立)，该对等位基因的基因型可以具有mj(mj+1)/2个可能的基因型多态性之一。在多个基因位点的基因型阵列被称为基因型分布。如果这些基因处于连锁平衡(le，即统计上独立)时，基因型分布的总数在理论上是其交叉乘积
将所观测对象的基因型分布作为一个整体是上述方法的一种补充，并通过系统方法或整体方法研究它们与结果的关联，即，将所观测的基因型分布作为范例，计算其他基因型分布与该范例的相似性，并评估与范例的相似性是否与疾病表型相关联。已知种群研究的样本量为n，从研究内部可能获得的范例总数最大为n，这样既减少了数据维度，又没有任何信息的缺失。如上所述，这些由八个ii类hla基因形成独特的基因型分布的实际数目实际上小于样本量n。如果将所有观测到的基因型分布作为范例，则可以直接评估所有这些范例的对象的相似性度量与t1d的关联。这些范例观测的规范化也促进了oor的提出。
从形式上看，对基因型分布表示为
其中，logit是对疾病的概率
1.2、面向对象的回归框架概述
oor的动机是直接的，而它的表现也非常简单。目前，要使用oor必须解决不同的三个方法学问题：1)相似性度量的选择，2)范例的选择，3)包含信息的范例的选择(即非零βk系数)，各种不同的选择会导致生成不同版本的oor框架。
1.2.1、相似性度量：单纯从理论上考虑，相似性度量的选择需要确保核函数是对称和半正定的。在实践中，大多数的相似性度量都比较合适，且与应用的场景相关。在此，使用适合于遗传分析的相似性度量。假设
其中i(.)是一个指示函数，并且每个
1.2.2、范例的选取：
范例的选取方法有很多，主要取决于所要分析的目标。第一种方法，鉴于这些hla基因可能具有不同的基因型分布，可通过对所有对象的聚类分析，以及采用特定的相似性度量来确定主要图谱。其中，可将每个集群内具有代表性的基因型分布选为范例。第二种方法，将每个独特的基因型分布选为范例。第三种方法，从文献中确定一组基因型分布，这样也可以确保结果可被合理解释。第四种方法，通过对某些联合关联或基因-基因的相互作用的研究，设计一定的基因型分布来作为范例。关于范例的选取方法，后文将有描述。
1.2.3变量选择：
除了降维，预测器的数量可以与oor的样本量n一样大。正如预测的那样，由于这些范例不与疾病表型相关联，许多回归系数等于零。因此，在oor中的一项重要任务即为选择包含信息的范例。目前，在所有用于变量选择的技术中，惩罚似然法是应用最为广泛的。在此，相对比传统的选择变量的逐步回归方法，本文考虑三种惩罚似然法：lasso、脊回归和弹性网络。
1.3、范例的选择
如前所述，oor方程确定的范例
1.3.1、聚类分析：作为编码人体先天免疫的必需基因，hla基因在整个人类的进化过程是高度选择的。如前所述，hla的基因型数据倾向于集群，这可以通过成对出现的相似性度量，在n×n相似性矩阵
1.3.2、“独特”的对象：不进行任何聚类分析，而对成对的相似性度量
1.4、变量选择
在确定范例
1.4.1、逐步选择：最有名的传统变量选择的策略大概是由预测器进行的逐步选择，无论是仅向前，仅向后或双向，均是基于信息准则(ic)的度量的，基于ic的度量可如akaike’sic(aic)或者bayesianic(bic)。基于大量文献对似然估算的描述，须注意的是，概率模型可以构建如下带有aic惩罚的对数似然函数：
其中，k′i为与q个范例的相似性度量的向量。基于aic的连续值，逐步选择过程采用了双向策略，选择那些与疾病结果(yi)相关联的预测器。其存在的一个问题是，当预测的数量接近样本量时，这种策略往往会失败。存在的另一个问题是，可能由于过拟合，独立的数据集的预测的可重复性往往比较差。
1.4.2、惩罚似然：当范例的数量接近样本量，首选的变量选择的方法是使用惩罚似然法，所述惩罚似然法包括三种被广泛使用的方法：lasso、脊回归和弹性网络。使用上述公式[13]中相同的符号和变量来表示，该惩罚对数似然函数可被写为
其中λ是用以确定惩罚水平的调谐参数，|β|1和|β|2分别是l1的范数和l2的范数，θ分别取值为0或1或0.5，分别对应lasso、脊回归和弹性网络。优选的，估算的调谐参数λ具有基于交叉验证的最小预测误差。
1.5、惩罚参数和变量部分
众所周知，在惩罚似然方法的文献中，调谐参数将估算回归系数的偏差与他们的估算的方差进行交换。通常，惩罚参数的估算是通过交叉验证进行的，然而，交叉验证过程是一个随机过程，并且因此估算的惩罚参数也是随机的，因此会不可避免地影响变量的选择。在这里，建议采用多次重复交叉验证过程，并基于随后会利用固定的惩罚参数进一步对变量选择的稳定性(参见下文)进行评估，估算它的经验分布。计算上，可用10倍交叉验证估算惩罚参数(在cv.glmnet默认推荐，glmnet的r实现)，并重复计算，比如100次。所有经验估算的参数随后被用于构建经验分布，以评估这些估算是否来自单一模式分布。
1.6、评估固定惩罚参数的变量选择的稳定性(λ)
实际中，所有处理复杂或高维数据的变量选择方法面临的主要挑战，是选择的变量的稳定性。oor的变量选择也不例外。在评估上述的惩罚参数估算的经验分布时，要关注选定的包含信息的范例是否稳定。为了解决这个问题，可使用自助法。简要地说，从研究群体随机抽取样本观测并放回，这样可以保持样本量不变。对于每个自助样本，进行具有两个或多个固定惩罚参数和/或使用不同的方法的惩罚似然分析。然后，计算kappa统计，度量变量是否由两种或更多的方法一致地选择。
1.7、t1d病例对照研究
正如上面提到的，青少年i型糖尿病(t1d)和hla基因的病例对照研究促进了oor研究的发展，其中的细节已被公开(zhaoetal.2015提交)。简单地说，这项研究确定了970个i型糖尿病患者作为病例，他们的年龄范围从1岁到18岁，且来自不同位置的诊所。并从相应的地区确定了448个未患i型糖尿病者作为对照。遵循人类受试者的审查和批准的要求，从所有研究对象中采集血液样本，并提取他们的dna。虽然测试多个分子靶点，本研究使用下一代测序技术以评估hla基因的高分辨率基因型(hladrb1*，*drb345，*dqa1，dqb1*，dpa1和dpb1)。这项研究的分析目标是研究i型糖尿病与hla基因的关联，并构建i型糖尿病特性与这些hla基因型的预测模型。为了建立验证集，随机选择了479个病例和226个对照作为训练集，其余部分作为验证集(222个对照和483个病例)。对照以及病例的所有基因的等位基因频率在训练集和验证集中很大程度是类似的(为了说明，补充表s2包括hla-drb1对照以及来自训练集和验证集的病例的等位基因频率)。
表s2在训练集和验证集中对照以及病例的hla-drb1等位基因频率
2、结果
2.1、应用于hla-drb1
为了对oor处理复杂的hla数据的过程进行说明，首先对t1d仅与hla-drb1基因的关联进行分析。表6的对角线的上方和下方分别列出了对照和病例中的hla-drb1的基因型分布。对于那些对角线以下的纯合基因型，对照和病例中的基因型频率分别用分子和分母(#/#)表示。该基因型频率表示出的直观印象是，只有44个等位基因的基因型分布是稀疏的，且只有159个独特的基因型，数量上比理论上根据hwe计算的可能的基因型数目990(＝44×45/2)要小得多。其次，需注意的是，某些基因型在病例和对照之间呈现出明显不同的频率，该频率意味着它们与t1d的关联情况。例如，纯合体04：01：01/04：01：01在病例和对照中分别具有0.6/9.3的频率，这意味着15.5的频率比。在另一个极端，杂合子15：01：01/07：01：01在病例和对照中分别具有0/3.4的频率，这意味着这个杂合子看起来可预防i型糖尿病。对于那些常见的基因型，基于当前的样本量对t1d关联的直接评价是实际可行的，且在文献中已经被研究。然而，对于许多不太常见的基因型，因为稀疏、样本量小，以及大量的比较，则很难进行严谨的评估。考虑到期望整体检测t1d与基因的关联，也在寻找可替代的其他分析方法。
表6：对照(对角线的上方)和病例(对角线的下方)的训练集中的hla-drb1的估算的基因型频率。对照和病例的纯合基因型的基因型频率分别对应分子/分母。
考虑通过公式[12]训练t1d与hla-drb1的关联的oor模型，而无需采用任何假设。由于某些等位基因的等位基因频率不同以及与hwe的偏差，理论上可能并不存在许多基因型，即，它们的频率为0(表6)，故，oor则可被简化为
其中fk对应第k个独特基因型的频率，
在本例中，在各对对象之间，相似性矩阵的元素采用值1表示为相同，采用值0.5表示为共享一个等位基因，采用值0表示不共享等位基因。图12示出了其中的705个对象的相似性矩阵的热图，其中示出了共享两个等位基因的对象(红)，共享一个等位基因的对象(黑色)和不共享等位基因的对象(绿色)。从hla-drb1的角度来看，可以识别出一组相同的对象(红色正方形落在对角线上)，以及另一组只共享一个等位基因的对象(绿色长方形)。
为了进一步深入了解范例特异性的边际关联，借由上述oor公式，还可以对t1d与每一个范例的相似性度量进行单变量关联。单变量分析的结果包括了补充表中列出的所估算的对数几率比、标准偏差、z分数和p值(表s2)，以及范例和相关联的基因型。为了更直观地分析，表7中呈现矩阵形式中的四舍五入为整数的z分数，且为简单起见，对应于0.05或更好的显著性水平(没有校正多重比较)，该z分数的绝对值设置为大于等于2。这些单变量分析的结果显示了：hla-drb1*03：01：01和*04：01：01与t1d正相关，其着色为红色条纹。另一方面，6个等位基因hla-drb1*07：01：01，*11：01：01，*11：04：01，12:01:01,13:01:01和15：01：01：01可预防t1d，着色为绿色条纹。要注意的是风险和保护等位基因的杂合基因型倾向于与t1d正相关。
表s3：范例特异性边际回归分析得到的估算的回归系数、标准偏差、z分数和p值。
表7通过oor从边际关联分析提取估算的z分数(四舍五入到整数，等于或大于2)。两个主要的等位基因(hla-drb1*03：01：01和*04：01：01)用于评估较大的风险关联(红色条)。6个等位基因(hla-drb1*07：01：01，*11：01：01，*11：01：01，*11：04：01，*12：01：01，*13：01：01和*15：01：01)用于评估与i型糖尿病的较大的保护关联。
在排除与i型糖尿病没有关联的范例
相比之下，脊回归方法生成所有范例的估算系数，且对任何范例都不取消选择。为了说明，在表s4的所有范例由相应的系数进行了排序。不同于lasso估算方法，脊回归的估算系数取零附近较小的值。其中，估算系数的方向性往往是与病例/对照源的所有范例相一致的。此外，对于那些由lasso选择的范例，脊估算在方向性上与那些通过lasso获得结果也是一致的。表s4的第三列示出了由弹性网络估算的系数，其中选择了39范例。选择的这39个范例大多数与lasso选择的范例重叠。从数量上看，弹性网络和lasso之间的估算系数是高度相关的(未显示)。而逐步回归方法选择了14个范例，其中10个与lasso选择的范例重合。尽管这看似有很高的一致性，但与lasso获得的范例所对应的系数相比，许多估算系数的值往往相当大。
表s4：使用4个不同选择方法得到的训练集中hla-drb1的估算回归系数
为了对通过这四个方法选定的范例的预测模型的性能加以了解，对接受者操作曲线(roc)进行了分析，并评估所有四个预测模型的敏感度、特异性和曲线下面积(auc)。图13示出了在训练集以及在验证集的roc曲线和相关联的auc值，具体包括lasso(图13a)、脊回归(图13b)、弹性网络(图13c)和逐步(图13d)。在训练集中，估算的roc曲线和auc的值约为0.9，上述四个方法在很大程度上都相类似。如预测的那样，在验证集，估算的auc值小幅减少至0.866。其中，前述三种方法的auc值的方差都小于0.001。前述三个方法的roc分析结果的类似，表明可能有许多具有不同的范例以及类似的预测性能的预测模型。
然而，与第四种逐步回归分析的结果对比，预计auc下降到0.5，为空值。这一结果表明，逐步的过程可能会因估算出较大的回归系数，过度拟合训练数据集。
2.2、应用于所有ii类hla基因
为了建立一个i型糖尿病的预测模型，将oor应用到所有8个ii类hla基因(hla-drb1，drb345，dqa1，dqb1，dpa1和dpb1)，使用相同的训练集研究范例，并建立预测模型，并验证验证集中的预测模型。相对于上述相似性度量，此处使用了等式中定义的未加权相似性度量，表示为
基于该相似性矩阵，将被观测的基因型分布的一个子集选为范例。鉴于样本量相对有限和基因型分布的神秘，将训练集中所有独特的基因型分布选为范例。换句话说，选择的所有范例均是独特的，并涵盖训练集中观测到的所有基因型分布。操作上，用于进行成对相似性度量的阈值设置为1，训练集中共有499个范例，且作为描述性关联分析的一部分，应用oor进行i型糖尿病与所有范例的单变量关联分析；并沿hla基因型列出(表s5)了估算系数、标准误差、z分数和它们的p值。其中，范例由z分数排序，并且z分数值与病例和对照状态相一致。
表s5：对hla-dr，-dq和-dp基因的全部499个范例进行边际关联分析，得到的估算回归系数、标准偏差、z分数和p值
目前的任务是要用lasso建立预测模型。在前面的讨论中，逐步方法适用于过度拟合预测模型，而并不适用建立预测模型。即使预测性的auc是所需要的，脊回归往往为所有范例提供“谨慎估算系数”，且对任何范例都不取消选择。而弹性网络相对于脊回归和lasso是一种折衷的方法，其具有与lasso相当的性能。为了对变量选择进行分析，选择lasso建立一个i型糖尿病的预测模型。表3列出了基于lasso的回归系数估算，其中该回归系数估算由回归系数排序。通过lasso选择的共有26个包含信息的范例。通过合并病例(d)/对照(n)和研究识别号码得到范例识别号码。显然，对那些从病例衍生的范例的估算系数倾向于为正，而对那些来自于对照的范例的估算系数倾向于为负。例如，与范例如d1612高度相似的对象，具有相对高的t1d的风险；与范例如n000982相似的对象，将有相对较低的t1d风险。
根据从训练集估算的作为权重的系数，可构造一个风险评分作为加权的和的公式，如下所示：
其中，对那些所有26个选定的范例求和，在表3中示出估计风险评分
为了验证上述预测模型，采用了固定的范例和上述模型中的相关联的加权，计算验证集所有样本的风险评分。参见箱图的左下图，示出了对照和病例中的风险评分的分布(图15)。显然，在验证集的风险评分的经验分布与在训练集中的风险评分的经验分布在很大程度上是类似的。此外，验证集的roc分析显示了相类似的灵敏度特异性曲线，该roc曲线的auc＝0.89(图15)。
2.3、选择范例的稳定性
已知的是，该惩罚参数(λ)的选择对变量的选择有直接和深刻的影响。常规的交叉验证通常用于确定出可实现最小偏差的惩罚值(或其他性能度量，如分类误差，或auc)。图16的顶图显示了偏差与不同的惩罚参数值(对数刻度)的xy坐标图。它示出了最低的估算惩罚参数的对数值，所述对数值取值在-6.0到-5.5之间。此函数的平坦性意味着对应于最小偏差的估算惩罚参数在很大程度上受交叉验证过程的影响。为了评估它的影响力，重复1000次估算惩罚参数，并估算相应的值。图16的下图显示估算惩罚参数的经验分布。可见，在训练集中估算的惩罚值是离散的15个不同的值，这可能是因为相似性矩阵的离散性造成的。
由于惩罚参数的值会影响变量的选择，需关注的是，所选择的变量在不同的惩罚参数值下是否是稳定的，其次，即使具有固定的惩罚参数，“选择”本身是否稳定。为了解决这个问题，对15个不同的惩罚参数值进行了自助分析。对于1000个自助分析样本的每个样本，分别设置固定的
如前所述，有多个类似性能的预测模型。现在的问题是，即使选择的范例以及相关联的系数不同，预测模型在惩罚参数值不同时是否也有类似的性能。为此，使用lasso，在固定的惩罚参数值下，选择包含信息的范例构建相应的预测模型。对每个预测模型，进行训练集以及验证集的roc分析。图18示出估算auc值的15个roc分析结果。显然，roc曲线基本上是类似的。在训练集中auc值从0.91变化到0.93，而在验证集中这些值约为0.89。
鉴于类似的性能和不同的惩罚参数值下选择的范例的高度一致性，选择了中等惩罚参数值
3、分析：
在本文中，描述了一种面向对象的回归(oor)的新方法，来建立关于生物大数据的共同特征，即高度多态性基因的预测模型。为了解决多态性基因的复杂性，首先，通过oor确定一组范例，其中，该范例的基因型分布在所观测到的基因型中具有代表性。然后，通过oor选择每个对象和范例之间的取决于场景的基因的相似性度量，作为一个新的“度量”来度量所有对象和范例的相似性，并创建协变量矩阵。然后，通过采用现代惩罚似然方法，通过oor选择一组包含信息的范例来构建预测模型。然后，作为“经典”的回归方法，使用oor分析“范例”与疾病的单变量关联以及多变量的关联。不同于常规侧重于单个基因的回归，oor的回归系数在量化疾病与范例相似性的关联时，需要结合上述新的度量来进行分析，即结合上述与范例的相似性来确定风险等级(见下文关于整体评估的详细讨论)。从这个角度来看，oor是对常规的回归方法的一种补充。
在给定的说明性的例子中，使用了在瑞典进行的一项病例-对照的i型糖尿病研究，探讨了疾病与hla基因的关联。简单地说，为了说明oor及其解释，初步研究的重点侧重于i型糖尿病与单个基因hla-drb1的关联上。作为以基因为中心的回归的补充方法，oor的单变量分析揭示了t1d与单个携带hla-drb1*03:01:01,*04:01:01,*07:01:01,*11:01:01,*11:04:01,*12:01:01,*13:01:01和*15:01:01基因的关联图谱。在该例中，使用了hla-drb1建立了i型糖尿病预测模型，其中通过四个不同的变量的方法来选择包含信息的范例。在针对上述例子的应用中，lasso选择了23个包含信息的范例，与风险升高相关联的范例趋向于来自于病例，而那些与风险降低相关联的范例往往是来自对照。另外，还发现，通过逐步方法选择的范例与那些由lasso选择的范例趋于重叠，但是相关联的系数估算的绝对值往往更大。令人惊奇的是，在评估预测模型的性能时还发现，auc会下降到0.5，表明预测模型完全失效。这大概与预测模型对训练数据过拟合相关。同时，脊回归方法会保留所有范例，并产生与所有范例相关联的谨慎回归系数。而弹性网络方法比lasso会选择更多包含信息的范例，但比脊回归方法少，相当于是这两种方法之间的妥协。另外可见，预测模型的性能与由三种方法选择的范例性能在很大程度上是类似的。为了在简约、诠释和性能之间保持平衡，在本本发明所述应用中选择使用lasso方法。
基于hla-drb1的初步研究的结果，对所有hla基因(drb1，drb345，dqa1，dqb1，dpa1和dpb1)建立了一个预测模型，随后评估其性能，以及评估在不同惩罚参数值下所选择的预测器的稳定性。在训练集中，oor选择了26个包含信息的范例作为预测器，该预测模型拥有极好的敏感度和特异性特性，对应的auc为0.93。固定范例和回归系数后，将预测模型应用在独立选择的验证集上，通过roc分析显示与那些训练集中类似的灵敏度和特异性，此时auc为0.89。如果由外部的数据集进一步验证后，这个预测模型可随时用于在一般人群中筛查t1d。
虽然oor有上述优势，但其也有局限性。通过构建这个模型，oor将基因为中心的回归问题转化为“对象与范例的基因型分布的相似性”的问题。因此，结果的解释取决于相似性度量。例如，如果建立与范例的正关联，结果意味着，任何人只要其基因型分布与范例相似，则处于疾病的高风险。因此，这样的正关联不能精确定位与疾病阶段相关联的特定的基因多态性或其组合。毕竟，oor不会解释哪些基因是重要的。当然，传统的回归方法更加合适用于精确定位病因基因。
另一个问题是与相似性度量的选择相关。在疾病与hla基因的关联分析的场景下，将等位基因身份计数的未加权平均值作为相似性度量。虽然这样是直观的，但可以考虑使用其他度量替代，如使用基于hla基因型血缘同源性推断得到的、对象间血缘同源性的加权平均值。正如预测的那样，相似性度量的选择影响了对结果的解释，实现了场景特异性的灵活性。
oor的特性之一是，当相似性度量构建后，oor需要确定“范例”作为用于构建模型的预测器。oor默认假定范例的数目比样品量小(q<n)。通常情况下，选择范例代表一组具有取决于相似性度量的基因型分布的一个或多个对象。对于8个hla基因，训练数据集的705个对象中有一些是相同的，但许多在各自集群内彼此相似(参见图14为例)。在上述例子的应用中，用于相似性度量的阈值选择1.0时，会从705个对象中选出499个范例。在不缺失信息的情况下，q个范例包含这些复杂的基因型的所有统计信息。假设样本量增加的速度比范例的数目增加的速度更快，仍然可以依靠常规的渐进解法进行统计评估。需要注意的是，范例特异性的预测器彼此间高度相关，例如，范例的相关矩阵。在实践中，用于相似性度量的阈值可以选择低于1.0，用于确定较小的一组范例，以进一步进行分析，这尤为符合较大样本量的需求。
oor的另一个重要特性是，oor结果对于等位基因特异性或基因型特异性的传统回归分析的结果是互补的。hla基因的基因型特异性回归分析，通常仅限于那些常见的基因型，诸如hla-drb1*03:01:01/03:01:01或*04:01:01/04:01:01，其中为了统计分析，还要求观测数量足够大。为了克服此限制，等位基因特异性回归分析假设了模型的额外效果，并量化疾病与个别等位基因的关联。但是，额外效果的假设可能不适合某些等位基因。当然，等位基因特异性回归分析(当包括多个基因时，等同于单倍型特异性回归分析)对于不常见等位基因也同样存在着挑战性。与此相反，oor则绕过上述限制，将分析目标侧重于评估疾病与对象和范例的基因型相似性的关联。
对于结果的解释，oor和协变量特异性回归方法有一定的不同。协变量特异性回归侧重于个别协变量的特定影响，以及如果统计学上显著，对于个别回归系数的解释为相应的协变量有显著的关联，即“简约”论。与此相反，oor评估疾病与对象和范例组的相似性的关联，如果发现一个或多个回归系数从零显著偏离，其结果意味着，与该范例的相似性指示了较高或较低的疾病的风险，即个人风险的“整体”论。事实上，正是oor的这种“整体性”，规避了传统回归分析的复杂性的问题。
oor使用的相似性度量的概念也与统计遗传学中开发并使用的多种方法相关。虽然本文并不意图追踪这些联系，但需注意，古典和现代遗传学旨在通过利用家族内对象的相关性发现结果相关联的易感基因，因为共享的疾病基因在被发现之前，可能会导致相关个体中相似性的增加。在遗传学的早期，隔离和联系方法用于描述和发现家族聚集性基因。在现代遗传学，一些研究小组提出，通过评估遗传标记的相似性并使用相似性回归来发现疾病基因。虽然有着相似的科学目标，但oor使用相似性评分来发现哪些范例有更高的疾病风险，而不是发现哪些snp(单核苷酸多态性)与疾病相关联。
oor还与一些现有的分析方法存在联系。在统计遗传学文献的背景下，oor与序列核关联测试(skat)共享相同的理论基础，即表现定理。最近开发的用于检测gwas基因-基因之间作用的方法中，skat在遗传分析方面受到巨大好评，因为它使用该定理来非参数化地表示snp的所有基因间作用的综合影响，并检测基因-基因之间的作用的存在，这是一个gwas遗传分析的挑战性的问题。最近，pan(2011)表明，skat测试与相似性回归方法本质是等同的。除了共享相同理论基础，oor还具有完全不同的分析目标，即评估疾病与“范例特异性相似性”的关联，并因此直接对范例的相似性度量建模，而不是为范例特异性系数假设一个随机分量。
在更广阔的背景下，oor与核逻辑回归和支持向量机密切相关。所有三种方法共享相同的表现定理，利用该定理的一般表达形式，统称为核机器。然而，oor通过相似性度量将核函数的选择形式化，利用集群战略确定范例，并通过惩罚似然方法选择那些包含信息的范例。可见，基于前人研究所取得的成就之上，oor提供了一种新的方法来分析疾病与复杂协变量的关联。
对于计算机科学文献的数据挖掘来说，oor与k近邻方法(knn)也有着密切的联系。k近邻方法的核心思想是由某些特性定义的相对“亲密邻居”的对象趋向于有类似的结果。从本质上说，可以用k近邻方法进行预测，而不用做任何建模假设，因此该方法也被称为非参数预测方法。然而，k近邻方法的效率没有其它建模方法高，其原因之一是它并没有考虑到这样一个事实，即许多邻居具有同等疾病关联(即结果关联)(无论是无效假设或备择假设)，而通过邻居的组合是可以提高预测精度的。相比之下，oor利用周边信息(即，相似性度量)与多个包含信息的范例关联。在概念层面，oor可以被看作是k近邻回归函数估算的延伸。
另一种密切相关的方法是成员等级分析法，简称为gom。从概念上讲，gom通过引入一组潜在成员变量，假设该潜在成员变量的分布是合理的，gom对结果的联合分布以及协变量建模，在整合了所有gom潜在成员变量后，可推导边际似然用于估算和参考，而不是单一协变量的特定边际解释。在此方面，oor类似gom，利用相似性度量获得分析目标，但其在建模假设和相关实施上是不同的。oor的主要优点在于，无需假设潜在成员的分布，而完全基于经验证据进行推断。
oor在下述两方面有很大的发展。首先，在逻辑回归模型下构建oor，用于二元疾病表型的应用。通过将逻辑回归扩展到广义线性模型，oor可以应用于与其它类型的表型的研究，如连续、分类或截尾的表型，并适当选择关联函数，对表型和协变量的关系进行建模。第二，在其他类型的复杂协变量，例如文本串(例如，来自网络搜索)、电子信号或二维图像方面的应用。此外，协变量可以是高维数据，其维度的数目可远远大于样本量。对于这些不同的应用，关键是要选择背景相关的相似性度量，来定义对象之间的关于其协变量特性的“相似性度量”。研究oor的长期目标是使其适用于大数据所产生的各种表型与各类协变量。
显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于本领域技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无法对所有的实施方式予以穷举。凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。