CN112270952B - 一种识别癌症驱动通路的方法 - Google Patents
一种识别癌症驱动通路的方法 Download PDFInfo
- Publication number
- CN112270952B CN112270952B CN202011185104.6A CN202011185104A CN112270952B CN 112270952 B CN112270952 B CN 112270952B CN 202011185104 A CN202011185104 A CN 202011185104A CN 112270952 B CN112270952 B CN 112270952B
- Authority
- CN
- China
- Prior art keywords
- matrix
- gene
- mutation
- value
- chromosome
- 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
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)构造加权的非二进制突变矩阵:
现有某个癌症的体细胞突变矩阵拷贝数变异矩阵和基因表达矩阵在体细胞突变矩阵拷贝数变异矩阵和基因表达矩阵三个矩阵中行表示一种癌症的相同样本集p,列分别表示基因集GS、GC和GE,在矩阵中,sij∈{0,1}(i=1,2,…,|p|,j=1,2,…,|GS|),i样本中j基因突变,sij值为1,反之值为0;矩阵中每个元素cij∈{-2,-1,0,1,2}(i=1,2,…,|p|,j=1,2,…,|GC|),表示i样本中j基因拷贝数变异值;在矩阵中eij∈R(i=1,2,…,|p|,j=1,2,…,|GE|),表示i样本中j基因表达量;令矩阵中的基因集GA=GS∪GC,样本集为p,令aij∈{0,1}(i=1,2,…,|p|,j=1,2,…,|GA|),其中为突变矩阵,当sij取值为1或i样本中j基因处于统计显著变异区域时,aij值为1,反之值为0,为了进一步整合突变矩阵和表达矩阵在突变矩阵和基因表达矩阵中取基因集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基因在正常样本中表达量的差异倍数用表示,其中则dij值为否则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的覆盖度对于矩阵M|p|×k中一行的互斥度,考虑这一行的离散程度,用变异系数计算每行的互斥度,每行互斥度之和为整个M|p|×k的互斥度,具体表示如公式(1)所示:
其中当趋近于0值时,对于变异系数值影响很大,如果矩阵M|p|×k中一行的最大权值mi≤0.5,则令该行的互斥度为放缩该行互斥度,避免该行互斥度过高对通路识别造成影响,根据和公式(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)所示:
4)设定交叉算子:
交叉算子,决定了GA算法的全局搜索能力,因此保证种群多样性,能有效提高搜索能力,本技术方案方法按照染色体适应度从大到小,给染色体Xi一个排名Ri,则每个染色体被选中的概率如公式(4)所示:
为了保证染色体的可行性,采用轮盘赌随机从父种群选取两条染色体,重复基因分别给子代的两条染色体,剩余基因放在一个集合,采用均匀交叉的方式,对于剩余基因集合中连续的每对基因,随机生成一个二值数据,如果该二值数据为1,则这对基因的第一个基因放入第一条子染色体,第二个基因放入第二条子染色体,反之,第一个基因放入第二条子染色体,第二个基因放入第一条子染色体,经过一次交叉,生成两个子染色体;
5)设定变异算子:
给定一条子染色体X={x1,x2,…,xk}(xi=1,2,…,|G|),确定候选基因集合从子染色体中随机删除一个基因,得到基因集X′,将HX中基因顺序打乱,遍历前个基因,选出基因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体细胞突变矩阵拷贝数变异矩阵和基因表达矩阵中突变矩阵拷贝数变异矩阵和基因表达矩阵中行表示一种癌症的相同样本集p,列分别表示基因集GS、GC和GE,在矩阵中,sij∈{0,1}(i=1,2,…,|p|,j=1,2,…,|GS|),i样本中j基因突变,sij值为1,反之值为0;矩阵中每个元素cij∈{-2,-1,0,1,2}(i=1,2,…,|p|,j=1,2,…,|GC|),表示i样本中j基因拷贝数变异值;在矩阵中eij∈R(i=1,2,…,|p|,j=1,2,…,|GE|),表示i样本中j基因表达量;令矩阵中的基因集GA=GS∪GC,样本集为p,令aij∈{0,1}(i=1,2,…,|p|,j=1,2,…,|GA|),其中为突变矩阵,当sij取值为1或i样本中j基因处于统计显著变异区域时,aij值为1,反之值为0,为了进一步整合突变矩阵和表达矩阵中突变矩阵和表达矩阵中取基因集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基因在正常样本中表达量的差异倍数用表示,其中则dij值为否则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的覆盖度对于矩阵M|p|×k中一行的互斥度,考虑这一行的离散程度,用变异系数计算每行的互斥度,每行互斥度之和为整个M|p|×k的互斥度,具体表示如公式(1)所示:
其中当趋近于0值时,对于变异系数值影响很大,如果矩阵M|p|×k中一行的最大权值mi≤0.5,则令该行的互斥度为放缩该行互斥度,避免该行互斥度过高对通路识别造成影响,根据的和公式(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)所示:
4)设定交叉算子:
交叉算子,决定了GA算法的全局搜索能力,因此保证种群多样性,能有效提高搜索能力,本例方法按照染色体适应度从大到小,给染色体Xi一个排名Ri,则每个染色体被选中的概率如公式(4)所示:
为了保证染色体的可行性,采用轮盘赌随机从父种群选取两条染色体,重复基因分别给子代的两条染色体,剩余基因放在一个集合,采用均匀交叉的方式,对于剩余基因集合中连续的每对基因,随机生成一个二值数据,如果该二值数据值为1,则这对基因的第一个基因放入第一条子染色体,第二个基因放入第二条子染色体,反之,第一个基因放入第二条子染色体,第二个基因放入第一条子染色体,经过一次交叉,生成两个子染色体;
5)设定变异算子:
给定一条子染色体X={x1,x2,…,xk}(xi=1,2,…,|G|),确定候选基因集合从子染色体中随机删除一个基因,得到基因集X′,将HX中基因顺序打乱,遍历前个基因,选出基因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)构造加权的非二进制突变矩阵:
现有某个癌症的体细胞突变矩阵拷贝数变异矩阵和基因表达矩阵在体细胞突变矩阵拷贝数变异矩阵和基因表达矩阵三个矩阵中行表示该癌症的相同样本集p,列分别表示基因集GS、GC和GE,在矩阵中,sij∈{0,1}(i=1,2,…,|p|,j=1,2,…,|GS|),i样本中j基因突变,sij值为1,反之值为0;矩阵中每个元素cij∈{-2,-1,0,1,2}(i=1,2,…,|p|,j=1,2,…,|GC|),表示i样本中j基因拷贝数变异值;在矩阵中eij∈R(i=1,2,…,|p|,j=1,2,…,|GE|),表示i样本中j基因表达量;令矩阵中的基因集为GA=GS∪GC,样本集为p,令aij∈{0,1}(i=1,2,…,|p|,j=1,2,…,|GA|),其中为突变矩阵,当sij取值为1或i样本中j基因处于统计显著变异区域时,aij值为1,反之值为0,为了进一步整合突变矩阵和表达矩阵在突变矩阵和表达矩阵中取基因集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基因在正常样本中表达量的差异倍数用表示,其中则dij值为否则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的覆盖度对于矩阵M|p|×k中一行的互斥度,考虑这一行的离散程度,用变异系数计算每行的互斥度,每行互斥度之和为整个M|p|×k的互斥度,具体表示如公式(1)所示:
其中当趋近于0值时,对于变异系数值影响很大,所以如果M|p|×k中一行的最大权值mi≤0.5,则令该行的互斥度为放缩该行互斥度,避免该行互斥度过高对通路识别造成影响,根据和公式(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)所示:
4)设定交叉算子:
按照染色体适应度从大到小,给染色体Xi一个排名Ri,则每个染色体被选中的概率如公式(4)所示:
为了保证染色体的可行性,采用轮盘赌随机从父种群选取两条染色体,重复基因分别给子代的两条染色体,剩余基因放在一个集合,采用均匀交叉的方式,对于剩余基因的集合中连续的每对基因,随机生成一个二值数据,如果该二值数据为1,则这对基因的第一个基因放入第一条子染色体,第二个基因放入第二条子染色体,反之,第一个基因放入第二条子染色体,第二个基因放入第一条子染色体,经过一次交叉,生成两个子染色体;
5)设定变异算子:
给定一条子染色体X={x1,x2,…,xk}(xi=1,2,…,|G|),确定候选基因集合从子染色体中随机删除一个基因,得到基因集X′,将HX中基因顺序打乱,遍历前个基因,选出基因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的驱动通路。
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 CN112270952A (zh) | 2021-01-26 |
CN112270952B true 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 (7)
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 | 皇家飞利浦有限公司 | 用于使用多组学癌症谱来识别功能性患者特异性体细胞畸变的整合方法和系统 |
CN112259163A (zh) * | 2020-10-28 | 2021-01-22 | 广西师范大学 | 基于生物网络和亚细胞定位数据识别癌症驱动模块方法 |
CN113228194A (zh) * | 2018-10-12 | 2021-08-06 | 人类长寿公司 | 用于癌症基因组和临床数据综合分析的多组学搜索引擎 |
CN114023383A (zh) * | 2021-11-04 | 2022-02-08 | 广西师范大学 | 一种识别癌症驱动通路的无参非线性智能优化方法 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US10655187B2 (en) * | 2013-04-18 | 2020-05-19 | Gencurix Inc. | Genetic marker for early breast cancer prognosis prediction and diagnosis, and use thereof |
-
2020
- 2020-10-30 CN CN202011185104.6A patent/CN112270952B/zh active Active
Patent Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
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)
Title |
---|
CoDP: Cooperative Driver Pathways Discovery With Matrix Factorization and Tri-Random Walk;Ziying Yang等;《IEEE Access》;20190605;第7卷;第77738-77749页 * |
一种基于突变基因网络的癌症驱动通路识别算法;郭炳等;《计算机科学》;20180823;第45卷(第7期);第230-236,242页 * |
基于加权突变矩阵的癌症驱动通路识别模型与算法研究;朱凯;《中国优秀硕士学位论文全文数据库 医药卫生科技辑》;20210915(第9期);第E072-7页 * |
基于多组学数据识别癌症驱动通路的模型和算法;蔡齐荣等;《计算机科学》;20190915;第46卷(第9期);第310-314页 * |
基于突变基因网络的致癌驱动通路检测算法;吴昊;《计算机学报》;20170224;第41卷(第6期);第1180-1194页 * |
癌症驱动通路的启发式识别方法研究;胡灿军;《中国优秀硕士学位论文全文数据库 医药卫生科技辑》;20190315(第3期);第E072-16页 * |
通过融合基因、miRNA和代谢通路的多重关系数据联合分析挖掘的肿瘤驱动代谢途径;喵偷鹰;《https://zhuanlan.zhihu.com/p/112537964》;20200311;第1页 * |
Also Published As
Publication number | Publication date |
---|---|
CN112270952A (zh) | 2021-01-26 |
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) | 一种基于遗传禁忌和贝叶斯网络的上位性位点挖掘方法 | |
CN114420211A (zh) | 一种基于注意力机制的rna-蛋白质结合位点预测方法 | |
Ruffieux et al. | A global-local approach for detecting hotspots in multiple-response regression | |
CN104573004B (zh) | 一种基于双阶遗传计算的基因表达数据的双聚类方法 | |
CN108509764B (zh) | 一种基于遗传属性约简的古生物谱系演化分析方法 | |
Oluoch et al. | A review on RNA secondary structure prediction algorithms | |
CN112270952B (zh) | 一种识别癌症驱动通路的方法 | |
Rusman et al. | Guide RNA repertoires in the main lineages of Trypanosoma cruzi: high diversity and variable redundancy among strains | |
Sheng et al. | A survey of computational methods and databases for lncRNA-miRNA interaction prediction | |
CN114093420B (zh) | 一种基于XGBoost的DNA重组位点预测方法 | |
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 | |
CN114065933B (zh) | 一种基于人工免疫思想的未知威胁检测方法 | |
CN116994645B (zh) | 基于交互式推理网络的piRNA与mRNA靶标对的预测方法 | |
Wu et al. | Multiple sequence alignment using ga and nn | |
Wu et al. | Haplotyping a single triploid individual based on genetic algorithm | |
CN116738322A (zh) | 一种基于GA-XGBoost的日降水等级分类方法 |
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 |