*癌症特征的K-均值和聚类模型

在这里，正如【Kakushadze和Yu，2016b】中所述，我们将遵循一种务实的方法，做一些简单的事情——数据中存在着大量的噪音，认为做复杂的事情根本不会带来回报。因此，作为第一个切入点，我们可以采用XIS=ln（1+Gis）（3）这考虑到Gis=0的情况；对于Gis 1我们有Ris≈ ln（Gis），根据需要。其次，【Kakushadze和Yu，2016b】的详细实证分析揭示了其中所称的“总体”模式在发生计数数据中的明确存在。这种“整体”模式被解释为体细胞突变噪声未相关这是因为事情几乎是随机的，手头上唯一的“分布”是flat。在融资方面，类似于这一模式的是所谓的“市场”模式（参见【Bouchaud和Potters，2011】及其参考文献），对应于广泛市场的整体运动，这是为了（事实上模糊）真正的潜在癌症特征，因此必须以某种方式加以考虑。下面是理解“总体”模式的简单方法。设相关矩阵ψij=Cor（Xis，Xjs），其中Cor（·，·）是序列相关。一、 e.，ψij=Cij/σIσj，其中σI=Ciiare方差，序列协方差矩阵xcij=Cov（Xis，Xjs）=d- 1dXs=1ZisZjs（4），其中Zis=Xis-夏尔连续降级，而“Xi=dPds=1Xis”的意思是“Xi”。平均成对相关ρ=N（N-1） PNi，j=1；不同突变类别之间的i6=jψij是非零的，事实上，对于我们研究的大多数癌症类型来说，i6=jψij是高的。这就是前面提到的体细胞突变噪音，必须加以考虑。如果我们按癌症类型（见下文）对样本进行聚合，并计算相关矩阵ψij以获得聚合数据（在我们研究的14种癌症类型中，见下文），则平均相关性ρ超过惊人的96%。

另一种思考方式是，不同样本（或癌症类型，如果我们按癌症类型聚合样本）中的发生计数并没有在所有样本（癌症类型）中统一标准化。因此，在vanilla matrix Gis上运行NMF、聚类或任何其他特征提取算法（或（3）中定义的“log”xis）将相当于苹果和桔子的混合，从而掩盖真正的潜在癌症特征。继【Kakushadze和Yu，2016b】之后，分解出“整体”模式（或“去噪”矩阵Gis），因此最简单的是对矩阵Xis的横截面（即，跨越96个突变类别）贬低。也就是说，我们使用Xis代替Xis，这是通过贬低Xis的列获得的：Xis=Xis- Xs=Xis-NNXj=1Xjs（5）我们应该注意，使用Xisin而不是Xisin（1）不会影响聚类。实际上，g in（1）在形式Xis的变换下是不变的→ Xis+s、 在哪里sis是任意的d向量，因此我们也有Yas→ Yas+s、 soXis公司- Yasis保持不变。事实上，这很好：这意味着去噪不会给聚类本身带来任何额外的错误。然而，矩阵中的实际权重会受到去噪的影响。我们讨论了fixingwiabelow的算法。然而，在确定权重之前，我们还需要一种成分，通过这种额外的成分去噪可以影响聚类。影响所有股票（不同程度）-流入（流出）市场的现金往往会推高（降低）股价。这是市场风险因素，为了缓解这种风险，可以持有中性股票组合（即持有相同的美元多头和空头头寸）。在本文中，“serial”指的是“over the index s”。Cij的总体标准化，即d- 1（无偏估计）vs。

对于我们的目的而言，定义慈进（4）的分母中的d（最大可能性估计）是无关紧要的。因此，在这种情况下，d=n=14 in（4）。基于上述原因，我们应该使用Xis，而不是Gis。2.4.1【Kakushadze和Yu，2016c】中讨论了标准化对数计数，聚类Xis（或等效Xis）将是次优的。问题是这样的。设σibe系列标准偏差，即（σi）=Cov（Xis，Xis），其中，如上所述，Cov（·，·）是系列协方差。这里我们假设样本是按癌症类型聚合的，因此s=1，d，d=n=14。现在，σi在横截面上并不均匀，并且在突变类别中有很大的变化。σiis的密度如图1所示，并且是倾斜的（有尾的）。σireads总结：最小值=0.2196，第一个Qu.=0.3409，中位数=0.4596，平均值=0.4984，第三个Qu.=0.6060，最大值=1.0010，SD=0.1917，MAD=0.1859，偏斜度=0.8498。如果我们简单地将Xis聚类，那么σi的这种可变性将无法解释。一个简单的解决方案是将标准化的demeaned log countseXis=Xis/σi而不是Xis进行聚类。通过这种方式，我们将不均匀（和倾斜）的标准偏差从对数计数中剔除。请注意，现在去噪确实会产生差异，包括。

事实上，如果我们使用exis=Xis/σi（回想一下σi=Cov（Xis，Xis））而不是ofexis=Xis/σiin（1）和（2），那么数量g（以及聚类）将不同。2.5固定聚类数既然我们知道要聚类什么（即eXis）以及如何获得“唯一”聚类，我们需要弄清楚如何确定（目标）聚类数K，这是我们上述算法的输入之一。在【Kakushadze和Yu，2016b】中，有人认为，在癌症特征的背景下，可以通过建立统计因子模型来确定它们的数量【Kakushadze和Yu，2017b】，即特征的数量只是统计因子的数量。因此，出于同样的原因，我们在这里确定了聚类算法中的（目标）聚类数，其中统计因子的数量通过[Kakushadze和Yu，2017b]的方法确定。2.5.1有效RankSo，继【Kakushadze和Yu，2017b】和【Kakushadze和Yu，2016b】之后，我们设定K=圆形（eRank（ψ））（6）更准确地说，【Kakushadze和Yu，2016c】的讨论是在金融背景下进行的，towit，量化交易，有其自身的细微差别（见下文）。然而，其中一些讨论是相当普遍的，可以适用于各种各样的应用。Qu.=四分位，SD=标准偏差，MAD=平均绝对偏差。在其他情况下，已经讨论了各种聚类数量的确定方法，例如，【Rousseeuw，1987】【Pelleg和Moore，2000】【Steinbach等人，2000】【Goutte等人，2001】【Sugarand James，2003】【Hamerly和Elkan，2004】【Lleit'i等人，2004】【De Amorim和Hennig，2015】。在金融领域，这些被称为统计风险模型【Kakushadze和Yu，2017b】。有关多因素风险模型的讨论和文献，请参见【Grinold和Kahn，2000】【Kakushadze和Yu，2016a】以及其中的参考文献。

有关统计风险因素数量的前期工作，请参见【Connor和Korajczyk，1993年】和【Bai和Ng，2002年】。这里，圆（·）可以被floor（·）=b·c代替。这里，eRank（Z）是对称半正定义矩阵Z的有效秩【Roy和Vetterli，2007年】。它被定义为秩（Z）=exp（H）（7）H=-LXa=1paln（pa）（8）pa=λ（a）PLb=1λ（b）（9），其中λ（a）是Z的L个正特征值，H具有（香农a.k.a.谱）熵的含义【Campbell，1960年】【Yang等人，2005年】。让我们强调一下，在（6）中，矩阵ψij是根据已减记的对数x is计算的。eRank（ψij）的含义是，它是矩阵ψij的有效维数的度量，它不一定与其正值的数目L相同，但通常更低。这是因为许多d向量X可以连续高度相关（这通过特征值中的大间隙表现出来），从而进一步降低了相关矩阵的有效维数。2.6如何计算权重？剩下要完成的一件事是弄清楚如何计算权重SWIA。令人高兴的是，在聚类的背景下，与NMF相比，我们有着显著的简化，一旦我们确定了聚类，即矩阵，计算权重就会变得非常简单OhmiA=δG（i），A（或等效地，映射G：{i}7→{A} ，i=1，N、 A=1，K、 其中，为了便于标记，我们使用Kto表示“最终”聚类中的聚类数（见上文）。

正如在NMF中一样，我们希望通过权重矩阵wia和曝光矩阵HAs的乘积来近似矩阵Gisvia，两者都必须是非负的。更准确地说，由于我们必须移除“整体”模式，即去除矩阵Gis的噪声，继【Kakushadze和Yu，2016b】，我们将近似于指数化的去噪对数矩阵Xis：Gis=exp（Xis）（10），我们可以通过取Gis=exp（平均值（Xis）+Xis，orGis=exp（中值（Xis）+Xis来包括整体归一化，或Gis=exp（中位数（Xs）+Xis）（回想一下，Xs是Xis的列平均值的向量-参见公式（5）），等等，使其看起来更像原始矩阵Gis；然而，这并不影响提取的签名。此外，从技术上讲，在重新指数化之后，我们应该“减去”定义（3）中添加的额外1（假设我们包括上述总体正常化之一）。然而，数据中固有的噪声使得这一点毫无意义。请注意，使用标准化的demeaned log countsexis得到相同的ψij。这是因为W的每一列都是权重，将其归一化为1。因此，我们希望将Gisvia近似为一个乘积W H。然而，通过聚类，我们得到了WiA=wiδG（i），a，即，我们有一个块（簇）结构，其中对于给定值a，除了i之外，所有WiAare为零∈ J（A）={J | G（J）=A}，即对于i标记的属于A标记的簇的置换类别。因此，我们的Gis到产品的矩阵分解现在简化为一组依赖性分解，如下所示：Gis≈ wiHAs，我∈ J（A），A=1，K（11）所以，不再需要运行NMF了！事实上，如果我们能够以某种方式确定HAsfor agiven集群，那么在该集群中，我们可以确定相应的权重wi（i∈ J（A））通过一系列线性回归：Gis=εis+wiHAs，i∈ J（A），A=1，K（12），其中ε是回归残差。一、 e.对于每个A∈ {1。

，K}，我们回归了d×nAmatrix[（G）T]si（i∈ d向量HAs（s=1，…，d）上的J（A），nA=| J（A）|），回归系数仅为nA向量wi（i∈ J（A）），而残差是d×nAmatrix[（ε）T]si。请注意，此回归是在没有截距的情况下运行的。现在，这一切都是有意义的∈ J（A））回归使二次误差termPds=1εis最小化。此外，如果Has为非负，则权重wi自动为非负，因为它们由以下公式给出：wi=Pds=1GisHG（i），sPds=1HG（i），s（13）。现在，我们希望将这些权重归一化：Xi∈J（A）wi=1（14）这始终可以通过重新缩放HAs来实现。或者，我们可以选择HAs而不用担心归一化，计算wivia（13），重新缩放它们以使其满足（14），同时相应地重新缩放HAs。任务完成！2.6.1固定曝光膨胀，几乎。。。我们仍然需要弄清楚如何确定风险敞口。最简单的方法是注意我们可以使用矩阵OhmiA=δG（i），Ato将GIS中的指数i替换为指数A，即我们可以取A=ηANXi=1OhmiAGis=eηAnAXi∈上标T表示矩阵变换。也就是说，在归一化常数eηA（通过（14）固定）之前，我们简单地在每个簇中采用GIS的横截面平均值。（回想一下，nA=J（A）。）与所有GI均为阳性一样，sode定义的HAs自动为阳性。因此，宽视野（13）也都是积极的。这是一个好消息–消失的WI将相当于一个完整的权重矩阵WiA（即，一些突变将不属于任何簇）那么，（15）为什么有意义呢？看看（12），我们可以观察到，如果用A标记的每个簇内的剩余ε是横截面的，则是随机的，那么我们期望pi∈J（A）ε为≈ 如果我们这里有一个精确的等式，那么我们会有（15），其中ηA=1（即ηA=nA），假设归一化（14）。

实际上，残余ε并不完全是“随机的”。首先，每个簇中的nAof突变类别数量不多。第二，如上所述，不同突变类型的序列标准差存在差异。这导致我们考虑变化。2.6.2上述变量我们认为，根据系列标准偏差σi中的横截面变异性（和偏度），将标准化去甲基对数计数x=Xis/σ进行聚类是有意义的。当计算Has和wias时，我们可能会担心GIS中的类似影响。这可以通过使用归一化数量egis=Gis/ωi来缓解，其中ωi=Cov（Gis，Gis）是序列方差。也就是说，我们可以定义=eηAνAXi∈J（A）eGis=eηAνAXi∈J（A）ωiGis（16）wi=ωiPds=1eGisHG（i），sPds=1HG（i），s=Pds=1GisHG（i），sPds=1HG（i），s（17），其中νA=Pi∈J（A）1/ωi。因此，1/ωi是集群平均值中的权重。2.6.3另一个变量人们可能会想，考虑到Gis和Gis的倾斜大致对数正态分布，将暴露量与退化对数计数的群内横截面平均值相关联，而不是与Gis的横截面平均值相关联，是否有意义？这很容易实现。因此，我们可以确定（这确保了HAs的积极性）：ln（HAs）=ln（eηA）+纳西族∈J（A）Xis（18）指数化we getHAs=eηA易∈J（A）Gis1/nA（19），即，这里我们假设ε在（12）中是/ωI近似随机的。一、 e.这里我们用几何平均值代替（15）中的算术平均值。如上所述，这里我们也可以引入非平凡权重。请注意，（17）的形式与（13）的形式相同，它只受重量的影响。因此，我们可以在几何平均中引入权重如下：ln（HAs）=ln（eηA）+uAXi∈J（A）eXis=ln（eηA）+uAXi∈J（A）σiXis（20），其中uA=Pi∈J（A）1/σi.回想一下（σi）=Cov（Xis，Xis）。因此，我们有：HAs=eηAYi∈J（A）（Gis）1/uAσi（21）因此，权重是指数1/uAσi。

也可能有其他变化。2.7实现我们现在准备讨论上述算法的实际实现，其中大部分R代码已在[Kakushadze和Yu，2016b]和[Kakushadze和Yu，2016c]中提供。R源代码见附录A。3实证结果3.1数据汇总在下面的实证分析中，我们使用了与[Kakushadze和Yu，2016b]中相同的基因组数据（仅来自已公布的样本）。表1总结了这些数据（借用自【Kakushadze和Yu，2016b】），其中给出了总计数、样本数和数据来源，如下所示：A1=【Alexandrov等人，2013b】，A2=【Love等人，2012年】，B1=【Tirode等人，2014年】，C1=【Zhang等人，2013年】，D1=【Nik Zainal等人，2012年】，E1=【Puente等人，2011年】，E2=【Puente等人，2015年】，F1=【Cheng等人，2016年】，G1=【Wang等人，2014年】、H1=【Sung等人，2012年】、H2=【Fujimoto等人，2016年】、I1=【Imielinksi等人，2012年】、J1=【Jones等人，2012年】、K1=【Patch等人，2015年】、L1=【Waddell等人，2015年】、M1=【Gundem等人，2015年】、N1=【Scelo等人，2014年】。【Kakushadze和Yu，2016b】的附录A中给出了具有相应出版物来源的样本ID。在我们下面的分析中，我们按14种癌症类型对样本进行了汇总。结果数据见表2和表3.3.1.1数据结构。基础数据由一个矩阵组成，称为Gis，其元素为突变类型的发生计数，标记为i=1，在s=1，…，标记的样品中，N=96，d、 更准确地说，我们可以使用一个矩阵GIS，该矩阵GIS结合了不同癌症类型的数据；或者，我们也可以选择使用单个矩阵[G（α）]，其中：α=1，n标记n种不同的癌症类型；如前所述，i=1，N=96；s=1，d（α）。这里d（α）是用α标记的癌症类型的样本数量。

组合矩阵GIS是通过简单地将矩阵[G（α）]按列相加（即自举）得到的。对于我们在这里使用的数据（见上文），这个“大矩阵”有1389列。通常，单个矩阵[G（α）]是，因此，“大矩阵”包含大量噪声。对于某些癌症类型，我们可以获得相对较少的样本。我们也可以有“稀疏填充”的数据，即对于某些突变类别有许多零。如上所述，不同的样品不一定是统一标准化的。等等。底线是数据有噪音。此外，直观性很明显，我们使用的矩阵越大，从统计学上讲，我们应该期望通过任何合理的算法获得更多的“签名”（或簇）。然而，如上所述，大量的签名基本上是无用的，并且违背了从一开始提取它们的全部目的——我们有96个变异类别，因此很明显，签名的数量不能超过96！如果我们最终得到，比如说，50多个特征码，这对我们了解潜在的癌症有什么新的或有用的意义？答案可能只是大多数癌症彼此没有太多共同之处，从治疗应用的角度来看，这将是一个令人失望的结果。为了至少在一定程度上缓解上述问题，继【Kakushadze和Yu，2016b】之后，我们可以按癌症类型聚合样本。这样我们得到了一个N×N矩阵，我们也将其称为Gis，其中索引s=1，d现在取对应于癌症类型的d=n值。