CN112270952A - 一种识别癌症驱动通路的方法 - Google Patents

一种识别癌症驱动通路的方法 Download PDF

Info

Publication number
CN112270952A
CN112270952A CN202011185104.6A CN202011185104A CN112270952A CN 112270952 A CN112270952 A CN 112270952A CN 202011185104 A CN202011185104 A CN 202011185104A CN 112270952 A CN112270952 A CN 112270952A
Authority
CN
China
Prior art keywords
matrix
gene
chromosome
mutation
value
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.)
Granted
Application number
CN202011185104.6A
Other languages
English (en)
Other versions
CN112270952B (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.)
Guangxi Normal University
Original Assignee
Guangxi Normal 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 Guangxi Normal University filed Critical Guangxi Normal University
Priority to CN202011185104.6A priority Critical patent/CN112270952B/zh
Publication of CN112270952A publication Critical patent/CN112270952A/zh
Application granted granted Critical
Publication of CN112270952B publication Critical patent/CN112270952B/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
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/10—Ploidy or copy number detection
    • 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
    • G16B20/50—Mutagenesis

Landscapes

  • Bioinformatics & Cheminformatics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Biotechnology (AREA)
  • Biophysics (AREA)
  • Genetics & Genomics (AREA)
  • Molecular Biology (AREA)
  • Chemical & Material Sciences (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Analytical Chemistry (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本发明公开了一种识别癌症驱动通路的方法，包括如下步骤：1）构造加权的非二进制突变矩阵；2）设定识别模型；3）设定适应度函数；4）设定交叉算子；5）设定变异算子；6）设定合作策略；7）设定参数；8）构造初始种群；9）执行迭代操作。这种方法可以提供更多有用信息，扩展性强，速率高，求解速度快，包含富集在重要驱动通路上的基因多。

Description

一种识别癌症驱动通路的方法
技术领域
本发明涉及癌症驱动通路识别领域，具体是一种识别癌症驱动通路的方法。
背景技术
随着深度测序技术的快速发展，近年来，癌症基因组图谱计划(The CancerGenome Atlas，简称TCGA)、国际肿瘤基因组协作组(the International Cancer GenomeConsortium，简称ICGC)等大规模的癌症工程，提供了大量癌症多组学数据。在过去的数项研究中，设计有效的计算方法来识别癌症发生过程中的“驱动突变”成为热点。然而，大多数方法都无法确定基因突变的异质性，即来自同一癌症的不同样本之间，得到的突变基因也可能不同。研究人员发现，不同的突变基因靶向同一生物通路几率很高，并且发现癌症的发展实质上是由生物通路控制。于是，有必要将观点从基因水平转移到通路水平，这对于捕获癌症中的异质现象至关重要。因此识别驱动通路的问题逐渐成为热点，该问题主要分为：单驱动通路识别、合作驱动通路识别以及泛癌驱动通路识别。
目前基于先验生物通路知识和从头识别方法识别癌症驱动通路是两种主要研究方法。由于当前可利用的先验知识不完备，并且存在生物数据噪声。因此，在不依赖任何现有背景知识的识别方法是不可或缺的，从头识别方法应运而生，本文主要针对该方法进行研究。在癌症突变的大多数组合模式下，研究人员发现两个基本特性。其一，突变驱动通路中的基因应覆盖大量癌症患者样本，即“高覆盖性”。其二，同一驱动通路中不同基因在同一样本以低概率共同突变，即“高互斥性”。基于以上两种特性。在2012年，首先基于这两个属性提出了一个权重函数，并将其最大化定义为最大权子矩阵问题，并用马尔可夫链蒙特卡洛MCMC方法来解决这个问题。同年，又提出了二进制线性规划算法和遗传算法GA，与MCMC相比GA更加具有时间上的优势，GA算法也很容易应用于解决他们提出的整合了基因表达谱的整合模型。2013年，通过整合体细胞突变，拷贝数变异和基因表达，提出了一种基于网络的方法iMCMC。2016年，设计了一种多目标优化算法MOGA，高覆盖范围和高互斥性之间的权衡使MOGA表现出更可靠的性能。在2019年，重新构造了最大权重子矩阵问题模型，通过使用路径中基因的平均权重来调节覆盖范围和互斥性。然后设计了单亲遗传算法PGA-MWS对该模型进行求解。在上述方法中，大多数方法都试图借助其他组学数据减轻突变数据中噪声的负面影响，并基于驱动程序路径的两个特征生成识别模型。尽管如此，在覆盖率和排他性的计算中通常只考虑突变数据，而其他组学数据则用于计算基因的权重以表明它们是否重要，这是很常见的。
因此需要一种新的组学整合方法对各个数据进行融合，而不是在计算驱动通路权重时将其分开计算，设计合理的计算模型对驱动通路进行权重计算，并且以往的算法在数据规模较大的情况下，效率并不是很高，因此设计一种更高效的算法对计算模型进行求解，以解决现有方法的不足。
发明内容
本发明的目的在于克服现有技术的不足，提供一种识别癌症驱动通路的方法。这种方法可以提供更多有用信息，扩展性强，速率高，求解速度快，包含富集在重要驱动通路上的基因多。
实现本发明目的的技术方案是：
一种识别癌症驱动通路的方法，包括如下步骤：
1)构造加权的非二进制突变矩阵：
现有某个癌症的体细胞突变矩阵
Figure BDA0002751210320000021
拷贝数变异矩阵
Figure BDA0002751210320000022
和基因表达矩阵
Figure BDA0002751210320000023
在体细胞突变矩阵
Figure BDA0002751210320000024
拷贝数变异矩阵
Figure BDA0002751210320000025
和基因表达矩阵
Figure BDA0002751210320000026
三个矩阵中行表示一种癌症的相同样本集p，列分别表示基因集GS、GC和GE，在矩阵
Figure BDA0002751210320000027
中，sij∈{0,1}(i＝1,2,…,|p|,j＝1,2,…,|GS|)，i样本中j基因突变，sij值为1，反之值为0；矩阵
Figure BDA0002751210320000028
中每个元素cij∈{-2,-1,0,1,2}(i＝1,2,…,|p|,j＝1,2,…,|GC|)，表示i样本中j基因拷贝数变异值；在矩阵
Figure BDA0002751210320000029
中eij∈R(i＝1,2,…,|p|,j＝1,2,…,|GE|)，表示i样本中j基因表达量；令矩阵
Figure BDA00027512103200000210
中的基因集GA＝GS∪GC，样本集为p，令
Figure BDA00027512103200000211
aij∈{0,1}(i＝1,2,…,|p|,j＝1,2,…,|GA|)，其中
Figure BDA00027512103200000212
为突变矩阵，当sij取值为1或i样本中j基因处于统计显著变异区域时，aij值为1，反之值为0，为了进一步整合突变矩阵
Figure BDA00027512103200000213
和表达矩阵
Figure BDA00027512103200000214
在突变矩阵
Figure BDA00027512103200000215
和基因表达矩阵
Figure BDA00027512103200000216
中取基因集G＝GA∩GE，重新得到两个矩阵A|p|×|G|和E|p|×|G|，对于基因表达数据，存在正常样本表达矩阵N|n|×|G|，n表示正常样本，在矩阵N|n|×|G|中，nij∈R(i＝1,2,…,|p|,j＝1,2,…,|G|)，表示i样本中j基因表达量，令差异倍数矩阵D|p|×|G|，dij∈R(i＝1,2,…,|p|,j＝1,2,×,|G|)，表示i样本中j基因表达量相比j基因在正常样本中表达量的差异倍数用
Figure BDA00027512103200000217
表示，其中
Figure BDA00027512103200000218
则dij值为
Figure BDA0002751210320000031
否则dij值为0，处理好差异倍数矩阵D|p|×|G|，进一步对突变矩阵A|p|×|G|进行加权处理，整合成加权突变矩阵，对于A|p|×|G|，如果aij＝1，并且dij≥λ1，则aij＝1.5，如果aij＝0，并且dij≥λ2，则aij＝(2·l)-1·dij，其中λ1和λ2是截取差异倍数的阈值，l是j基因对应所有样本中差异倍数的最大值，针对突变基因，λ1取较低值，使aij∈{1,1.5}，以提高该突变基因的突变可信值；针对不突变基因，λ2取较高值，使aij∈[0,0.5]，以提高该不突变基因的突变可信值，使其可能成为潜在基因，经过加权重新得到加权突变矩阵A|p|×|G|，aij∈[0,1.5](i＝1,2,…,|p|,j＝1,2,…,|G|)；
2)设定识别模型：
针对加权突变矩阵A|p|×|G|，基于高覆盖和高互斥两个特性，重新构建新的整合模型，假设M|p|×k为矩阵A|p|×|G|的任一子矩阵，令Γ(m)＝{mi|mi＝max{aim|m∈M},i＝1,2,…,|p|}记录矩阵M|p|×k每行中最大权值，令矩阵M的覆盖度
Figure BDA0002751210320000032
对于矩阵M|p|×k中一行的互斥度，考虑这一行的离散程度，用变异系数计算每行的互斥度，每行互斥度之和为整个M|p|×k的互斥度，具体表示如公式(1)所示：
Figure BDA0002751210320000033
其中
Figure BDA0002751210320000034
当
Figure BDA0002751210320000035
趋近于0值时，对于变异系数值影响很大，如果矩阵M|p|×k中一行的最大权值mi≤0.5，则令该行的互斥度为
Figure BDA0002751210320000036
放缩该行互斥度，避免该行互斥度过高对通路识别造成影响，根据
Figure BDA0002751210320000037
和公式(1)，对整合数据后的最大权重子矩阵问题重新定义模型：给定突变加权矩阵A|p|×|G|和正整数k(k＜|G|)，在矩阵A|p|×|G|中确定矩阵M|p|×k，使函数值W(M)最大，如公式(2)所示：
W(M)＝α(M)+ω(M) (2)，
其中α(M)表示矩阵M|p|×k的覆盖度，α(M)由M|p|×k中每行最大突变权值相加得到，α(M)越大表示覆盖样本越多，并且突变可信值也越大；ω(M)表示矩阵M|p|×k的互斥度，ω(M)由矩阵M|p|×k中每行变异系数值相加得到，变异系数越大，则离散程度越高，其互斥度就越大；
3)设定适应度函数：
每个染色体对应一个问题解，因此需要对该解进行评估，给定染色体Xi(i＝1，2，…P)，P为种群大小，适应度函数Fitness(Xi)的定义如公式(3)所示：
Figure BDA0002751210320000041
其中，
Figure BDA0002751210320000042
表示染色体Xi对应的子矩阵；
4)设定交叉算子：
交叉算子，决定了GA算法的全局搜索能力，因此保证种群多样性，能有效提高搜索能力，本技术方案方法按照染色体适应度从大到小，给染色体Xi一个排名Ri，则每个染色体被选中的概率如公式(4)所示：
Figure BDA0002751210320000043
为了保证染色体的可行性，采用轮盘赌随机从父种群选取两条染色体，重复基因分别给子代的两条染色体，剩余基因放在一个集合，采用均匀交叉的方式，对于剩余基因集合中连续的每对基因，随机生成一个二值数据，如果该二值数据为1，则这对基因的第一个基因放入第一条子染色体，第二个基因放入第二条子染色体，反之，第一个基因放入第二条子染色体，第二个基因放入第一条子染色体，经过一次交叉，生成两个子染色体；
5)设定变异算子：
给定一条子染色体X＝{x1，x2，…，xk}(xi＝1，2，…，|G|)，确定候选基因集合
Figure BDA0002751210320000045
从子染色体中随机删除一个基因，得到基因集X′，将HX中基因顺序打乱，遍历前
Figure BDA0002751210320000044
个基因，选出基因g，使适应值Fitness(MX′∪{g})最大，对应于子矩阵MX′∪{g}的基因集X′∪{g}为新子染色体，即X＝X′∪{g}；
6)设定合作策略：
采用种群间相互合作策略，在种群交叉、变异和选择操作后，比较两个种群适应度最好的染色体和对方适应度最差的染色体，如果最好适应度高于对方最差适应度，则将该条染色体替换对方适应度最差的染色体；
7)设定参数：
输入加权的非二进制突变矩阵A|p|×|G|，公式(2)中的模型，参数k用于限制找到的驱动通路大小，然后输入CGA-MWS相关参数：种群规模P、变异概率Pm、最大演化代数maxstep、最优值保持恒定的阈值maxt；
8)构造初始种群：
染色体编码采用十进制编码方式，一个解为k个基因构成的集合，即X＝{x1，x2，…，xk}(xi＝1，2，…，|G|)，将|G|个基因顺序随机打乱，然后取前k个基因构成初始染色体，生成两个初始种群pop0和pop1，每个种群大小为P/2，计算两个种群各自染色体的适应值，将pop0和pop1中各自最优的染色体相比较，保存最好的个体到变量best中，初始迭代次数step＝0，最优值保持恒定的代数t＝0；
9)执行迭代操作：
(1)若step＞max step或t＞maxt，转入步骤9)的(4)，得到大小为k的驱动通路，否则转入步骤9)的(2)；
(2)种群pop0和pop1各自通过基于排名的概率，采用轮盘赌随机选两条父代染色体，通过交叉算子进行交叉，生成两条子染色体，并各自放入子种群pop0′和pop1′中，重复P/4次，针对子种群pop0′和pop1′，对每条染色体随机一个突变概率Pm′，如果Pm′＜Pm，则对该条染色体进行变异操作，将变异操作中得到的适应值最高的染色体替换该条染色体，将pop0和pop0′里所有染色体按适应值从高到低排序，取前P/2条染色体放入下一代种群popstep+1中，pop1和pop1′进行相同操作得到下一代种群popstep+2；
(3)对种群popstep+1和popstep+2进行合作策略操作，比较两个种群的最优适应值，取两者中适应值最高的染色体，若该染色体适应值大于best染色体的适应值，则更新best染色体，t＝0；否则t＝t+1，step＝step+1，返回步骤9)的(1)；
(4)将best染色体转换为基因集，由此得到子矩阵M，并子矩阵M其输出，输出的子矩阵M即为大小为k的驱动通路。
本技术方案具有以下优点：
(1)新的多组学数据融合方法，整合了组学数据之间的特点，可提供更多有用信息。
(2)新的模型可针对非二元矩阵进行模型求解，更具有扩展性。
(3)新的识别算法，对大规模数据，有更高效的速率，可更快对模型进行求解。
(4)整个方法找出的单癌种驱动通路，包含更多富集在同一个通路的基因。
这种方法可以提供更多有用信息，扩展性强，速率高，求解速度快，富集在重要驱动通路上的基因多。
附图说明
图1为实施例1中胶质母细胞瘤GBM，驱动通路规模为3时运行图；
图2为实施例2中胶质母细胞瘤GBM，驱动通路规模为10时运行图；
具体实施方式
下面结合附图和实施例对本发明做进一步阐述，但不是对本发明的限定。
本例针对单驱动通路识别问题进行阐述。
实施例1：
一种识别癌症驱动通路的方法，包括如下步骤：
1)构造加权的非二进制突变矩阵：
现有胶质母细胞瘤GBM体细胞突变矩阵
Figure BDA0002751210320000061
拷贝数变异矩阵
Figure BDA0002751210320000062
和基因表达矩阵
Figure BDA0002751210320000063
中突变矩阵
Figure BDA0002751210320000064
拷贝数变异矩阵
Figure BDA0002751210320000065
和基因表达矩阵
Figure BDA0002751210320000066
中行表示一种癌症的相同样本集p，列分别表示基因集GS、GC和GE，在矩阵
Figure BDA0002751210320000067
中，sij∈{0，1}(i＝1，2，…，|p|，j＝1，2，…，|GS|)，i样本中j基因突变，sij值为1，反之值为0；矩阵
Figure BDA0002751210320000068
中每个元素cij∈{-2，-1，0，1，2}(i＝1，2，…，|p|，j＝1，2，…，|GC|)，表示i样本中j基因拷贝数变异值；在矩阵
Figure BDA0002751210320000069
中eij∈R(i＝1，2，…，|p|，j＝1，2，…，|GE|)，表示i样本中j基因表达量；令矩阵
Figure BDA00027512103200000610
中的基因集GA＝GS∪GC，样本集为p，令
Figure BDA00027512103200000611
aij∈{0，1}(i＝1，2，…，|p|，j＝1，2，…，|GA|)，其中
Figure BDA0002751210320000071
为突变矩阵，当sij取值为1或i样本中j基因处于统计显著变异区域时，aij值为1，反之值为0，为了进一步整合突变矩阵
Figure BDA0002751210320000072
和表达矩阵
Figure BDA0002751210320000073
中突变矩阵
Figure BDA0002751210320000074
和表达矩阵
Figure BDA0002751210320000075
中取基因集G＝GA∩GE，重新得到两个矩阵A|p|×|G|和E|p|×|G|，对于基因表达数据，存在正常样本表达矩阵N|n|×|G|，n表示正常样本，在矩阵N|n|×|G|中，nij∈R(i＝1，2，…，|p|，j＝1，2，…，|G|)，表示i样本中j基因表达量，令差异倍数矩阵D|p|×|G|，dij∈R(i＝1，2，…，|p|，j＝1，2，…，|G|)，表示i样本中j基因表达量相比j基因在正常样本中表达量的差异倍数用
Figure BDA0002751210320000076
表示，其中
Figure BDA0002751210320000077
则dij值为
Figure BDA0002751210320000078
否则dij值为0，处理好差异倍数矩阵D|p|×|G|，进一步对突变矩阵A|p|×|G|进行加权处理，整合成加权突变矩阵，对于A|p|×|G|，如果αij＝1，并且dij≥λ1，此时λ1＝3，则aij＝1.5，如果aij＝0，并且dij≥λ2，此时λ2＝7，则aij＝(2·l)-1·dij，其中λ1和λ2是截取差异倍数的阈值，l是j基因对应所有样本中差异倍数的最大值，针对突变基因，λ1取较低值，使aij∈{1，1.5}，以提高该突变基因的突变可信值；针对不突变基因，λ2取较高值，使aij∈[0，0.5]，以提高该不突变基因的突变可信值，使其可能成为潜在基因，经过加权重新得到加权突变矩阵A|p|×|G|，aij∈[0，1.5](i＝1，2，…，|p|，j＝1，2，…，|G|)，其中|p|＝90，|G|＝920；
2)设定识别模型：
针对加权突变矩阵A|p|×|G|，基于高覆盖和高互斥两个特性，重新构建新的整合模型，假设M|p|×k为矩阵A|p|×|G|的任一子矩阵，令Γ(m)＝{mi|mi＝max{aim|m∈M}，i＝1，2，…，|p|}记录矩阵M|p|×k每行中最大权值，令矩阵M的覆盖度
Figure BDA0002751210320000079
对于矩阵M|p|×k中一行的互斥度，考虑这一行的离散程度，用变异系数计算每行的互斥度，每行互斥度之和为整个M|p|×k的互斥度，具体表示如公式(1)所示：
Figure BDA00027512103200000710
其中
Figure BDA00027512103200000711
当
Figure BDA00027512103200000712
趋近于0值时，对于变异系数值影响很大，如果矩阵M|p|×k中一行的最大权值mi≤0.5，则令该行的互斥度为
Figure BDA0002751210320000081
放缩该行互斥度，避免该行互斥度过高对通路识别造成影响，根据的
Figure BDA0002751210320000082
和公式(1)，对整合数据后的最大权重子矩阵问题重新定义模型：给定突变加权矩阵A|p|×|G|和正整数k(k＜|G|)，在矩阵A|p|×|G|中确定矩阵M|p|×k，使函数值W(M)最大，如公式(2)所示：
W(M)＝α(M)+ω(M) (2)，
其中α(M)表示矩阵M|p|×k的覆盖度，α(M)由M|p|×k中每行最大突变权值相加所得，α(M)越大表示覆盖样本越多，并且突变可信值也越大；ω(M)表示矩阵M|p|×k的互斥度，ω(M)由矩阵M|p|×k中每行变异系数值相加所得，变异系数越大，则离散程度越高，其互斥度就越大；
3)设定适应度函数：
每个染色体对应一个问题解，因此需要对该解进行评估，给定染色体Xi(i＝1，2，…P)，P为种群大小，适应度函数Fitness(Xi)的定义如公式(3)所示：
Figure BDA0002751210320000083
其中，
Figure BDA0002751210320000084
表示染色体Xi对应的子矩阵；
4)设定交叉算子：
交叉算子，决定了GA算法的全局搜索能力，因此保证种群多样性，能有效提高搜索能力，本例方法按照染色体适应度从大到小，给染色体Xi一个排名Ri，则每个染色体被选中的概率如公式(4)所示：
Figure BDA0002751210320000085
为了保证染色体的可行性，采用轮盘赌随机从父种群选取两条染色体，重复基因分别给子代的两条染色体，剩余基因放在一个集合，采用均匀交叉的方式，对于剩余基因集合中连续的每对基因，随机生成一个二值数据，如果该二值数据值为1，则这对基因的第一个基因放入第一条子染色体，第二个基因放入第二条子染色体，反之，第一个基因放入第二条子染色体，第二个基因放入第一条子染色体，经过一次交叉，生成两个子染色体；
5)设定变异算子：
给定一条子染色体X＝{x1，x2，…，xk}(xi＝1，2，…，|G|)，确定候选基因集合
Figure BDA0002751210320000092
从子染色体中随机删除一个基因，得到基因集X′，将HX中基因顺序打乱，遍历前
Figure BDA0002751210320000091
个基因，选出基因g，使适应值Fitness(MX′∪{g})最大，对应于子矩阵Mx′∪{g}的基因集X′∪{g}为新子染色体，即X＝X′∪{g}；
6)设定合作策略：
采用种群间相互合作策略，在种群交叉、变异和选择操作后，比较两个种群适应度最好的染色体和对方适应度最差的染色体，如果最好适应度高于对方最差适应度，则将该条染色体替换对方适应度最差的染色体；
7)设定参数：
输入加权的非二进制突变矩阵A|p|×|G|，其中|p|＝90，|G|＝920，公式(2)中的模型，参数k＝3用于限制找到的驱动通路大小，然后输入CGA-MWS相关参数：种群规模P＝460、变异概率Pm＝0.3、最大演化代数maxstep＝1000、最优值保持恒定的阈值maxt＝10；
8)构造初始种群：
染色体编码采用十进制编码方式，一个解为k＝3个基因构成的集合，即X＝{x1，x2，…，xk}(xi＝1，2，…，|G|)，将|G|＝920个基因顺序随机打乱，然后取前k＝3个基因构成初始染色体，生成两个初始种群pop0和pop1，每个种群大小为230，计算两个种群各自染色体的适应值，将pop0和pop1中各自最优的染色体相比较，保存最好的个体到变量best中，初始迭代次数step＝0，最优值保持恒定的代数t＝0；
9)执行迭代操作：
(1)若step＞max step或t＞maxt，转入步骤9)的(4)，得到大小为k＝3的驱动通路，否则转入步骤9)的(2)；
(2)种群pop0和pop1各自通过基于排名的概率，采用轮盘赌随机选两条父代染色体，通过交叉算子进行交叉，生成两条子染色体，并各自放入子种群pop0′和pop1′中，重复P/4次，针对子种群pop0′和pop1′，对每条染色体随机一个突变概率Pm′，如果Pm′＜Pm，则对该条染色体进行变异操作，将变异操作中得到的适应值最高的染色体替换该条染色体，将pop0和pop0′里所有染色体按适应值从高到低排序，取前P/2条染色体放入下一代种群popstep+1中，pop1和pop1′进行相同操作得到下一代种群popstep+2；
(3)对种群popstep+1和popstep+2进行合作策略操作，比较两个种群的最优适应值，取两者中适应值最高的染色体，若该染色体适应值大于best染色体的适应值，则更新best染色体，t＝0；否则t＝t+1，step＝step+1，返回步骤9)的(1)；
(4)将best染色体转换为基因集，由此得到子矩阵M，并子矩阵M其输出，输出的子矩阵M即为大小为k＝3的驱动通路，运行图如图1所示。
实施例2：
本例步骤1)中设定λ1＝3和λ2＝7，构造加权非二进制突变矩阵A|p|×|G|，其中|p|＝90，|G|＝920；
本例步骤7)输入加权的非二进制突变矩阵A|p|×|G|，其中|p|＝90，|G|＝920，公式(2)中的模型，需要寻找的驱动通路大小k＝10，CGA-MWS相关参数，设置种群规模P＝460、变异概率Pm＝0.3、最大演化代数maxstep＝1000、最优值保持恒定的阈值maxt＝10；
本例步骤9)得到大小为k＝10的驱动通路，运行图如图2所示。
其余步骤同实施例1。

Claims (1)

1.一种识别癌症驱动通路的方法，其特征在于，包括如下步骤：
1)构造加权的非二进制突变矩阵：
现有某个癌症的体细胞突变矩阵
Figure FDA0002751210310000011
拷贝数变异矩阵
Figure FDA0002751210310000012
和基因表达矩阵
Figure FDA0002751210310000013
在体细胞突变矩阵
Figure FDA0002751210310000014
拷贝数变异矩阵
Figure FDA0002751210310000015
和基因表达矩阵
Figure FDA0002751210310000016
三个矩阵中行表示该癌症的相同样本集p，列分别表示基因集GS、GC和GE，在矩阵
Figure FDA0002751210310000017
中，sij∈{0,1}(i＝1,2,…,|p|,j＝1,2,…,|GS|)，i样本中j基因突变，sij值为1，反之值为0；矩阵
Figure FDA0002751210310000018
中每个元素cij∈{-2,-1,0,1,2}(i＝1,2,…,|p|,j＝1,2,…,|GC|)，表示i样本中j基因拷贝数变异值；在矩阵
Figure FDA0002751210310000019
中eij∈R(i＝1,2,…,|p|,j＝1,2,…,|GE|)，表示i样本中j基因表达量；令矩阵
Figure FDA00027512103100000110
中的基因集为GA＝GS∪GC，样本集为p，令
Figure FDA00027512103100000111
aij∈{0,1}(i＝1,2,…,|p|,j＝1,2,…,|GA|)，其中
Figure FDA00027512103100000112
为突变矩阵，当sij取值为1或i样本中j基因处于统计显著变异区域时，aij值为1，反之值为0，为了进一步整合突变矩阵
Figure FDA00027512103100000113
和表达矩阵
Figure FDA00027512103100000114
在突变矩阵
Figure FDA00027512103100000115
和表达矩阵
Figure FDA00027512103100000116
中取基因集G＝GA∩GE，重新得到两个矩阵A|p|×|G|和E|p|×|G|，对于基因表达数据，存在正常样本表达矩阵N|n|×|G|，n表示正常样本，在矩阵N|n|×|G|中，nij∈R(i＝1,2,…,|p|,j＝1,2,…,|G|)，表示i样本中j基因表达量，令差异倍数矩阵D|p|×|G|，dij∈R(i＝1,2,…,|p|,j＝1,2,…,|G|)，表示i样本中j基因表达量相比j基因在正常样本中表达量的差异倍数用
Figure FDA00027512103100000117
表示，其中
Figure FDA00027512103100000118
则dij值为
Figure FDA00027512103100000119
否则dij值为0，处理好差异倍数矩阵D|p|×|G|，进一步对突变矩阵A|p|×|G|进行加权处理，整合成加权突变矩阵，对于A|p|×|G|，如果aij＝1，并且dij≥λ1，则aij＝1.5，如果aij＝0，并且dij≥λ2，则aij＝(2·l)-1·dij，其中λ1和λ2是截取差异倍数的阈值，l是j基因对应所有样本中差异倍数的最大值，针对突变基因，λ1取较低值，使aij∈{1,1.5}，以提高该突变基因的突变可信值；针对不突变基因，λ2取较高值，使aij∈[0,0.5]，以提高该不突变基因的突变可信值，使其可能成为潜在基因，经过加权重新得到加权突变矩阵A|p|×|G|，aij∈[0,1.5](i＝1,2,…,|p|,j＝1,2,…,|G|)；
2)设定识别模型：
针对加权突变矩阵A|p|×|G|，基于高覆盖和高互斥两个特性，重新构建新的整合模型，假设M|p|×k为矩阵A|p|×|G|的任一子矩阵，令Γ(m)＝{mi|mi＝max{aim|m∈M},i＝1,2,…,|p|}记录矩阵M|p|×k每行中最大权值，令矩阵M|p|×k的覆盖度
Figure FDA0002751210310000021
对于矩阵M|p|×k中一行的互斥度，考虑这一行的离散程度，用变异系数计算每行的互斥度，每行互斥度之和为整个M|p|×k的互斥度，具体表示如公式(1)所示：
Figure FDA0002751210310000022
其中
Figure FDA0002751210310000023
当
Figure FDA0002751210310000024
趋近于0值时，对于变异系数值影响很大，所以如果M|p|×k中一行的最大权值mi≤0.5，则令该行的互斥度为
Figure FDA0002751210310000025
放缩该行互斥度，避免该行互斥度过高对通路识别造成影响，根据
Figure FDA0002751210310000026
和公式(1)，对整合数据后的最大权重子矩阵问题重新定义模型：给定突变加权矩阵A|p|×|G|和正整数k(k＜|G|)，在矩阵A|p|×|G|中确定矩阵M|p|×k，使函数值W(M)最大，如公式(2)所示：
W(M)＝α(M)+ω(M) (2)，
其中α(M)表示矩阵M|p|×k的覆盖度，α(M)由矩阵M|p|×k中每行最大突变权值相加所得，α(M)越大表示覆盖样本越多，并且突变可信值也越大；ω(M)表示矩阵M|p|×k的互斥度，ω(M)由M中每行变异系数值相加所得，变异系数越大，则离散程度越高，其互斥度就越大；
3)设定适应度函数：
每个染色体对应一个问题解，因此需要对该解进行评估，给定染色体Xi(i＝1,2,…P)，P为种群大小，适应度函数Fitness(Xi)的定义如公式(3)所示：
Figure FDA0002751210310000027
其中，
Figure FDA0002751210310000028
表示染色体Xi对应的子矩阵；
4)设定交叉算子：
按照染色体适应度从大到小，给染色体Xi一个排名Ri，则每个染色体被选中的概率如公式(4)所示：
Figure FDA0002751210310000031
为了保证染色体的可行性，采用轮盘赌随机从父种群选取两条染色体，重复基因分别给子代的两条染色体，剩余基因放在一个集合，采用均匀交叉的方式，对于剩余基因的集合中连续的每对基因，随机生成一个二值数据，如果该二值数据为1，则这对基因的第一个基因放入第一条子染色体，第二个基因放入第二条子染色体，反之，第一个基因放入第二条子染色体，第二个基因放入第一条子染色体，经过一次交叉，生成两个子染色体；
5)设定变异算子：
给定一条子染色体X＝{x1,x2,…,xk}(xi＝1,2,…,|G|)，确定候选基因集合
Figure FDA0002751210310000032
从子染色体中随机删除一个基因，得到基因集X′，将HX中基因顺序打乱，遍历前
Figure FDA0002751210310000033
个基因，选出基因g，使适应值Fitness(MX′∪{g})最大，对应于子矩阵MX′∪{g}的基因集X′∪{g}为新子染色体，即X＝X′∪{g}；
6)设定合作策略：
采用种群间相互合作策略，在种群交叉、变异和选择操作后，比较两个种群适应度最好的染色体和对方适应度最差的染色体，如果最好适应度高于对方最差适应度，则将该条染色体替换对方适应度最差的染色体；
7)设定参数：
输入加权的非二进制突变矩阵A|p|×|G|和公式(2)中的模型，参数k用于限制找到的驱动通路大小，然后输入CGA-MWS相关参数：种群规模P、变异概率Pm、最大演化代数maxstep、最优值保持恒定的阈值maxt；
8)构造初始种群：
染色体编码采用十进制编码方式，一个解为k个基因构成的集合，即X＝{x1,x2,…,xk}(xi＝1,2,…,|G|)，将|G|个基因顺序随机打乱，然后取前k个基因构成初始染色体，生成两个初始种群pop0和pop1，每个种群大小为P/2，计算两个种群各自染色体的适应值，将pop0和pop1中各自最优的染色体相比较，保存最好的个体到变量best中，初始迭代次数step＝0，最优值保持恒定的代数t＝0；
9)执行迭代操作：
(1)若step＞maxstep或t＞maxt，转入步骤9)的(4)，得到大小为k的驱动通路，否则转入步骤9)的(2)；
(2)种群pop0和pop1各自通过基于排名的概率，采用轮盘赌随机选两条父代染色体，通过交叉算子进行交叉，生成两条子染色体，并各自放入子种群pop0′和pop1′中，重复P/4次，针对子种群pop0′和pop1′，对每条染色体随机一个突变概率Pm′，如果Pm′＜Pm，则对该条染色体进行变异操作，将变异操作中得到的适应值最高的染色体替换该条染色体，将pop0和pop0′里所有染色体按适应值从高到低排序，取前P/2条染色体放入下一代种群popstep+1中，pop1和pop1′进行相同操作得到下一代种群popstep+2；
(3)对种群popstep+1和popstep+2进行合作策略操作，比较两个种群的最优适应值，取两者中适应值最高的染色体，若该染色体适应值大于best染色体的适应值，则更新best染色体，t＝0；否则t＝t+1，step＝step+1，返回步骤9)的(1)；
(4)将best染色体转换为基因集，由此得到子矩阵M，并子矩阵M其输出，输出的子矩阵M即为大小为k的驱动通路。
CN202011185104.6A 2020-10-30 2020-10-30 一种识别癌症驱动通路的方法 Active CN112270952B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011185104.6A CN112270952B (zh) 2020-10-30 2020-10-30 一种识别癌症驱动通路的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011185104.6A CN112270952B (zh) 2020-10-30 2020-10-30 一种识别癌症驱动通路的方法

Publications (2)

Publication Number Publication Date
CN112270952A true CN112270952A (zh) 2021-01-26
CN112270952B CN112270952B (zh) 2022-04-05

Family

ID=74345286

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011185104.6A Active CN112270952B (zh) 2020-10-30 2020-10-30 一种识别癌症驱动通路的方法

Country Status (1)

Country Link
CN (1) CN112270952B (zh)

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2016095093A1 (zh) * 2014-12-15 2016-06-23 天津华大基因科技有限公司 肿瘤筛查方法、目标区域变异检测方法和装置
CN107003300A (zh) * 2014-12-12 2017-08-01 凯尔科迪公司 测量erbb信号传导通路活性以诊断和治疗癌症患者的方法
WO2017181134A2 (en) * 2016-04-15 2017-10-19 F. Hoffman-La Roche Ag Detecting cancer driver genes and pathways
CN108292326A (zh) * 2015-08-27 2018-07-17 皇家飞利浦有限公司 用于使用多组学癌症谱来识别功能性患者特异性体细胞畸变的整合方法和系统
US20190040473A1 (en) * 2013-04-18 2019-02-07 Gencurix Inc. Genetic marker for early breast cancer prognosis prediction and diagnosis, and use thereof
CN112259163A (zh) * 2020-10-28 2021-01-22 广西师范大学 基于生物网络和亚细胞定位数据识别癌症驱动模块方法
CN113228194A (zh) * 2018-10-12 2021-08-06 人类长寿公司 用于癌症基因组和临床数据综合分析的多组学搜索引擎
CN114023383A (zh) * 2021-11-04 2022-02-08 广西师范大学 一种识别癌症驱动通路的无参非线性智能优化方法

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20190040473A1 (en) * 2013-04-18 2019-02-07 Gencurix Inc. Genetic marker for early breast cancer prognosis prediction and diagnosis, and use thereof
CN107003300A (zh) * 2014-12-12 2017-08-01 凯尔科迪公司 测量erbb信号传导通路活性以诊断和治疗癌症患者的方法
WO2016095093A1 (zh) * 2014-12-15 2016-06-23 天津华大基因科技有限公司 肿瘤筛查方法、目标区域变异检测方法和装置
CN108292326A (zh) * 2015-08-27 2018-07-17 皇家飞利浦有限公司 用于使用多组学癌症谱来识别功能性患者特异性体细胞畸变的整合方法和系统
WO2017181134A2 (en) * 2016-04-15 2017-10-19 F. Hoffman-La Roche Ag Detecting cancer driver genes and pathways
CN113228194A (zh) * 2018-10-12 2021-08-06 人类长寿公司 用于癌症基因组和临床数据综合分析的多组学搜索引擎
CN112259163A (zh) * 2020-10-28 2021-01-22 广西师范大学 基于生物网络和亚细胞定位数据识别癌症驱动模块方法
CN114023383A (zh) * 2021-11-04 2022-02-08 广西师范大学 一种识别癌症驱动通路的无参非线性智能优化方法

Non-Patent Citations (7)

* Cited by examiner, † Cited by third party
Title
ZIYING YANG等: "CoDP: Cooperative Driver Pathways Discovery With Matrix Factorization and Tri-Random Walk", 《IEEE ACCESS》 *
吴昊: "基于突变基因网络的致癌驱动通路检测算法", 《计算机学报》 *
喵偷鹰: "通过融合基因、miRNA和代谢通路的多重关系数据联合分析挖掘的肿瘤驱动代谢途径", 《HTTPS://ZHUANLAN.ZHIHU.COM/P/112537964》 *
朱凯: "基于加权突变矩阵的癌症驱动通路识别模型与算法研究", 《中国优秀硕士学位论文全文数据库 医药卫生科技辑》 *
胡灿军: "癌症驱动通路的启发式识别方法研究", 《中国优秀硕士学位论文全文数据库 医药卫生科技辑》 *
蔡齐荣等: "基于多组学数据识别癌症驱动通路的模型和算法", 《计算机科学》 *
郭炳等: "一种基于突变基因网络的癌症驱动通路识别算法", 《计算机科学》 *

Also Published As

Publication number Publication date
CN112270952B (zh) 2022-04-05

Similar Documents

Publication Publication Date Title
Caye et al. TESS3: fast inference of spatial population structure and genome scans for selection
Otwinowski et al. Genotype to phenotype mapping and the fitness landscape of the E. coli lac promoter
US20220301658A1 (en) Machine learning driven gene discovery and gene editing in plants
EP3611799A1 (en) Array element arrangement method for l-type array antenna based on inheritance of acquired characteristics
CN109448794B (zh) 一种基于遗传禁忌和贝叶斯网络的上位性位点挖掘方法
Yadav et al. An overview of genetic algorithm and modeling
CN114420211A (zh) 一种基于注意力机制的rna-蛋白质结合位点预测方法
CN112734051A (zh) 一种针对分类问题的进化集成学习方法
CN104573004B (zh) 一种基于双阶遗传计算的基因表达数据的双聚类方法
Oluoch et al. A review on RNA secondary structure prediction algorithms
CN108509764B (zh) 一种基于遗传属性约简的古生物谱系演化分析方法
Banka et al. Evolutionary biclustering of gene expressions
CN112270952B (zh) 一种识别癌症驱动通路的方法
CN114093420B (zh) 一种基于XGBoost的DNA重组位点预测方法
Sheng et al. A survey of computational methods and databases for lncRNA-miRNA interaction prediction
CN114093426B (zh) 基于基因调控网络构建的标志物筛选方法
Peliti Fitness landscapes and evolution
KR20230043071A (ko) 변이체 병원성 채점 및 분류 그리고 이의 사용
CN114023383A (zh) 一种识别癌症驱动通路的无参非线性智能优化方法
Chipman et al. Using stochastic causal trees to augment Bayesian networks for modeling eQTL datasets
Huang Sparse model learning for inferring genotype and phenotype associations
Wu et al. Multiple sequence alignment using ga and nn
CN114065933B (zh) 一种基于人工免疫思想的未知威胁检测方法
Wu et al. Haplotyping a single triploid individual based on genetic algorithm
Dheenathayalan et al. Identifying significant genes from DNA microarray using genetic algorithm

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