Nat. Commun | 用于全基因组药物重定位的系统网络算法

 

1. 背景

DNA/RNA测序的最新进展实现了通过“精确”定位个性化疾病模块来快速识别新靶标并重新利用已批准的药物治疗异质性疾病。基因组学时代，药物开发已成为高度集成的系统性问题，互补多组学与计算方法成为新的研究范式，由于基因组学和系统生物学最新技术和计算方式的进步，使得利用导致人类疾病的癌症类型特异性机制来识别新靶向药物与治疗药物成为可能。基于网络的方法通过度量药物靶标与人类蛋白相互作用组中疾病蛋白的接近度，为药物重新定位靶标和联合疗法提供了可能性。

为了利用基因组信息实现个性化药物重定位，凯斯西储大学Cheng课题组开发了一种全基因组定位系统网络（GPSnet）算法，通过针对特定患者的DNA和RNA测序图谱的特异性靶标,获取疾病模块来重新利用药物。Cheng课题组研究了来自癌症基因组图谱的15种癌症类型，约5000名患者的全外显子组测序和转录组谱。研究表明，GPSnet预测出疾病模块可以预测药物反应，并对140种已批准药物的新适应症进行优先排序。Cheng课题组通过实验验证了已批准的心律不齐和心力衰竭药物ouabain特异性靶向HIF1α/LEO1-mediated细胞代谢途径，在肺腺癌中显示出潜在的抗肿瘤活性。

2. 方法

 

2.1全外显子组测序数据

Cheng课题组从三个来源下载了患者的正常肿瘤成对体细胞突变数据(1).Elledge Lab website at Harvard University  (2).Sanger website  (3).COSMIC: Catalog of Somatic Mutations in Cancer .为了减少冗余并确保本研究中体细胞突变数据的质量，本研究仅关注上述三个数据集中TCGA肿瘤正常匹配样本体细胞突变,并且从GDC网站下载了15种癌症类型的RNA-seq数据。

 

2.2构建药物-目标网络

Cheng课题组从六个常用数据源中为FDA批准的药物集成了高质量的药物靶标相互作用，并使用报告的亲和力数据定义了药物靶标相互作用,抑制常数/效价、解离常数、中值有效浓度、中值抑制浓度。从DrugBank数据库,治疗目标数据库和PharmGKB数据库获得了药物-靶标相互作用。

 

2.3建立人类蛋白质相互作用组

为了构建当前可用的全面的人类蛋白质相互作用组，Cheng课题组建立了15个常用数据库，其中包含多个实验证据与内部系统的人类蛋白质相互作用组。

 

2.4GPSnet模型简介

GPSnet算法包含两个主要部分，癌症类型特异性疾病模块识别与计算药物再利用。疾病模块应该与从RNA-seq序列中衍生出的共表达癌症类型蛋白质相互作用网络高度相关，其中疾病模块中的基因应倾向为高度突变的基因，通过随机搜索方法来识别GPSnet中的疾病模块。

首先对每一种癌症类型中的每一个基因i初始化一个分数，其中，是在相应的癌症类型中基因的突变频率，是基因中cDNA的长度，为了消除稀疏的体细胞突变数据的影响，使用网络平滑方法在整个网络上传输分数。其次，重新启动随机行走过程（RWR）用于计算平滑基因得分，随机游走开始于基因i,随机游走将有的概率随机移动到基因i的邻居基因，或者在每个迭代时间步中有概率回到原来的出发的基因i，其中是随机游走过程中的驱动重启概率。研究中使用来平衡度偏差，t时间步的向量分数与传播过程定义如下：

（1）

其中是每个基因初始化分数向量，是转移矩阵，如果基因i与基因j有相互作用则，如果没有则，其中是网络中基因j的邻居节点数，传播过程的理论解是:

（2）

其中向量中第i个元素是基因i中的平滑分数。

识别出的疾病模块定义为特定癌症类型共表达PPI网络中的一个子图，模块M的分数根据式(2)计算，其中m代表模块M中基因的数据量，是对应癌症类型的整个基因集的平均得分。以下步骤用于执行随机搜索过程以生成原始模块。最初，随机选择一个基因作为“种子”基因。

在人类蛋白质相互组中，使用表示模块M中相互作用的基因集。对于每一个属于集合中的基因i，使用式(4)计算其在特定癌症类型PPI网络中的连接度

（3）

（4）

对于每一属于集合的基因i，如果其基因i被添加到模块中，通过式(5)计算器扩展模块分数，重复搜索步骤，直到没有更多的基因可以添加到相应的模块中，从而获得一个原始模块。这项研究中，Cheng课题组为每种癌症类型构建了大约60000个原始模块，使用这样的方式，人类蛋白质交互作用网络的每一个基因被随机选择了五次。删除了少于10个基因的原始模块，并收集了得分最高的前1％的模块（约300个）。通过计算每个基因出现在这些模块中的次数来表明基因的可信度。最后，选择了基因置信度值大于0.5％的基因，在特定癌症类型特异性共表达PPI网络中，将这些基因中最大的连接元件组装起来，作为每种癌症类型的最终疾病模块。

（5）

 

2.5药物基因组学模型

Cheng课题组从GDSC数据库下载了癌细胞株中基因批量表达谱和药物反应数据.使用默认参数与线性核的LIBSVM R软件包构建回归模型，预测药物的IC 50值。

 

2.6基因集富集分析（GSEA）

Cheng课题组从Connectivity Map中收集药物基因标签。CMap包含来自于培养的人类细胞株的7000多种基因表达谱，在不同浓度下，这些人类细胞株用不同的生物活性小分子进行了处理，总共涵盖6100个个体。因此，CMap提供了针对给定探测集差异表达程度的度量。幅度（a）定义公式(6)所示,其中t是药物治疗组的标度和阈值平均差值，c是对照组的阈值平均差值,因此，a = 0表示无差异表达，a> 0表示治疗后表达增加，而a <0表示治疗后表达减少。

（6）

 

2.7网络邻近度

给定疾病蛋白集S与药物靶标集T，，在人类蛋白质相互组中，结点s与最接近疾病蛋白t之间的平均最短路径长度定义为最近的距离 测得的最接近距离d（s，t），表达式（7）式:

（7）

为了评估药物与特定疾病之间网络距离的重要性，研究人员构建了一个参考距离分布，在网络中，该距离分布与原始疾病蛋白和药物靶标大小和程度分布相同的两组随机选择的蛋白之间的预期距离相对应。

 

2.8途径富集分析

Cheng课题组使用ClueGO对经典KEGG途径中基因进行了富集分析。

 

2.9.细胞培养

所有细胞均在37°C热氧化碳浓度在5%的潮湿培养箱中培养。NSCLC细胞株A549，H522，H596，H1975，HCC827和PC9从American Type Culture Collection 获得，并在Roswell Park Memorial Institute (RPMI) 1640 Medium培养基中培养，培养基中添加了10％胎牛血清（FBS， Gibco）和青霉素链霉素。肺正常细胞株MRC5和WI38获自the Shanghai Cell Bank of the Chinese  Academy of Sciences ，并维持在Eagle最低基本培养基中，该培养基补充了10％胎牛血清（FBS，Gibco）和青霉素- 链霉素。对细胞系株进行支原体检测测试，并通过短串联重复序列（STR）进行鉴定。

 

2.10细胞活力测定

以3000–5000个细胞每孔的接种方式，将细胞接种到96孔板中12 小时，然后与所示化合物培育48 小时。根据相关协议，使用Promega的CellTiter 96 AQueous One Solution检测细胞活力。使用Graphpad Prism 7从剂量反应曲线计算的值。

 

2.11菌落形成

将A549或H522细胞以每孔1500个细胞的密度接种到6孔板中，加入2 ml的1640培养基中，其中添加了10％FBS和青霉素链霉素。细胞粘附后，将不同浓度的ouabain与细胞结合。每2天更换一次培养基，7天后，将菌落固定在4％多聚甲醛中，并用0.2％结晶紫染色。

 

2.12免疫印迹和抗体

用含有100 mM Tris-HCl、150 mM NaCl、1 mM EDTA、0.1％SDS、1％脱氧胆酸钠、1％Triton X-100蛋白酶抑制剂与磷酸酶抑制剂的混合溶液溶解细胞。使用BCA测定蛋白质浓度。将等量的蛋白质样品在SDS聚丙烯酰胺凝胶电泳上电泳，并转移到硝酸纤维素膜上。将膜在含5％BSA的TBST中密封，并用一抗和相应的荧光偶联二抗探测。

 

2.13实时定量PCR（RT-qPCR）

根据相关协议，使用Trizol试剂分离总RNA。cDNA的合成使用ReverTra Ace qPCR RT Master Mix进行。使用SYBR Green Real-Time PCR Master Mixs在CFX-96TM上进行qPCR反应。检测每个基因的量并通过GAPDH量进行归一化。

 

2.14转录因子网络分析

Cheng课题组从之前的两项研究中收集719个据报道的HIF1A转录因子靶标。为检查HIF1A转录因子在LUAD中的潜在功能，研究人员利用LUAD患者肿瘤样品中RNA-seq读数计数数据与DESeq2数据中TCGA匹配正常样品进行比较，计算719个基因的差异表达。研究人员使用调整后的P-value<0.05来定义差异表达基因。

 

2.15代谢物丰度与基因表达的相关性

Cheng课题组从先前的研究中收集了在约70个NSCLC细胞株中测试的丝氨酸或甘氨酸的丰度。接下来，研究人员从CCLE数据库中收集了NSCLC细胞株的全基因组基因表达谱。代谢物丰度与基因表达水平之间相关性通过Pearson相关系数计算，P-value用R平台F-statistics计算(v3.01, http://www.r-project.org/)。

 

2.16组织特定子网分析

Cheng课题组从GTEx V6（https://gtexportal.org/home/）下载了32个组织的RNA-seq数据（RPKM值）。对于每个组织（例如，肺），研究人员将超过80％样本中RPKM≥1的基因视为组织表达基因，其余基因视为组织未表达的基因。为了量化组织表达基因i在组织t中的表达意义，研究人员计算了在所有考虑的组织中，该基因的平均表达和标准差，基因表达在组织t中的重要性定义为。对于LUAD，研究人员通过将肺与GTEx中其他31种不同组织全基因组表达谱进行比较，建立了肺特异性蛋白质相互作用网络。

 

2.17生存分析

Cheng课题组下载了226种人类原发性I–II期肺腺癌的微阵列数据和生存概况。根据归一化表达水平，将患者分为前50％低表达（蓝色）组和前50％高表达（红色）组。使用GraphPad Prism 7软件中log-rank test确定Kaplan–Meier生存分析的P-value。

 

2.18统计分析

研究中的数据是从至少三个独立的实验中获得，不同实验组中所有数据均表示平均值±平均值的标准误差。使用Student’s t-test和P-values分析了两组之间的差异。当P<0.05时，认为具有统计学意义，本研究使用GraphPad Prism 7软件进行统计分析。

3. 结果

 

3.1  基于GPSnet识别癌症类型特异性疾病模块

Cheng课题组通过将患者全外显子组测序和RNA-seq数据反复整合到人类蛋白相互作用组中，计算确定所有15种癌症类型的最终癌症类型特异性疾病模块,如图1所示。

图1.基于GPSnet网络发现硅基药物再利用流程

 

3.2  基于GPSet识别新药物基因组学生物标志

为了研究GPSnet网络潜在药物基因组学应用，Cheng课题组从癌症药物敏感性基因组学(GDSC)数据库17下载了1065个细胞株的健壮多阵列(RMA)基因表达谱和药物反应数据。Cheng课题组根据GPSnet网络识别的疾病模块的RMA基因表达谱为特征，构建回归模型，对三种特定的癌症（BRCA, LUAD,SKCM），进行药物反应预测。基于多种因素的专业知识，集中研究了这三种癌症类型的7种药物，其AUC性能指标如图2所示。

图2.基于gpsnet网络预测的疾病模块药物基因组学验

 

3.3  基于GPSnet疾病模块的药物高可吸收性

相较于传统的统计学方法SMG,基于GPSnet网络预测模基因产物更有可能被批准的药物或临床试验药物的靶向目标，图3.在两种特殊的癌症中基GPSnet识别的疾病模块药物目标富集分析。

图3.在两种特殊的癌症中基于GPSnet识别的疾病模块药物目标富集分

 

3.4  基于GPSnet发现已批准药物的新适应症状

研究使用了网络距离法与基因集合富集分析法验证在已批准药物中识别的新的适应症，结果如4所示

图4.基于网络的有机硅药物再利用与实验验

 

3.5  Ouabain在NSCLC细胞中抑制HIF1α/LEO1路径

Cheng课题组通过在肺特异性人类蛋白质相互作用组的网络分析，研究了药物ouabain在NSCLC中的作用机制。结果如图5所示

图5.网络的实验验证药物ouabain在NSCLC的作用机制

4. 结论

本研究提供了一种基于网络的药物再使用方法，并通过实验验证了ouabain与 niclosamide作为治疗NSCLC的潜在药物是有效的，并且Cheng课题组通过实验发现药物ouabain的作用目标是NSCLC细胞中蛋白质HIF1α/LEO1-mediated 细胞代谢路径，GPSnet是一种基于网络的方法，集成了大规模患者DNA/RNA-seq数据与人类蛋白质相互组数据，在癌症和其他疾病领域有望能够加速药物靶向识别与药物开发的过程。

 

Data availability

http://ccsb.dana-farber.org/interactome-data.html.

 

Code availability

thttps://github.com/ChengF-Lab/GPSnet.

参考资料

Cheng F, Lu W, Liu C, et al. A genome-wide positioning systems network algorithm for in silico drug repurposing[J]. Nature communications, 2019, 10(1): 1-14.