CN111524554B - 基于lincs-l1000扰动信号的细胞活性预测方法 - Google Patents

基于lincs-l1000扰动信号的细胞活性预测方法 Download PDF

Info

Publication number
CN111524554B
CN111524554B CN202010331009.6A CN202010331009A CN111524554B CN 111524554 B CN111524554 B CN 111524554B CN 202010331009 A CN202010331009 A CN 202010331009A CN 111524554 B CN111524554 B CN 111524554B
Authority
CN
China
Prior art keywords
drug
cell activity
data
genes
cell
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010331009.6A
Other languages
English (en)
Other versions
CN111524554A (zh
Inventor
陈明
秦玉芳
陆家兴
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shanghai Ocean University
Original Assignee
Shanghai Ocean University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shanghai Ocean University filed Critical Shanghai Ocean University
Priority to CN202010331009.6A priority Critical patent/CN111524554B/zh
Publication of CN111524554A publication Critical patent/CN111524554A/zh
Application granted granted Critical
Publication of CN111524554B publication Critical patent/CN111524554B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding
    • G16B40/20—Supervised data analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N20/00—Machine learning
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N20/00—Machine learning
    • G06N20/20—Ensemble learning
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N5/00—Computing arrangements using knowledge-based models
    • G06N5/01—Dynamic search techniques; Heuristics; Dynamic trees; Branch-and-bound
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B50/00—ICT programming tools or database systems specially adapted for bioinformatics
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Medical Informatics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Software Systems (AREA)
  • Data Mining & Analysis (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Computation (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Computing Systems (AREA)
  • General Physics & Mathematics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Biophysics (AREA)
  • Mathematical Physics (AREA)
  • General Engineering & Computer Science (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Databases & Information Systems (AREA)
  • Bioethics (AREA)
  • Computational Linguistics (AREA)
  • Epidemiology (AREA)
  • Public Health (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Genetics & Genomics (AREA)
  • Molecular Biology (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本发明公开了基于LINCS‑L1000扰动信号的细胞活性预测方法，包括LINCS‑L1000、Achilles和CTRP的数据集预处理、建立模型、基因集富集分析、CTRP‑L1000与Achilles‑L1000跨数据集交互式验证和抗癌药物诱导下的细胞活性预测的步骤。本发明基于随机森林与弹性网的特征重要性加权融合算法模型，并结合FEBPSO‑XGBoost算法将基因组特征与癌细胞系中的药物敏感性相关联，建立药物反应的模型预测抗癌药物诱导下的细胞活性，克服现有的细胞活性预测方法存在预测精度过低和拟合较差的问题，用于新细胞系的细胞活性预测可为指导临床用药提供理论指导。

Description

基于LINCS-L1000扰动信号的细胞活性预测方法
技术领域
本发明属于生物信息学和机器学习技术领域，具体涉及基于LINCS-L1000扰动信号的细胞活性预测方法。
背景技术
癌症是全球第二大死亡原因，恶性肿瘤的产生基于多种因素，近年来，细胞凋亡与肿瘤的关系成为研究的热点之一，其中细胞肿瘤的相关基因，如癌基因凋亡抑制基因、细胞凋亡基因等的异常波动是引起某些生物学特征变化的重要原因。当基因调控出现异常时，可使细胞凋亡不足或增多，从而引起多样性疾病，其过程是在多基因严格调控下进行的，比如细胞凋亡功能的抑制将导致肿瘤的发生及免疫功能的异常，而细胞凋亡功能亢进与组织损伤和器官功能衰竭的发病密切相关，各种外界因素和遗传因素可能导致细胞失去活性。有研究表明，肿瘤细胞凋亡抑制因子活性的增强，能够抵御细胞凋亡信号通路的激活，细胞凋亡比例变化和细胞增殖的异常行为与化合物浓度时间和特征高度相关，以及在胚胎发育、组织器官塑形及衰老和病态细胞中起重要作用，是肿瘤形成与发展的关键因素之一。对于大多数刚进入临床试验阶段的化合物反应疗效，我们仍然缺乏正确选择针对特定肿瘤的靶向药物的方法，从而阻碍了癌症药物的发展，应用细胞活性机制通过化合物作用激发和恢复肿瘤细胞发生凋亡，是了解肿瘤发生机制和研发抗癌药物的有效途径。
为了寻找有效的抗癌药物，通过细胞信号通路中的功能及其相互作用对癌症特异性的改变，构建一个细胞信号模型来解释这些数据并做出准确的预测仍然是一个挑战，随着越来越多的癌症组学数据的出现，采用机器学习模型来模拟和发现化学物质对细胞机制的潜在生化过程变得越来越重要。
由于细胞扰动信号与细胞生存能力紧密相关，近年来，在药物敏感性和抗癌药物反应预测的研究中，利用现代机器学习算法开发基于每个肿瘤的分子谱测量的药物反应预测因子，通过从不同高覆盖分子数据预测细表型，有效的识别已知的药物基因组预测生物标志物，采用化合物控制凋亡通路中的关键蛋白或酶类的表达或功能，诱导异常细胞凋亡，从大规模药物筛选数据库预测细胞活性已成为个性化医疗的最佳方法。在最近一些关于药物敏感性和抗癌药物反应预测的研究中，机器学习模型被得到广泛应用，一种常见的方法是考虑使用实验测量的基因组特征(RNA表达、蛋白质表达、甲基化、SNP等)和对不同药物与细胞系之间的反应作为训练集，并针对每种单独药物设计一个或多个监督预测模型。多任务学习方法由于在训练时采用药物关联方法而备受青睐，依赖于基于内核的降维和多任务学习的贝叶斯(KBMTL)方法表现出显著的药物反应预测性能。
目前，现有的细胞活性研究工作大多是针对单个数据集建模，把细胞系的基因组特征作为输入变量构建回归模型，模型预测的准确度始终不能令人满意，我们发现在使用传统的回归模型预测时，由于药物浓度与时间并未考虑到预测分析中，使用传统的预测方法还未真正捕捉到细胞活性反应机理的实质。在使用线性回归方法预测差异表达信号对药物敏感性的影响中，发现药物引起的细胞凋亡与剂量和时间有关，但其由于只提取了差异表达基因与细胞活性之间的线性特征，无法同时提取差异表达基因与细胞活性之间的线性与非线性特征，存在预测精度过低和拟合较差的问题。
发明内容
本发明的目的在于针对现有的技术不足，提出一种基于LINCS-L1000扰动信号的细胞活性预测方法，包括以下步骤：
S1：数据预处理
对包括LINCS-L1000扰动转录组学信号、癌症治疗反应门户(CTRP)中药物治疗后的细胞活性信息、癌症依赖性图谱数据库(Achilles)中shRNA治疗前后的效应改变量的数据集进行匹配和筛选，获取差异表达基因的转录水平表达量与细胞表型信息；
S2：建立模型
建立随机森林与弹性网的特征重要性加权融合算法模型，以准确性函数最大化为标准，根据步骤S1预处理后的信息筛选得到的差异表达基因DEGs选择关键基因，并通过FEBPSO-XGBoost算法将关键基因的特征与癌细胞系中的药物治疗反应相关联，建立药物反应的模型预测细胞活性；
S3：基因集富集分析
包括对所述关键基因分别进行GO富集分析和KEGG通路分析的步骤；
S4：CTRP-L1000与Achilles-L1000跨数据集交互式验证
包括Achilles-L1000数据训练的系列模型预测CTRP-L1000数据中的细胞活性和CTRP-L1000数据训练的系列模型预测Achilles-L1000数据中的细胞活性的步骤，上述跨数据集交互式验证所述模型的有效性；
S5：抗癌药物诱导下的细胞活性预测
基于新细胞系的基因组特征，经步骤S4验证有效的模型预测抗癌药物诱导下的细胞活性。
根据本发明的一些实施方式，步骤S1中，所述数据预处理的过程包括分步骤：
S101：获取两阶段的LINCS-L1000扰动谱数据，并进行合并，得到在多种扰动状态下全基因组的基因表达情况；
S102：按照细胞系、化合物和浓度信息，将LINCS-L1000中的扰动转录组学信号与CTRP中药物治疗后的细胞活性信息、Achilles中shRNA治疗前后的效应改变量进行匹配和筛选，得到在多点浓度和多个时间点下的差异表达基因的转录水平表达量与细胞表型信息的数据集；
S103：步骤S102所得数据集按照不同的扰动时间和扰动类型分为3小时、6小时、24小时、96小时、120小时和144小时六大数据子集。
根据本发明的一些实施方式，步骤S2中，所述细胞活性预测分析的过程包括：通过FEBPSO-XGBoost算法评估多个不同扰动时间和不同药物浓度作用下关键基因的表达水平，预测不同细胞系在药物或者shRNA治疗后的细胞活性。
根据本发明的一些实施方式，步骤S2中，包括分步骤：
S201：建立随机森林与弹性网的特征重要性加权融合算法模型，分别依据所述关键基因在随机森林和弹性网中的贡献度进行排序，并通过以下公式(1)进行加权求和，获取的加权平均值与所需要选取的最佳基因数再次按基因贡献度依次排列：
Figure BDA0002464954760000031
其中，RFPearson和ENPearson均为随机森林与弹性网算法在验证集PertDT-V上的Pearson相关系数，该值越高说明计算所得的预测值与实验测量值之间的相关程度越高；
Figure BDA0002464954760000042
为差异表达基因DEGs在随机森林算法中的特征重要性排序或者选取的基因个数；
Figure BDA0002464954760000043
为差异表达基因DEGs在弹性网算法中的特征重要性排序或者选取的基因个数；
S202：在每棵决策树的叶子节点上有一个预测分数，通过如公式2所示的XGBoost算法多次迭代构建若干个弱评估器，所有决策树的预测分数之和为预测的细胞活性：
Figure BDA0002464954760000041
其中，fk(samplei[DEGs])表示第i个样本sample在选取的差异表达基因集DEGs上，在第k棵决策树上的预测分数，K表示决策树的数量；
S203：基于转录组学差异表达基因的表达水平，通过FEBPSO算法对参数进行寻优，得到每个模型的最优参数组合，依据所述关键基因与药物治疗反应之间的基因调控关系，通过上述最优参数组合和XGBoost算法对细胞活性预测。
根据本发明的一些实施方式，步骤S3中，所述GO富集分析与所述KEGG通路分析均与细胞凋亡过程相关。
根据本发明的一些实施方式，还包括通过癌细胞系百科全书CCLE和NCI60数据集中的数据验证模型的有效性，具体为：
a)在NCI60数据集中，使用GI50值作为药物敏感性的评价指标，当药效在50％生长抑制浓度范围内，对应药物敏感性评价指标GI50值，当药效在50％生长抑制浓度范围内无效时，记为最高浓度值，定义药物浓度差异变量(如式3所示)：
Δdrug_conc(dr，cl)＝drug_sens(dr，cl)-test_max_conc(dr，cl) (3)
其中，drug_sens(dr，cl)为使用药物dr治疗细胞系cl时记录的药物敏感性值GI50，test_max_conc(dr，cl)为使用药物dr治疗细胞系cl时所使用药物的最大测试浓度量；
b)在CCLE数据集中，采用活性面积作为药物敏感性的评价指标，并进行二值化处理。
与现有技术相比，本发明的有益效果是：
不同细胞在特定化合物作用下具有不同的扰动信号，为了评估药物敏感性预测在临床前的可行性，大量的药物诱导基因表达数据得以公开使用。通过这些数据，将扰动信号与细胞表型研究联系起来，建立机器学习模型，分析不同药物在不同疾病中的扰动关系，揭示隐藏在表型之下的生物学分子作用机理，预测个体患者的临床效果。由于癌细胞系对化学物质的敏感性不同，本发明通过对LINCS-L1000、Achilles和CTRP三大数据集匹配和筛选，采取多浓度和多时间点测量细胞全基因组的转录水平表达量，并基于随机森林与弹性网加权融合算法选择关键基因，利用特征基因的差异表达水平，结合FEBPSO-XGBoost算法将基因组特征与癌细胞系中的药物敏感性相关联，建立药物反应的预测模型，预测药物诱导下的细胞活性，克服现有预测方法存在预测精度过低和拟合较差的问题。
附图说明
图1为本发明细胞活性预测方法的工艺流程图。
图2为CTRP-L1000-24h数据集中的差异表达基因富集分析图。
图3为Achilles-L1000-96h数据集中的差异表达基因富集分析图。
图4为跨数据集验证结果图，使用Achilles-L1000数据训练的系列模型预测CTRP-L1000数据中的细胞活性，反之亦然。
图5为模型在LINCS-L1000-NCI60-24h上验证的ROC曲线(A)与PR曲线(B)。
图6为模型在LINCS-L1000-CTRP-NCI60-24h上验证的ROC曲线(A)与PR曲线(B)。
图7为多种药物在细胞系MCF7和PC3上的最低细胞活性值(A)和多种药物在细胞系A375和HT29上的最低细胞活性值(B)。
具体实施方式
本发明数据来源于互联网公开的数据集：扰动转录组学信号，包括LINCS-L1000-PhaseI和LINCS-L1000-PhaseII，选自LINCS-L1000，从GEO数据库(https://www.ncbi.nlm.nih.gov/geo/)下载获得；CTRP可从https://ocg.cancer.gov/programs/ctd2/data-portal下载获得；癌症依赖性图谱数据库Achilles可从https://portals.broadinstitute.org/achilles下载获得；NCI-60和CCLE分别从https://dtp.cancer.gov/discovery_development/nci-60和https://portals.broadinstitute.org/ccle下载获得。
如图1所示，本发明提供了一种基于LINCS-L1000扰动信号的细胞活性预测方法，先在扰动转录组学信号LINCS-L1000与癌症治疗反应门户CTRP数据集上完成差异表达基因的选择与细胞活性的预测分析，并进行模型结果的评估；再在癌症依赖性图谱数据库Achilles上进行独立集测试，得出模型在跨数据集上的性能(此处只列出CTRP-L1000系列模型在跨数据集Achilles-L1000上的测试过程，反之亦然)；然后根据癌细胞系百科全书CCLE数据集中的活性区域面积值与NCI-60数据集中的药物敏感性指标，完成模型的验证；最后用训练完成的模型预测未知抗癌药物治疗特定细胞系中产生的细胞活性。具体包括：
步骤一：数据预处理
对包括LINCS-L1000扰动转录组学信号、癌症治疗反应门户(CTRP)中药物治疗后的细胞活性信息、癌症依赖性图谱数据库(Achilles)中shRNA治疗前后的效应改变量的数据集进行匹配和筛选，获取差异表达基因的转录水平表达量与细胞表型信息；具体可以为：
S101：获取两阶段的LINCS-L1000扰动谱数据，并进行合并，得到在多种扰动状态下全基因组的基因表达情况；
S102：按照细胞系、化合物和浓度信息，将LINCS-L1000中的扰动转录组学信号与CTRP中药物治疗后的细胞活性信息、Achilles中shRNA治疗前后的效应改变量进行匹配和筛选，得到在多点浓度和多个时间点下的差异表达基因的转录水平表达量与细胞表型信息的数据集；
S103：步骤S102所得数据集按照不同的扰动时间和扰动类型分为3小时、6小时、24小时、96小时、120小时和144小时六大数据子集。
步骤二：建模和细胞活性预测分析
建立随机森林与弹性网的特征重要性加权融合算法模型，以准确性函数最大化为标准，根据步骤一预处理后的信息筛选得到的差异表达基因DEGs选择关键基因，并通过FEBPSO-XGBoost算法依据所述关键基因与药物治疗反应之间的基因调控关系对细胞活性预测分析；具体可以为：
S201：建立随机森林与弹性网的特征重要性加权融合算法模型，分别依据所述关键基因在随机森林和弹性网中的贡献度进行排序，并通过以下公式(1)进行加权求和，获取的加权平均值与所需要选取的最佳基因数再次按基因贡献度依次排列：
Figure BDA0002464954760000071
其中，RFPearson和ENPearson均为随机森林与弹性网算法在验证集PertDT-V上的Pearson相关系数，该值越高说明计算所得的预测值与实验测量值之间的相关程度越高；
Figure BDA0002464954760000073
为差异表达基因DEGs在随机森林算法中的特征重要性排序或者选取的基因个数；
Figure BDA0002464954760000074
为差异表达基因DEGs在弹性网算法中的特征重要性排序或者选取的基因个数。
S202：在每棵决策树的叶子节点上有一个预测分数，通过XGBoost算法(如公式2所示)多次迭代逐一构建多个弱评估器，所有树的预测分数之和即为预测的细胞活性：
Figure BDA0002464954760000072
其中，fk(samplei[DEGs])表示第i个样本sample在选取的差异表达基因集DEGs上，在第k棵决策树上的预测分数，K表示决策树的数量；
S203：基于转录组学差异表达基因的表达水平，通过FEBPSO算法对参数进行寻优，得到每个模型的最优参数组合，通过FEBPSO-XGBoost算法结合评估多个不同扰动时间和不同药物浓度作用下关键基因的表达水平，预测不同细胞系在药物或者shRNA治疗后的细胞活性。
步骤三：基因集富集分析
对步骤二中选择的关键基因分别进行GO富集分析和KEGG通路分析，两者均与细胞凋亡过程相关。
步骤四：CTRP-L1000与Achilles-L1000跨数据集交互式验证
包括：Achilles-L1000数据训练的系列模型预测CTRP-L1000数据中的细胞活性，CTRP-L1000数据训练的系列模型预测Achilles-L1000数据中的细胞活性。
步骤五：抗癌药物诱导下的细胞活性预测
基于新细胞系的基因组特征，用步骤S1至S4构建的WRFEN-XGBoost模型预测抗癌药物诱导下的细胞活性。
根据一些实施方式，还包括通过癌细胞系百科全书CCLE和NCI60数据集中的数据验证模型的有效性，具体为：
a)在NCI60数据集中，使用GI50值作为药物敏感性的评价指标，当药效在50％生长抑制浓度范围内，对应药物敏感性评价指标GI50值，当药效在50％生长抑制浓度范围内无效时，记为最高浓度值，定义药物浓度差异变量：
Δdrug_conc(dr，cl)＝drug_sens(dr，cl)-test_max_conc(dr，cl)
其中，drug_sens(dr，cl)为使用药物dr治疗细胞系cl时记录的药物敏感性值GI50，test_max_conc(dr，cl)为使用药物dr治疗细胞系cl时所使用药物的最大测试浓度量；
b)在CCLE数据集中，采用活性面积作为药物敏感性的评价指标，并进行二值化处理。
以下通过具体实施例对本发明进一步解释说明。
实施例1
按照上述数据来源，获取LINCS-L1000扰动转录组学信号、癌症治疗反应门户CTRP中化合物治疗后的细胞活性信息、癌症依赖性图谱数据库Achilles项目中shRNA治疗前后的效应改变量和癌细胞系百科全书CCLE和NCI60数据集中的药物敏感性的数据集，并对其进行预处理，具体为：
(1)LINCS-L1000与CTRP关联。合并两阶段LINCS-L1000扰动谱数据，基于相同的细胞系和Broad研究所提供的药物识别号，与CTRP中药物治疗后的细胞活性数据相关联，根据如下公式4控制浓度差异，针对同一种信号，在相同浓度条件、不同实验批次中测得的细胞活性值，取细胞活性平均值。最后将整理完成的数据按照不同的扰动时间进一步划分成3小时、6小时和24小时三个数据子集(CTRP-L1000-3h、CTRP-L1000-6h和CTRP-L1000-24h)进行学习：
doseDifference＝|log10(CTRPdose)-log10(LINCSdose)|≤0.2 (4)
其中，CTRPdose为CTRP数据集中癌症治疗药物对应的浓度值，LINCSdose为LINCS-L1000数据集中不同扰动信号对应的浓度值。
(2)LINCS-L1000与Achilles关联。将合并完成后的两阶段LINCS-L1000扰动谱数据，依据相同的细胞系和shRNA治疗剂，与Achilles项目进行关联，采用shRNA治疗前后的效应改变量大小作为细胞表型信息。
实施例2
构建WRFEN-XGBoost模型，具体步骤如下：
(1)建立随机森林与弹性网的特征重要性加权融合算法模型，分别依据关键基因在随机森林和弹性网中的贡献度进行排序，并使用公式进行加权求和，对获取的加权平均值与所需要选取的最佳基因数，再次按基因贡献度依次排列：
Figure BDA0002464954760000091
其中，RFPearson与ENPearson为使用随机森林与弹性网算法在验证集PertDT-V上的Pearson相关系数，该值越高说明计算所得的预测值与实验测量值之间的相关程度越高，
Figure BDA0002464954760000093
为差异表达基因DEGs在随机森林算法中的特征重要性排序或者选取的基因个数，/>
Figure BDA0002464954760000094
为差异表达基因DEGs在弹性网算法中的特征重要性排序或者选取的基因个数。
(2)使用XGBoost算法对细胞活性进行预测，每棵决策树的叶子节点上设一个预测分数，通过多次迭代逐一构建多个弱评估器，细胞活性预测结果是所有树的预测分数之和：
Figure BDA0002464954760000092
其中，fk(samplei[DEGs])表示第i个样本sample在选取的差异表达基因集DEGs上，在第k棵决策树上的预测分数，K表示决策树的数量。
(3)使用FEBPSO算法对参数进行寻优，解决调节参数过多导致训练时间过长的问题，得到每个模型的最优参数组合，并使用最优参数和XGBoost预测模型进行训练与测试。
以下通过与其他现有方法进行分析比较评估该模型的有效性，分别使用预测值与真实值的Pearson相关性、R2和均方误差进行衡量模型的预测性能，如表1所示。以CTRP-L1000-24h(S1)数据集为例，分别以0.8321、0.6922、0.025优于其他算法。具体而言，相比于PCA-Lasso、PCA-SVR、FTest-RF和MI-KNN算法，在Pearson相关系数方面，分别提高了5.50％、5.33％、4.77％和3.32％；在R2方面，分别提高了11.45％、11.45％、9.80％和7.92％；在均方误差方面，分别下降了19.35％、19.35％、16.67％和16.67％。结果表明，本发明的预测结果得到进一步的提升，其他模型的评估结果如表1、2、3所示。
表1本发明与其他算法的比较(Pearson相关性)
所用数据集 本发明 PCA-Lasso PCA-SVR Ftest-RF MI-KNN
CTRP-L1000-3h(浓度) 0.8392 0.7291 0.7166 0.7392 0.7791
CTRP-L1000-6h(浓度) 0.6770 0.5985 0.5962 0.6343 0.6105
CTRP-L1000-24h(浓度) 0.8859 0.8208 0.8176 0.8534 0.8706
CTRP-L1000-3h 0.7705 0.7210 0.7242 0.6488 0.7189
CTRP-L1000-6h 0.6239 0.5507 0.5545 0.5588 0.4991
CTRP-L1000-24h 0.8321 0.7887 0.7900 0.7942 0.8054
Achilles-L1000-96h 0.5893 0.5392 0.5266 0.5215 0.4334
Achilles-L1000-120h 0.5275 0.5036 0.5000 0.4800 0.3328
Achilles-L1000-144h 0.5348 0.5176 0.4975 0.4699 0.4066
表2本发明与其他算法的比较(R2)
所用数据集 本发明 PCA-Lasso PCA-SVR Ftest-RF MI-KNN
CTRP-L1000-3h(浓度) 0.6884 0.5179 0.5040 0.5362 0.5781
CTRP-L1000-6h(浓度) 0.4578 0.3582 0.3423 0.3969 0.3148
CTRP-L1000-24h(浓度) 0.7847 0.6736 0.6675 0.7262 0.7543
CTRP-L1000-3h 0.5803 0.5137 0.5153 0.4191 0.4828
CTRP-L1000-6h 0.3884 0.3012 0.2927 0.3098 0.1566
CTRP-L1000-24h 0.6922 0.6211 0.6211 0.6304 0.6414
Achilles-L1000-96h 0.3468 0.2905 0.2682 0.2684 0.1068
Achilles-L1000-120h 0.2782 0.2496 0.2379 0.2236 0.0240
Achilles-L1000-144h 0.2848 0.2657 0.2350 0.2176 0.1133
表3本发明与其他算法的比较(均方误差)
所用数据集 本发明 PCA-Lasso PCA-SVR Ftest-RF MI-KNN
CTRP-L1000-3h(浓度) 0.021 0.033 0.034 0.032 0.029
CTRP-L1000-6h(浓度) 0.054 0.064 0.066 0.060 0.069
CTRP-L1000-24h(浓度) 0.018 0.027 0.027 0.023 0.020
CTRP-L1000-3h 0.029 0.033 0.033 0.040 0.035
CTRP-L1000-6h 0.064 0.073 0.074 0.072 0.088
CTRP-L1000-24h 0.025 0.031 0.031 0.030 0.030
Achilles-L1000-96h 1.477 1.604 1.655 1.654 2.019
Achilles-L1000-120h 1.391 1.446 1.468 1.496 1.880
Achilles-L1000-144h 1.307 1.342 1.398 1.429 1.620
效果说明
(1)关键基因选择对细胞活性预测的影响
药物反应下的基因表达变化将直接影响预测性能，使用WRFEN算法进行关键基因的选择，达到最大化预测准确性的函数来选择出最佳基因子集，将扰动谱与细胞活性数据空间PertDT划分为70％的训练集PertDT-S、15％的验证集PertDT-V与15％的测试集PertDT-T，分别依据关键基因在随机森林和弹性网中的贡献度进行排序，如表4列出了每个数据子集中所选择的特征基因的数量和15个关键特征基因的贡献度排名。
表4算法WRFEN筛选出的前15个特征基因
Figure BDA0002464954760000111
(2)基因集富集分析
以CTRP-L1000-24h与Achilles-L1000-96h两大数据子集为例，对提取的特征基因进行GO富集分析与KEGG通路分析，从图2和3中可以发现两者均与细胞凋亡过程密切相关，其中最显著富集的通路R-HAS-1640170和GO:0007346是对细胞周期与细胞凋亡的调控，这也证实了本发明中使用药物或者shRNA治疗后差异表达基因构成了细胞凋亡的途径。
(3)基于FEBPSO-XGBoost预测药物诱导下的细胞活性
药物在不同时间和不同浓度范围内具有的活性不尽相同，考虑模型可能会涵盖更大的动态范围，为使模型具有更高的准确性和可行性，将训练数据集按照不同时间和不同扰动类型划分为CTRP-L1000-3h、CTRP-L1000-6h、CTRP-L1000-24h、Achilles-L1000-96h、Achilles-L1000-120h、Achilles-L1000-144h六大数据子集，分别进行WRFEN核心基因选择和FEBPSO-XGBoost算法结合，形成一个完整的预测模型，以解释细胞系在特定药物与浓度作用下观察到的完整的细胞凋亡水平，在理论上提供了一种合理且易于实现的方法，通过测试集评估了算法的性能，使模型的预测能力得到显著提升，如表5所示。
表5 WRFEN-XGBoost算法在不同数据子集上的评估结果
所用数据子集 Pearson相关性 R2 均方误差
CTRP-L1000-3h 0.7705 0.5803 0.029
CTRP-L1000-6h 0.6239 0.3884 0.064
CTRP-L1000-24h 0.8321 0.6922 0.025
Achilles-L1000-96h 0.5893 0.3468 1.477
Achilles-L1000-120h 0.5275 0.2782 1.391
Achilles-L1000-144h 0.5348 0.2848 1.307
(4)CTRP-L1000与Achilles-L1000跨数据集验证
如图4所示，在CTRP-L1000数据集中，以24小时的扰动时间最佳，在Achilles-L1000数据集中，以96小时的扰动时间最佳，在独立集验证方面，CTRP-L1000-24h数据集在CTRP-L1000-6h模型、CTRP-L1000-24h模型和Achilles-L1000-96h模型上分别以0.7416、0.8321和0.7319的Pearson相关系数优于其他模型，进一步证实了药物在较长的扰动时间后能够取得优良的预测性能。
(5)基于NCI60与CCLE的药物敏感性推断
使用GI50值作为药物敏感性评价指标，当药效在50％生长抑制浓度取值范围内，对应药物敏感性评价指标GI50值，当药效在50％生长抑制浓度范围内无效时，记为最高浓度值。定义药物浓度差异变量，其计算方法如公式(4)所示，当药物浓度差异变量的值小于零时，代表该药物为有效药物，反之为无效药物：
Δdrug_conc(dr，cl)＝drug_sens(dr，cl)-test_max_conc(dr，cl) (4)
其中，Δdrug_conc(dr,cl)为使用药物dr治疗细胞系cl时产生的药物浓度差异变量，drug_sens(dr,cl)为使用药物dr治疗细胞系cl时记录的药物敏感性值GI50，test_max_conc(dr,cl)为使用药物dr治疗细胞系cl时所使用药物的最大测试浓度量。
如图5所示，使用ROC曲线与PR曲线来衡量算法在评判药物有效性方面所做的贡献，图(A)所示的ROC曲线可以发现，在LINCS-L1000-NCI60-24h数据集的预测过程中，以Achilles-L1000-96h模型所做的预测最为准确。图(B)所示的准确率-召回率评估曲线中，Achilles-L1000-96h模型仍以曲线下的面积AUC＝0.94超越其他模型，证实了Achilles-L1000-96h模型在LINCS-L1000-NCI60-24h数据集的预测过程中是有效的，可进一步用于后续其他药物的有效性检测。
如图6所示，将LINCS-L1000、CTRP数据和NCI60数据按照前面所叙述的匹配方法进行关联，仍然使用扰动时间最长的药物，即筛选扰动时间为24小时对应的药物，记匹配完成的数据为LINCS-L1000-CTRP-NCI60-24h。使用Achilles-L1000-96h模型对药物有效性进行预测时，ROC曲线下的面积AUC值为0.78，PR曲线下的面积AUC为0.98。同样使用CCLE中的药物敏感性数据对预测的细胞活性值进行评估时，Achilles-L1000-96h模型在跨数据集验证上也表现出了卓越的性能。
(6)药物诱导下的细胞活性预测
通过新细胞系的基因组特征对其进行细胞活性预测分析。如图7所示，使用散点图绘制多种药物在不同细胞系中的细胞活性最低预测值，Tamoxifen在乳腺癌MCF7中的抗雌激素活性使存活肿瘤区域的内皮密度发生了显著下降，导致肿瘤消退；在治疗PC3前列腺癌中由于细胞增殖效率的降低，Mitoxantrone和Doxorubicin显示出坏死作用；Evodiamine可诱导人类黑色素瘤A375细胞死亡；Disulfiram可介导黑色素瘤A375的凋亡，同时对正常细胞产生最小毒性。
综上所述，本发明的细胞活性预测方法可以广泛用于新细胞系细胞活性预测当中，为指导临床用药提供理论指导。
上述对实施例的描述是为了便于该技术领域的普通技术人员能理解和使用本发明。熟悉本领域技术人员显然可以容易的对这些实施例作出各种修改，并把在此说明的一般原理应用到其他实施例中，而不必经过创造性的劳动。因此，本发明不限于上述实施例。本领域技术人员根据本发明的原理，不脱离本发明的范畴所做出的改进和修改都应该在本发明的保护范围之内。

Claims (3)

1.基于LINCS-L1000扰动信号的细胞活性预测方法，其特征在于，所述方法包括以下步骤：
S1：数据预处理：
对包括LINCS-L1000扰动转录组学信号、癌症治疗反应门户(CTRP)中药物治疗后的细胞活性信息、癌症依赖性图谱数据库(Achilles)中shRNA治疗前后的效应改变量的数据集进行匹配和筛选，获取差异表达基因的转录水平表达量与细胞表型信息；
S2：建立模型：
建立随机森林与弹性网的特征重要性加权融合算法模型，以准确性函数最大化为标准，根据步骤S1预处理后的信息筛选得到的差异表达基因DEGs选择关键基因，并通过FEBPSO-XGBoost算法将关键基因的特征与癌细胞系中的药物治疗反应相关联，建立药物反应的模型预测细胞活性；
S3：基因集富集分析：包括对所述关键基因分别进行GO富集分析和KEGG通路分析的步骤；
S4：CTRP-L1000与Achilles-L1000跨数据集交互式验证：
包括Achilles-L1000数据训练的系列模型预测CTRP-L1000数据中的细胞活性和CTRP-L1000数据训练的系列模型预测Achilles-L1000数据中的细胞活性的步骤，上述跨数据集交互式验证所述模型的有效性；
S5：抗癌药物诱导下的细胞活性预测：
基于新细胞系的基因组特征，经步骤S4验证有效的模型预测抗癌药物诱导下的细胞活性；
步骤S2中，通过FEBPSO-XGBoost算法评估多个不同扰动时间和不同药物浓度作用下关键基因的表达水平，预测不同细胞系在药物或者shRNA治疗后的细胞活性；
步骤S2中，包括分步骤：
S201：建立随机森林与弹性网的特征重要性加权融合算法模型，分别依据所述关键基因在随机森林和弹性网中的贡献度进行排序，并通过以下公式(1)进行加权求和，获取的加权平均值与所需要选取的最佳基因数再次按基因贡献度依次排列：
Figure FDA0003960628160000021
其中，RFPearson和ENPearson均为随机森林与弹性网算法在验证集PertDT-V上的Pearson相关系数，该值越高说明计算所得的预测值与实验测量值之间的相关程度越高；
Figure FDA0003960628160000022
为差异表达基因DEGs在随机森林算法中的特征重要性排序或者选取的基因个数；
Figure FDA0003960628160000023
为差异表达基因DEGs在弹性网算法中的特征重要性排序或者选取的基因个数；
S202：在每棵决策树的叶子节点上有一个预测分数，通过如公式2所示的XGBoost算法多次迭代构建若干个弱评估器，所有决策树的预测分数之和为预测的细胞活性：
Figure FDA0003960628160000024
其中，fk(samplei[DEGs])表示第i个样本sample在选取的差异表达基因集DEGs上，在第k棵决策树上的预测分数，K表示决策树的数量；
S203：基于转录组学差异表达基因的表达水平，通过FEBPSO算法对参数进行寻优，得到每个模型的最优参数组合，依据所述关键基因与药物治疗反应之间的基因调控关系，通过上述最优参数组合和XGBoost算法对细胞活性预测；
步骤S3中，所述GO富集分析与所述KEGG通路分析均与细胞凋亡过程相关。
2.根据权利要求1所述的细胞活性预测方法，其特征在于，所述的数据预处理过程包括分步骤：
S101：获取两阶段的LINCS-L1000扰动谱数据，合并，得到在多种扰动状态下全基因组的基因表达情况；
S102：按照细胞系、化合物和浓度信息，将LINCS-L1000中的扰动转录组学信号与CTRP中药物治疗后的细胞活性信息、Achilles中shRNA治疗前后的效应改变量进行匹配和筛选，得到在多点浓度和多个时间点下的差异表达基因的转录水平表达量与细胞表型信息的数据集；
S103：步骤S102得到的数据集分为3小时、6小时、24小时、96小时、120小时和144小时的数据子集。
3.根据权利要求1或2所述的细胞活性预测方法，其特征在于，还包括通过癌细胞系百科全书CCLE和NCI60数据集中的数据验证模型的有效性，具体为：
a)在NCI60数据集中，使用GI50值作为药物敏感性的评价指标，当药效在50％生长抑制浓度范围内，对应药物敏感性评价指标GI50值，当药效在50％生长抑制浓度范围内无效时，记为最高浓度值，定义药物浓度差异变量；
b)在CCLE数据集中，采用活性面积作为药物敏感性的评价指标，并进行二值化处理。
CN202010331009.6A 2020-04-24 2020-04-24 基于lincs-l1000扰动信号的细胞活性预测方法 Active CN111524554B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010331009.6A CN111524554B (zh) 2020-04-24 2020-04-24 基于lincs-l1000扰动信号的细胞活性预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010331009.6A CN111524554B (zh) 2020-04-24 2020-04-24 基于lincs-l1000扰动信号的细胞活性预测方法

Publications (2)

Publication Number Publication Date
CN111524554A CN111524554A (zh) 2020-08-11
CN111524554B true CN111524554B (zh) 2023-03-24

Family

ID=71903982

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010331009.6A Active CN111524554B (zh) 2020-04-24 2020-04-24 基于lincs-l1000扰动信号的细胞活性预测方法

Country Status (1)

Country Link
CN (1) CN111524554B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112599218A (zh) * 2020-12-16 2021-04-02 北京深度制耀科技有限公司 药物敏感预测模型的训练方法、预测方法及相关装置
CN112466403B (zh) * 2020-12-31 2022-06-14 广州基迪奥生物科技有限公司 一种细胞通讯分析方法及系统
CN112820417B (zh) * 2021-01-26 2022-12-23 四川大学 一种基于转录组学的前列腺癌药物组合预测的方法
CN114373550B (zh) * 2022-03-21 2022-06-21 普瑞基准科技（北京）有限公司 基于分子结构及基因表达的药物ic50深度学习模型预测方法
CN114743685B (zh) * 2022-04-01 2024-01-05 中国医学科学院北京协和医院 一种基于人工智能的子宫内膜癌风险筛查方法及系统

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109055476A (zh) * 2018-08-07 2018-12-21 河南大学 一种生物活性肽p11在治疗甲状腺癌的药物中效果分析方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20080118576A1 (en) * 2006-08-28 2008-05-22 Dan Theodorescu Prediction of an agent's or agents' activity across different cells and tissue types
JP6382459B1 (ja) * 2015-06-15 2018-08-29 ナントミクス，エルエルシー 細胞系ゲノミクスからの薬物応答の患者特異的予測のためのシステムおよび方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109055476A (zh) * 2018-08-07 2018-12-21 河南大学 一种生物活性肽p11在治疗甲状腺癌的药物中效果分析方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Optimal_Intervention_in_Markovian_Gene_Regulatory_Networks_With_Random-Length_Therapeutic_Response_to_Antitumor_Drug;Mohammadmahdi R. Yousef;《IEEE》;20130711;全文 *
可诱导转基因植物雄性不育的研究;陈明;《CNKI》;20041215;全文 *

Also Published As

Publication number Publication date
CN111524554A (zh) 2020-08-11

Similar Documents

Publication Publication Date Title
CN111524554B (zh) 基于lincs-l1000扰动信号的细胞活性预测方法
Gulhan et al. Detecting the mutational signature of homologous recombination deficiency in clinical samples
Peng et al. Predicting drug response based on multi-omics fusion and graph convolution
Keegan et al. Meta-analysis of Drosophila circadian microarray studies identifies a novel set of rhythmically expressed genes
EP2387758B1 (en) Evolutionary clustering algorithm
CN109072309A (zh) 癌症进化检测和诊断
CN111933211B (zh) 癌症精准化疗分型标志物筛选方法、化疗敏感性的分子分型方法和应用
CN109952611A (zh) 达沙替尼响应预测模型及其方法
CN113555070A (zh) 机器学习算法构建急性髓系白血病药敏相关基因分类器
Segura-Lepe et al. Predictive modelling using pathway scores: robustness and significance of pathway collections
CN110349633A (zh) 一种基于辐射响应生物学通路筛选辐射生物标志物及预测辐射剂量的方法
Lin et al. A model-based constrained deep learning clustering approach for spatially resolved single-cell data
KR20200105069A (ko) 빅테이터를 이용한 조건별 마이크로 rna 표적 조사 방법
Kim et al. Prediction of acquired taxane resistance using a personalized pathway-based machine learning method
Su et al. Integrated analysis of ovarian cancer patients from prospective transcription factor activity reveals subtypes of prognostic significance
Chiang et al. Inferring the transcriptional regulatory mechanism of signal‐dependent gene expression via an integrative computational approach
Wang et al. Identification of cancer trait genes and association analysis under pan-cancer
Sethi et al. Deciphering common temporal transcriptional response during powdery mildew disease in plants using meta-analysis
Donakonda et al. System analysis identifies distinct and common functional networks governed by transcription factor ASCL1, in glioma and small cell lung cancer
US20140288846A1 (en) System and method to identify dysregulated pathways and related interactions
Chin et al. Cytokine-expression patterns reveal coordinated immunological programs associated with persistent MRSA bacteremia
Mégret et al. Precision machine learning to understand micro-RNA regulation in neurodegenerative diseases
Hossain et al. Estimation of weighted log partial area under the ROC curve and its application to MicroRNA expression data
CN114203255B (zh) 一种基于机器学习的中药抗癌关键靶标预测方法
Zhang et al. P53 pathway activate detection based on machine learning: The modified XGBoost-based method of pan-cancer pathway activity detection in the cancer genome atlas

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant