最近《山河令》牵动了数亿人的心，wink透过易容术，看透了周美人的骨 相，这一眼万年，让人好生留恋，故曰： 山河不足重，重在遇知己。 今天也带大家透过药物反应数据库和肿瘤数据库表面的样子，深层刨析，看看其骨相究竟如何清秀。看完这篇烧脑的文章,就放假继续追剧啦!
期刊名：Nature Communications2041 -1723 2020 年影响因子/JCR分区：12.121/ Q12010 4316
给定一张图 x，神经网络可以告诉你其中是什么内容，或者说属于什么类别 c。 生成模型分为 VAE 和 GAN 两条路（其实还有一条值得继续的路是 PixelCNN/RNN……），十分有趣的是，它们拥有互补的长处和短处，如下图: 以上都加上了类别作为输入，所以称之为 CVAE 和 CGAN。直观地看： 所以刚好是互补。这其实来自于它们各自的构造。首先，在构造生成模型时，人们渐渐发现，最常用的是下面四个网络： E：编码器 Encoder。给定一张图 x，可以将其编码为隐变量 z（且希望满足高斯分布）。如果还给定了类别 c，那么生成的隐变量就会质量更高（更随机）。 G：生成器 Generator。给定隐变量 z（如随机噪声），就可以生成像模像样的图像。如果还给定了类别 c，那么就会生成像模像样的属于类别 c 的图像。 C：分类器 Classifier。给定一张图 x，输出所属类别 c。这是大家的老朋友。 D：辨别器 Discriminator。给定一张图 x，判断它是真实的图片，还是“电脑乱想出来的”。这是 GAN 首先引入的网络，它可以和 G 左右互搏，互相进步。 回归到这篇文章，这里应用VAE和Generator结合对 药物处理肿瘤数据的分类学习 。抛开图像处理不谈，这个VAE被称为 变分自动编码器 ，VAE是以非监督的方式使用神经网络学习复杂数据分布的最流行的方法之一，例如使用神经网络。这是一个植根于贝叶斯推断的概率图模型，即该模型旨在学习训练数据的基本概率分布，以便它可以容易地从该学习分布中采样新的数据。这个想法是学习一个被称为潜变量的训练数据的低维潜在表示（我们假设已经产生了我们的实际训练数据的潜变量（这些变量不是直接观察到的，而是通过数学模型推断的）。这些潜在变量可以存储有关模型需要生成的输出类型的有用信息。这里数据分布形式采用高斯分布模型作为先验，学习时采取训练数据可能性最大化。在推理时，我们可以简单地从N（0,1）中采样z并将其馈送到解码器网络以生成新的数据点。 VAE最厉害的地方就是既能学习生成模型和又可以推理模型 。 ------------------------------ ------------------- 由于肿瘤间和肿瘤内的异质性，癌症患者的药物反应显著不同。并且就像之前我们推文中提到的，肿瘤微环境对治疗效果起着很重要的作用。作者 采用VAE模型压缩低维空间的数千个基因特征向量，并且这个编码向量可以准确推算药物反应 ，表现优于基于signature-gene分析方法，并适当地控制过拟合问题。基于该模型体现的药物反应，发现了表达相关组分(the expression-regulated component ,EReX)高度相关性。利用这个模型，研究人员对TCGA数据库肿瘤数据进行药物反应评估，并找到包括细胞系起源、体细胞突变和肿瘤突变负荷、肿瘤微环境和混杂因素特征与药物反应的相关性关系。文章信息量很足，一起来学习一下吧！ 基线数据 -- The Cancer Cell Line Encyclopedia (CCLE) 获取1156个细胞系的基因表达数据（RNA-seq），每个细胞系都有其原始的细胞系。检测到的细胞系可以与多种癌症类型相匹配，包括实体癌症、造血和淋巴组织。少于20个样本的细胞系将被排除，然后最终得到19个细胞系的1100株细胞系建立VAE模型。利用表达最显著的基因构建VAE模型，最终选取他们对应的细胞系进行CCLE和GDSC药物反应预测。 药物反应数据 -- CCLE提供了504个细胞系的24种抗癌化合物，包括突变和药物反应（“ActArea“），GDSC提供了同一细胞系评估251种化合物的药物反应。 在Elastic Net 模型 中使用log-transformed IC50（LN_IC50）计算。 TCGA 多组学数据 -- 数据从USCS Xena获取，分别下载了有体细胞突变，CNV和mRNA表达的TCGA样本，CNV数据来自Affymetrix 6.0平台，有5个级别代表不同的CNV状态: deep deletion（CN=-2）， copy loss (CN=−1)，neutral (CN=0)，copy gain (CN=1)和amplification (CN=2)。RNA-seq的mRNA表达数据以log2(RPKM + 1)值的形式下载。 基因表达数据的预处理 -- 所有RNA-seq数据都使用log2(RPKM + 1)格式进行转换。采用基于rank-reversed percentile的变换对RNA-seq和微阵列数据进行预处理。 VAEN模型训练 -- 如图1A，这里实现了一个三层的VAE模型，包括输入层-编码器、中间层-解码器和输出层(图1A)。VAE模型这里使用了Python的深度学习库- Keras (version 2.1.)辅助TensorFlow backend。编码器是用一个平均向量和一个标准偏差向量分别对输入向量进行编码的过程，然后进行非线性变换（ReLU;S activation）。对于VAE模型，每个药物回归模型都依据Elastic Net进行lambda的5倍交叉验证。采用平均R2 在剩余样本中选择模型，使用CCLE药物反应和CCLE基因组学数据对CCLE模型进行训练。使用GDSC药物反应和CCLE基因组学数据训练GDSC模型。 相关性测试 -- 利用Wilcoxon秩和检验评估基于体细胞突变的药物-基因相关性，比较每个基因突变样本与野生型样本的预测药物反应。基于基因表达的药物-基因相关性采用Pearson相关系数评估。 CCLE ：https://portals.broadinstitute.org/ccle.GDSC：https://www.cancerrxgene.org/.TCGA：UCSC Cancer Genome Browser Xena文章源代码: https://github.com/bsml320/VAEN/
图1，整体研究流程框架。A，VAE模型；B，训练药物反应的回归模型流程；C，t-SNE plot展示细胞系极其相关转录组；D，分解观察到的(测量到的)药物反应；E，CCLE数据库细胞系和癌症类型测量和估算药物反应图解
如上图1展示，首先，利用基线基因表达数据建立VAE模型（图1A，C）对于深度神经网络生成模型利用CCLE和GDSC测定的IC50数据建立回归模型来估算药物反应（图1B），并用该模型应用于TCGA和其他临床数据集进行验证（图1E）。这些预测模型是通过使用Elastic Net strategy (VAE model followed by Elastic Net, or VAEN)实现。从概念上讲,观察到(测量)药物响应可以分解为expression-regulated 药物（EReX）或者其他的解释元素（基因变异，甲基化，miRNA ，lncRNA等）和不确定元素（如批次效应，试验条件等）（图1D）。VAEN可以用来估算EReX，理论上，无论测量平台或预测方法如何，每种药物的EReX都是相同的。 来自CCLE和GDSC项目的原始细胞系代表了广泛的细胞系谱系，如上皮细胞、间充质细胞和造血细胞等。根据提到的筛选条件，获取20个样本以上的细胞系数据，和细胞系中表达差异最大的基因用于构建VAE模型（n=6163）。参数扫描用来评估参数的准确性，这里测试了潜在大小、学习速率、批效应和代数。利用原始数据和重建数据之间的损失来评估拟合模型。图1C所示即为t-SNE分布的数据和观察到的主要簇的细胞系，其中造血和淋巴组织来源、成纤维细胞、皮肤、胃和乳腺细胞系可以形成明显的簇。然而，这里没有观察到与任何药物显著相关的聚类。 接下来谈一下这个VAEN模型。VAE模型本身的拟合是为了实现输入数据的降维，但本身并不适合药物反应预测。与传统PCA不同，数据降维不是唯一的表示目的，对于相同的输入矩阵，可能存在很多矩阵代表输入数据，其中每一个矩阵都可能与其他矩阵略有不同，它们分别适用于一种或某些药物，而不是所有的药物。因此这里生成100个VAE模型为接下来的药物反应归因建立一个库。Elastic Net用来对每种药物进行降维，利用每种VAE模型的潜在向量训练VAEN预测模型。他们使用固定的alpha=0.5，并进行已有的lambda 5倍的交叉验证。然后，他们通过标准多元线性回归的十倍交叉验证来评估选定的潜在向量(即EN模型中系数非零的向量)。模型效率以样本中样本的平均决定系数R2 来衡量。在100个VAE模型(即100个潜伏矩阵)中，接下来的分析中，会选择每种药物的平均滞留R2 最高的模型。 图2，评估模型效率。A，所有药物的样本内皮尔森相关系数(PCC)分布；B，评估CCLE和GDSC小组共同的14种药物；C，利用TCGA癌症数据评价14种共同药物；D，使用共同化合物17-AAG(在GDSC中也称为坦奈斯比霉素)的例子演示观察到的和预测的药物反应。
在将原始基因表达数据输入成VAE模型之前，他们测试了几种 预处理方法。(1)样本间每个特征的z-score归一化，并将每个特征缩小到[0,1]范围。(2)每个样本所有基因的z-score归一化。(3)基于秩的逆变换。综合评估后，基于秩的逆变换和s型激活的VAE模型。图2A中显示，皮尔逊相关系数(PCC)测定的24种CCLE化合物的自预测精度在0.38 (LBW242)和0.77(伊立替康)之间。GDSC的251种化合物的自预测精度在0.26 (Avagacestat)和0.82 (AZ628)，其中203/251(80.88%)化合物PCC >0.5。 然后使用CCLE和GDSC检测的14种药物对预测结果进行评估。发现，CCLE IC50与GDSC IC50的相关性不如CCLE ActArea与GDSC LN_IC50的相关性，因此，他们选择CCLE ActArea来表示药物反应。在比较各细胞系的预测药物反应时，发现观察到的和预测到的药物反应呈正相关(图2B)。当一种药物在CCLE和GDSC之间表现出高度一致的药物反应(例如，在CCLE和GDSC中观察到的药物反应在227个细胞系，PCC (nilotinib) = 0.75之间)，预测数据也呈现出高一致性的趋势（PCC(nilotinib)= 0.8预测药物反应）。值得注意的是，对于每种药物，在CCLE预测模型和GDSC预测模型中选择的效应估算的VAE模型不同。对于构象不一致的药物结构，不同的VAE模型在不同的模型中有不同的投影。所以图2B中的结果表明，这个VAEN模型在各个研究小组中具有可重复性。正如上面提到的(图1D)， EReX是可以被转录组解释的组分。因此，图2B的结果表明，这个基于基因表达的低维模型，在估计EReX方面取得了较高的性能。图2C使用TCGA评估CCLE和GDSC中药物反应，原药浓度高的药物，预测的药物反应趋于高度一致，而图2D正是显示这个实例。当使用拟合的CCLE和GDSC模型比较每种癌症类型的预测药物反应时，发现50.65% (234/462) 的癌症-药物配对具有>0.5的相关性。Bonferroni校正后，在≥20种肿瘤类型中所有药物均表现出显著相关性(P < 0.05/462 = 1.08 × 10−4) (图2E)。
图3，33种癌症类型中药物反应预测。A，观察到的和预测的药物反应在CCLE中的分布；B，每种药物的敏感和不敏感样品的定义；C，对30种癌症(不包括三种免疫相关癌症:DLBC、LAML和THYM)敏感或不敏感样本使用CCLE药物进行富集试验结果；D-E，敏感样本和不敏感样品；F，GDSC具有代表性药物对30种癌症(不包括三种免疫相关癌症类型:DLBC、LAML和THYM)的敏感或不敏感样本进行富集试验结果。 接着，将VAEN模型应用于CCLE细胞株(n = 1100)和TCGA样本(33种癌症类型，n = 10459)。预测了每个样品对24种CCLE化合物和251种GDSC的药物化合物反应。这里使用CCLE演示模型性能，如图3A所示，预测的药物反应与原始数据相似。17-AAG、伊立替康、紫杉醇、PD-0325901和拓扑替康等化合物的反应范围相对较大，而其他化合物的反应范围相对较小(如AEW541、厄洛替尼、Nutlin-3、PLX4720)分布较窄。用CCLE模型将预测药物反应最强的前5% TCGA样本定义为敏感组，将预测药物反应最弱的5% 样本定义为不敏感组(图3B)。使用ActArea作为药物反应的测量。对敏感组和不敏感组的每种癌症类型的样本进行了富集分析(Fisher精确检验)。如图3C所示每种癌症类型在敏感组和不敏感组均观察到富集模式。如图3D所示(利用CCLE模型预测ActArea治疗TCGA样品)和3E(利用GDSC模型预测-LN(IC50)治疗TCGA样品)，有较高比例的样品在SKCM (CCLE: 96.20%;GDSC: 97.01%)和在UVM(CCLE: 98.75%;GDSC: 98.75%)对PLX4720敏感。
图4，验证TCGA数据。A，采用HER2免疫组化方法，比较不同亚型TCGA-BRCA患者对拉帕替尼的预测反应；B，5-氟尿嘧啶处理的TCGA-STAD样本的生存分析；C，紫杉醇处理TCGA-BRCA样品的生存分析；D，按反应状态分层的TCGA-BRCA样品对紫杉醇的预测反应比较；E，对两种MET抑制剂(PF2341066和PHA665752)的预测反应与MET扩增(顶部面板)或MET表达的关联(下面板)使用CCLE模型；F，利用GDSC模型预测三种MET抑制剂(crizotinib(也称为PF2341066)、PHA665752和foretinib(仅在GDSC中测试))的反应与MET扩增(上panel)或MET表达(下panel)的关联; G，利用先前报道的基因标记对AZD6244、厄洛替尼和AZD0530进行验证。 接着使用注释良好的靶点药物进行进一步验证，证实了药物反应预测具有很高的准确性。 这里使用药物有拉帕替尼，MET抑制剂(克唑替尼/PF2341066，福替尼，和PHA665752)， BRAF抑制剂(PLX4720和RAF265)和几种MEK抑制剂。 如图4A所示，使用CCLE模型(P=2.60 × 10 −19 )和GDSC模型(P=2.52 × 10 −9 )预测lapatinib (ERBB2抑制剂)药物反应与TCGA-BRCA中的HER2免疫组化状态有显著相关性。 在CCLE,测定了两个MET抑制剂(PF2341066和PHA665752)，两者均与MET表达增加显著相关(PF2341066的P=9.84 × 10 −11 ，而PHA665752的P=1.75 × 10 −4 ，图4E)。 在GDSC模型 中,三个测定MET抑制剂(PF2341066，称为克唑替尼GDSC, foretinib和PHA665752)。 这些化合物的药物预测响应也与MET表达或MET扩增显著相关(图4F)。 总的来说，EGFR/KRAS/HRAS/NRAS/BRAF突变的样本在大多数癌症类型中倾向于增加对相应的激酶抑制剂的敏感性，尽管一些相关性没有达到显著性。 这可能是由于复杂的肿瘤微环境或者如突变、CNVs和其他组学水平并发症等混杂因素造成的。 图4B所示，接受5-氟尿嘧啶治疗的STAD样本，预测对该药物反应高的组比预测反应低的组有显著更好的生存结局(P = 0.0023)。 综合图4C来看，尽管紫杉醇处理的完全或部分药效反应的组生存显著高于病情稳定或临床进展组(P = 0.03，图4C)但是BRCA样本还是具有一定的边际效应(P = 0.20，图4C)。 使用gene expression signatures验证 为了进一步验证输入的药物反应预测情况，这里还收集了以前报道的药物基因特征。使用线性回归拟合模型估算了已知的基因特征和药物反应关系。对于每种癌症中的每个基因将样本分成三类，Q25,Q25-Q75和Q75用来表示表达降低四分之一，表达降低一半和表达更高的范围。拟合线性回归模型，将组变量作为定量类型，使模型评估药物反应在三组样本中变化的趋势。以AZD6244为例，之前的一项研究报告了18个基因特征对AZD6244敏感。图4G的GSEA显示，这些基因与对AZD6244的反应显著相关(P = 1.44 × 10−4 ，归一化富集评分(NES) = 2.30)。对于AZD0530(也被称为saracatinib，一种Src和Abl抑制剂)，在CCLE中与药物显著相关(P= 3.89 × 10−4 ,NES= 1.91)，在GDSC模型中与药物显著相关 (P= 0.059, NES= 1.48)。
图5，独立数据集验证。A，使用厄洛替尼治疗的GSE33072样本的生存分析；B，使用数据集GSE32989比较上皮样细胞系(n = 44)或间质样细胞系(n = 25)对厄洛替尼的预测反应；C-E，vemurafenib治疗的黑色素瘤样本生存分析；D-F，比较PLX4720在亲本细胞系和衍生抗性细胞系中的预测反应；G，采用pCR方法比较BRCA亚组对紫杉醇的预测反应。 为了验证VAEN模型，研究人员还收集了6个带有治疗和生存注释的数据集。如图5A所示，用估算药物反应的中位数对样本进行分层，发现高反应组与低反应组之间存在显著差异。数据集是由微阵列平台生成的。然而，由于模型使用基于秩的归一化，可以方便地转换数据，并使用训练过的VAEN模型进行预测。研究表明这些细胞系具有74个基因的EMT特征。重复了原来的研究，将这些细胞系分为上皮样亚组(n = 44)和间叶细胞样亚组(n = 25)。如图5B所示，所有四种VAEN模型(CCLE A-model: P = 6.88 × 10−11 ,CCLE Smodel: P = 3.61 × 10−8 ,GDSC A-model: P = 6.59 × 10−5 ，和GDSC S-model: P = 4.39 × 10−4 )显示两亚组间差异显著，上皮样组反应较好，与文献报道一致。如图5所示，对PLX4720反应较高的样本生存状态较好(CCLE S-model P = 0.031, GDSC A-model P = 0.0049，图5C, E)。亲本细胞系的反应显著高于配对抗性细胞系(CCLE s模型P = 8.032 × 10−4,GDSC模型P = 0.017，paired t-test，图5D, F)。图5G应用紫杉醇药物对于PCR组的响应，显著高于RD组。综上所述，这些结果证明我们预测的药物反应在临床数据上是可靠的。
图6，药物相似性。A, 四组样本的24种CCLE药物的层次聚类结果；B，热图显示了各种药物与CCLE中肿瘤突变负荷(TMB)之间的关联模式；C-D，在使用两种细胞毒性药物后，有反应者(25%反应最高的样本)和无反应者(剩下的75%样本)中显示TMB分布；E，GDSC药物与TMB的相关性分布；F-G，药物种类在药物-癌症类型相关组中的富集与负相关模式(E中的蓝点)和正相关模式(E中的红点)。 对于24种CCLE化合物，使用三组药物反应进行了层次聚类，预测在CCLE和在TCGA样本药物反应(ActArea)。 如图6A中，两种EGFR抑制剂，两种MEK抑制剂，两种Raf激酶B抑制剂和三种细胞毒性化合物具有一致的聚类模式。 研究发现具有细胞毒性的药物与较高的肿瘤突变负荷(TMB)相关。 在CCLE中，细胞毒性化合物(伊立替康和拓扑替康)具有相似的药物反应谱，通常聚集在一起。 通过比较每一种癌症类型中有反应和无反应的TMB (Wilcoxon rank-sum test, two-sided)，发现在多种癌症类型中，对细胞毒性化合物(伊立替康和拓扑替康)有反应的患者有增加TMB的倾向(图6B)。 这一趋势在其中癌症中也观察到（如图6C,D）。 而GDSC中的化合物有 22个化合物被分为一组，另一个为单独一组，但均发现与TMB相关（图6E）。 与TMB增加相关的化合物，富集发现与DNA复制(P = 8.88 × 10 −6 )和细胞周期(P = 0.019)相关，相反，与TMB降低相关的化合物，富集发现EGFR信号(P = 8.84 × 10 −3 )和PI3K/mTOR信号(P = 0.037)相关（图6G,F）。
图7，体细胞突变与药物反应的关系。A，药物-基因关联火山图；B，热图显示每种癌症类型中显著的药物-基因关联；C，热图显示由突变簇决定的显著药物-基因关联；D，对药物AZD6244的反应在两个Ras基因突变簇的样本中的分布；E，药物反应与拷贝数增加的关系；F，药物反应与17号染色体拷贝数增加的关系。 为了识别与输入药物反应相关的基因组特征，这里使用TCGA中的多个组学数据进行关联分析。对于体细胞突变，研究了每种癌症类型中有害的SNV与每种药物的潜在关联。此外，考虑到一些靶向药物对激酶结构域的特定突变具有较强的特异性，研究人员根据突变簇在蛋白序列中的位置进行了关联测试。图7A所示，观察到24个CCLE化合物(benjami - hochberg (BH))中有2343个敏感基因-药物关联和288个不敏感关联（adjusted P或PBH <0.05)。在体细胞突变与药物反应增加显著相关的基因中，我们发现了几个来自PAM通路基因FGFR3,HRAS, KRAS, BRAF, NRAS, TP53, KIT, NFE2L2 （图7B）。 接着探究了突变簇和药物反应之间的关系。发现了723种簇-癌-药物的关联(PBH < 0.2，定义为> 200k碱基，排除长基因的关联)，其中24种癌症类型中有32种unique基因(图7C)。HRAS、KRAS和NRAS的Ras蛋白中G12/G13位点的突变对MEK抑制剂AZD6244和PD-0325901的敏感性增加(BLCA, CESC, HNSC, PAAD, THYM, UCEC, UCS)(图7 D)。这可能是由于癌症样本的异质性，如亚型，突变及混杂因素导致的。 图7E显示基因-药物关联是CNV发生在特定的染色体上，特别是染色体3,7,8,10,12,17染色体。7号染色体上EGFR所在区域的扩增主要与LGG的激酶抑制剂相关。而EGFR抑制剂拉帕替尼与ERBB2扩增(chr17, P=2.30 × 10−27 图7F)相关。
图8，肿瘤异质性的药物反应。A，所有药物的基因t得分分布，突出显示间质基因；B，基因t-score在TSI特征；C，SKCM中BRAF突变与PLX4720应答关系的研究。 肿瘤组织被细胞包围并浸润，统称为肿瘤微环境(TME)。TME由细胞外基质(ECM)、成纤维细胞、神经内分泌细胞、脂肪细胞、免疫细胞和炎症细胞组成，几乎参与了肿瘤发生的每一个步骤。这里发现一些TME与要药物反应相关，首先，癌相关成纤维细胞(CAF)基因和细胞外基质(ECM)基因与拉帕替尼的敏感性呈负相关，即CAF/ECM基因表达增加与拉帕替尼的不敏感性相关。即HER2抑制剂的反应取决于肿瘤样本的微环境。如图8所示，这里研究了18个之前报道过的T细胞炎症预测基因，几种药物与TIS基因呈强阳性相关性(伊立替康、尼洛替尼、PHA665752, PLX4720, and RAF265) 而其他药物(如:AZD6244和PD-0325901)与TIS基因无相关性。由于TIS基因的表达高度提示免疫热或免疫冷状态，这些结果表明具有高水平t细胞炎症的样品(如免疫热)可能对PLX4720和RAF265更敏感。此外图8C中，在SKCM中，当活动模型中包含“下调ERBB2 ERBB3信号通路”，交互作用项具有统计学意义(P = 1.79 × 10−5 )，且仅在信号通路表现低活性的Q25组，即突变组BRAF对BRAF抑制剂敏感(P = 4.92 × 10−3 )。而在THCA中，对PLX4720的应答是由癌症亚型驱动的，Braf样亚型的应答显著高于Ras样亚型(A-model: p = 2.72 × 10−16 ;S-model:P = 8.31 × 10−4 ;paired t-test，图8D)。这些结果表明，癌症样本中的药物反应比细胞系模型中复杂得多，许多混杂因素可能会影响药物反应。 ---------------------------------------------- 至此文章的研究结果和分析就已经结束。由此可见热门的机器学习分类预测方法对于挖掘现有公共数据资源的数据潜藏规则，具有重要意义。文章重新探索报道过的基因特征和药物反应关系，并结合细胞系预测模型，进而拓展到复杂肿瘤环境，验证药物反应的预测情况。 目前发现突变，TMB等仍然是该模型预测的一个壁垒，所以下一个方向，是否可以模型优化，如调整参数模型，加入协变量等方式，非线性拟合反应条件，获取更好的结果呢？ 如此计算模型，便像是手握白衣剑，驾驭流云九宫步如蜻蜓点水，坐看云舒……你不妨来试试。 有生信相关问题联系：18501230653
