CN113903395A - 一种改进粒子群优化的bp神经网络拷贝数变异检测方法与系统 - Google Patents

一种改进粒子群优化的bp神经网络拷贝数变异检测方法与系统 Download PDF

Info

Publication number
CN113903395A
CN113903395A CN202111265633.1A CN202111265633A CN113903395A CN 113903395 A CN113903395 A CN 113903395A CN 202111265633 A CN202111265633 A CN 202111265633A CN 113903395 A CN113903395 A CN 113903395A
Authority
CN
China
Prior art keywords
neural network
output
algorithm
layer
hidden layer
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.)
Pending
Application number
CN202111265633.1A
Other languages
English (en)
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.)
Liaocheng University
Original Assignee
Liaocheng 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 Liaocheng University filed Critical Liaocheng University
Priority to CN202111265633.1A priority Critical patent/CN113903395A/zh
Publication of CN113903395A publication Critical patent/CN113903395A/zh
Pending legal-status Critical Current

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
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/004—Artificial life, i.e. computing arrangements simulating life
    • G06N3/006—Artificial life, i.e. computing arrangements simulating life based on simulated virtual individual or collective life forms, e.g. social simulations or particle swarm optimisation [PSO]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/08—Learning methods
    • G06N3/084—Backpropagation, e.g. using gradient descent
    • 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

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Theoretical Computer Science (AREA)
  • General Health & Medical Sciences (AREA)
  • Molecular Biology (AREA)
  • Biophysics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Artificial Intelligence (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Evolutionary Biology (AREA)
  • Analytical Chemistry (AREA)
  • Medical Informatics (AREA)
  • Genetics & Genomics (AREA)
  • Chemical & Material Sciences (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Biomedical Technology (AREA)
  • Computational Linguistics (AREA)
  • Data Mining & Analysis (AREA)
  • Evolutionary Computation (AREA)
  • Computing Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Software Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明本文提出了一种基于bp神经网络和自适应粒子群优化算法的拷贝数变异检测方法，流程如图1所示。该方法有三个特点：(1)选取了影响拷贝数变异检测的多个特征，并充分考虑了多个特征的相互作用；(2)它使用粒子群算法优化BP神经网络的权值和阈值，解决了BP神经网络已陷入局部最优的问题；(3)它使用了一个变异算子，解决了粒子群算法搜索精度低的问题。

Description

一种改进粒子群优化的BP神经网络拷贝数变异检测方法与 系统
技术领域
本发明属于生物信息学领域，具体涉及一种基于自适应粒子群算法优化的BP神经网络拷贝数变异检测方法。
技术背景
拷贝数变异是指DNA拷贝数偏离了正常数值范围，一般表现为扩展或者缺失。人类基因组中包含的拷贝数区域约占基因序列的1％，因此CNV与人类疾病的发生和发展息息相关，例如自闭症和阿兹海默症。这些疾病诊断的关键一步就是准确检测目标基因组的拷贝数变异，与之相应的基因组测序技术在迅速发展。下一代测序技术由于测序速度快、数据吞吐量高、成本低等优势，迅速发展成为主流的测序技术。
下一代测序技术又称为高通量测序技术，可以产生大量的测序数据，是生命科学领域迅猛发展的分支之一，也使得生物信息学进入了大数据时代。由此也产生了大量的基于下一代测序数据的拷贝数变异检测算法。
基于下一代测序技术的拷贝数变异检测方法主要分为四类：基于双端比对，基于分裂读段，基于序列组装，基于读深(Read Depth,RD)。基于双端比对的方法可以有效的从不一致映射的成对读取中识别扩增、缺失和倒置的CNVs，但是在低复杂度区域中的检测能力较差。基于分裂读段的方法通过分析断点位置检测CNVs，适用于检测读段长度较短的结构变异。基于序列组装的方法通过比较组装后的序列与参考基因组的差异检测CNVs，但在序列长度较短时不利于拼接。
在这四类策略中，使用最广的是基于RD策略设计的方法。由于基因组不同位置上的读深信号存在差异，基于RD的方法通过分析这种差异检测拷贝数变异。早期基于RD的方法的方法只利用RD信号这一个特征值，这些方法包括G-CNV，CNV-seq，CNAseg，CNVnator，BIC-seq，Control-FREEC，CNAnorm等。其中，Control-FREEC使用RD信号作为输入，构建和分析拷贝数和B等位基因频率分布图谱，以检测基因组的拷贝数状态。然而上述方法只能检测出显著性的CNV，且不能准确定位CNV的断点位置。近年来，越来越多的拷贝数变异检测方法开始考虑除了RD信号以外的影响因素。这些方法包括m-HMM，CNV-RF，CoNVaDING，seqCNV，ICopyDAV，CONDEL，modSaRa2，BagGMM，CNV_IFTV，SM-RCNV，CNV-LOF，AITAC，DINTD，CRSCNV，DCC，RKDOSCNV，IndivCNV等。
基于RD的方法主要按照以下几个步骤检测CNV：(1)原始基因组文件的提取和处理(2)参考基因组比对，(3)GC含量偏差校正，(4)RD信号对比和计算，(5)建立统计模型以检测CNV。其中，CNV-IFTV建立了一个包括隔离树和全变分的统计模型，这个模型首先为每个基因组bin设置一个分数，然后使用隔离树算法对基因组bin进行训练。CNV-LOF没有从全局的角度考虑RD信号，而是对RD信号进行局部观察，并使用离群值因子来反映CNV。RKDOSCNV采用基于扩展的最近邻算法和核密度估计方法评估RD信号的核密度分布情况，并选择合理的阈值来筛选CNV。CRSCNV使用包含交叉模型的统计策略，并将基因组箱拆分成多个部分进行分析。IndivCNV使用全变分模型对RD信号进行平滑处理，然后利用潜变量模型构建权值矩阵，最后对信号进行分层，根据分层矩阵能量谱在每层的占比鉴别个体性拷贝数变异。然而上述方法在低纯度和低肿瘤覆盖率的样本中并未取得十分显著的效果。
近年来，机器学习被广泛应用于分类和预测问题，取得了显著的效果。例如，图形神经网络被用于预测新冠肺炎疫情中单个传染性个体在可能不再完全易感的人群中传播新感染的数量，有效的估计新冠肺炎的流行程度。MFCNV将多个特征作为BP神经网络模型的输入以预测拷贝数变异的状态。但是bp神经网络仍有很多缺点，例如陷入局部最优，过拟合。
基于上述研究，本文提出了一种基于bp神经网络和自适应粒子群优化算法的拷贝数变异检测方法。该方法有三个特点：(1)选取了影响拷贝数变异检测的多个特征，并充分考虑了多个特征的相互作用；(2)它使用粒子群算法优化BP神经网络的权值和阈值，解决了BP神经网络已陷入局部最优的问题；(3)它使用了一个变异算子，解决了粒子群算法搜索精度低的问题。
发明内容
本发明的目的在于解决上述现有技术中存在的难题，提供一种基于自适应粒子群算法优化的BP神经网络拷贝数变异检测方法，提升检测性能，提高预测精度。本发明是通过以下技术方案实现的：
一种基于自适应粒子群算法优化的BP神经网络拷贝数变异检测方法，包括：
S1：研究了基于下一代测序数据的拷贝数变异检测问题；
S2：对数据进行预处理；
S3：定义与CNV相关的四个特征值，计算并归一化；
S4：利用四个特征值构建BP神经网络；
S5：使用自适应粒子群算法优化BP神经网络的权值和阈值，训练BP神经网络，根据训练完成的BP神经网络模型预测每个窗的CNV状态(正常、扩增或缺失)。
S6：所提出的算法对研究的问题进行了有效的验证。
所述S1拷贝数变异检测问题中，采用下一代测序数据(即高通量测序数据)，与第一代测序技术合成终止后再测序不同的是，高通量测序技术加入了可逆终止末端，因此可边合成边测序，极大的降低了测序成本，缩短了测序用时，提高了测序通量。
所述S2对数据进行预处理时，采用经典的比对方法BWA处理输入的两个文件(即FASTA格式的参考文件和FASTQ格式的测序样本文件)，产生一个BAM格式的对齐文件,然后使用SAMtools软件从BAM文件中获取读取计数(RC)配置文件。这个RC配置文件的模板是参考基因组，其中大部分位置是常规碱基(“A”、“T”、“C”和“G”)，而其中一小部分位置无法确定。未确定的位置通常用字母“N”填充，因此无法与这些区域匹配测序读数。为了获得合理的RC文件，舍弃了参考基因组中包含“N”的区域。
所述S3定义了与CNV相关的四个特征值，即RD、GC含量、基因组中相邻位点之间的相关性和碱基质量。RD为短读映射后的密度，是NGS数据中的重要信号。GC含量是鸟嘌呤和胞嘧啶占DNA碱基总数的比值。由于测序深度的原因，GC含量容易出现偏差，影响CNV的检测。虽然很多方法预先校正了GC含量偏差，但没有消除偏差对CNV检测的影响。因此，PBCNV将GC含量作为分析和检测的一个特征。基因组中相邻位点之间的相关性反映了CNV区域与相邻区域的关联度。碱基质量反映了测序的准确性，也是影响CNV检测的因素之一。由四个特征组成的特征向量如公式所示：
Bi＝(Ri,Qi,Ci,Gi)
其中Ri值代表RD值，Qi代表碱基质量，Ci代表基因组中相邻位点之间的相关性，Gi代表GC含量，Bi代表第i个滑动窗口的特征向量。
为了平衡特征值之间的数值差异，我们对特征值进行归一化。以RD值为例，归一化方法如公式所示：
Figure BDA0003326889280000041
其中，Rmax代表RD最大值，Rmin代表RD最小值。
所述S4利用四个特征值构建BP神经网络，BP神经网络由输入层、隐含层和输出层组成。不同的层完全互连，并且可能存在一个或多个隐藏层。在PBCNV中，采用了三层神经网络。BP神经网络的结构图如图2所示，输入层节点数由特征值的个数决定。PBCNV选择四个特征值，因此网络有四个输入层节点。输出层的节点数由CNV的类别决定。CNV有三种类别(正常、扩增或缺失)，因此网络在输出层有三个节点。一般认为增加隐层数可以减少网络误差，提高精度，但也会使网络变得更加复杂，从而增加网络训练时间，也可能导致过拟合。因此，可以通过增加隐藏层中的节点数来获得较低的误差。在PBCNV中，隐含层节点的数量是15个，这是经过大量实验确定的。
输入层和隐含层之间的激活函数采用Sigmoid函数，隐含层和输出层之间的函数采用线性函数。Sigmoid函数如公式所示：
Figure BDA0003326889280000051
其中x表示特征向量。学习率衡量权重更新的大小。如果学习率较低，则训练进展缓慢；如果学习率较高，则可能会给损失函数带来不良后果。在PBCNV中，学习率被设置为0.1。
所述S5使用自适应粒子群算法优化BP神经网络的权值和阈值，训练BP神经网络，根据训练完成的BP神经网络模型预测每个窗的CNV状态(正常、扩增或缺失)。所提PBCNV算法由两部分组成：BP神经网络部分和自适应粒子群优化部分。
(1)自适应粒子群优化算法部分。APSO算法首先在解空间中初始化一组粒子，每个粒子代表一个潜在的最优解。粒子特征由位置、速度和适应值表示。适应值由适应度函数计算，其值的好坏表示粒子的质量。粒子在解空间中运动，通过跟踪个体极值和种群极值来更新个体位置。每次粒子更新其位置时，都会计算一个适应值。通过比较新粒子与个体极值和群体极值的适应值，更新个体极值和群体极值的位置。具体步骤如下：
步骤1，种群初始化
步骤2，寻找初始极值
步骤3，更新粒子
步骤4，变异操作
所述步骤1是这样实现的：
随机初始化粒子位置和粒子速度，根据适应度函数计算粒子适应值。取预期输出和预测输出之间的绝对误差作为适应值，如公式所示。
Figure BDA0003326889280000061
其中，m表示输出节点的总数，n表示系数，di和pi分别表示神经网络的第i个节点的期望输出和预测输出。
所述步骤2是这样实现的：
根据初始粒子适应度值找到个体极值和群体极值。
所述步骤3是这样实现的：
根据公式更新粒子的位置和速度，然后根据新粒子的适应度值更新个体极值和群体极值。
Figure BDA0003326889280000062
Figure BDA0003326889280000063
其中ω是权重，D是空间维度，d＝1,2,...,D；I是粒子数，i＝1,2,...,I；k是当前迭代次数；Vid是粒子速度，c1和c2是非负常数，称为加速因子，r1和r2是分布在[0，1]之间的随机数。
所述步骤4是这样实现的：
为了提高算法找到更优值的可能性，使用变异算子以一定的概率初始化粒子。突变概率为0.1。
(2)BP神经网络部分。将APSO算法优化后的初始权值和阈值代入BP神经网络进行训练。算法二给出了BP神经网络的实现过程，并采用粒子群优化算法对个体进行适应度测试。对不符合适应度标准的个体进行变异，并对BP神经网络训练进行重新训练，直到个体达到适应度标准或精度要求。具体步骤如下：
步骤1，神经网络初始化
步骤2，计算隐含层的输出
步骤3，计算输出层的输出
步骤4，计算误差
步骤5，更新权值
步骤6，更新阈值
所述步骤1是这样实现的：
初始化输入层、隐含层和输出层的节点数。输入层、隐含层和输出层的节点数分别为n,l,m。输入层与隐含层之间的权值为ωij，隐含层与输出层之间的权值为ωjk，输入层与隐含层之间的偏差为aj，隐含层与输出层之间的偏差为bk，学习率为η。
所述步骤2是这样实现的：
根据输入变量X，输入层和隐含层之间的连接权值ωij以及隐含层阈值a，计算隐含层输出H。
Figure BDA0003326889280000071
其中，l为隐含层节点数，s为隐含层激励函数。
所述步骤3是这样实现的：
根据隐含层输出H，隐含层与输出层之间的权值为ωjk和阈值b计算神经网络预测输出O。
Figure BDA0003326889280000072
所述步骤4是这样实现的：
根据网络预测输出O和期望输出Y，计算网络预测误差e。
Figure BDA0003326889280000073
所述步骤5是这样实现的：
根据网络预测误差e更新网络连接权值ωij，ωjk。
Figure BDA0003326889280000074
ωjk＝ωjk+ηHjek
所述步骤6是这样实现的：
根据网络预测误差e更新网络节点阈值a，b。
Figure BDA0003326889280000081
bk＝bk+ηek
所述S6使用所提出的算法对研究的问题进行了有效的验证。在临床实验前，通过仿真实验验证了算法的合理性和可靠性。为了将PBCNV应用于临床诊断，用实际数据样本验证了该方法的有效性。我们将HICA算法与其他当前流行的算法进行了比较，包括CNVnator，FREEC，iCopyDAV，ReadDepth，MFCNV，PBCNV，来验证提出的算法的优势。
附图说明
图1.PBCNV工作流程图。
图2.BP神经网络结构图。
图3.APSO优化BP神经网络流程图。
图4.仿真样本中PBCNV及其他五种方法性能表现图。
图5.仿真样本中PBCNV及其他五种方法边界偏差图。
图6.仿真样本中PBCNV及其他五种方法ROC曲线图。
图7.真实样本中PBCNV及其他五种方法性能表现图。
具体实施方式
下面结合附图对本发明做进一步详细描述：
拷贝数变异是一种基因组结构性变异，是人类疾病特别是癌症诊断的重要因素之一。下一代测序技术的出现，带来了大量的测序数据，同时产生了相应的拷贝数变异检测算法。然而，现有检测算法受测序误差和GC含量偏差的影响，在低样本纯度和低肿瘤覆盖率样本中的检测能力不足。针对这种问题，我们提出了一种基于BP神经网络和自适应粒子群优化算法的拷贝数变异检测方法，PBCNV。该方法具有三个特点：(1)选取了影响拷贝数变异检测的多个特征，并充分考虑了多个特征的相互作用；(2)针对BP神经网络模型易陷入局部最优的缺点，使用粒子群优化算法优化BP神经网络的权值和阈值；(3)针对粒子群算法搜索精度较低的缺点，引入变异操作，提高算法的检测能力。
1.PBCNV工作流程。
PBCNV是一种基于RD策略的设计方法，工作流程图如图1所示。首先，算法从比较样本和参考基因组开始，使用BWA算法比较FASTQ格式和FASTA格式的文件。接下来，使用Samtools软件获得读取计数(RC)文件，该文件用于计算RD值。然后，PBCNV通过以下步骤识别CNV：(1)定义和计算与CNV检测相关的4个特征，包括RD、GC含量、基因组相邻位点间的相关性和碱基质量；(2)确定BP神经网络的结构；(3)使用APSO算法优化BP神经网络；(4)BP神经网络的训练；(5)CNV的预测。
2.特征的定义和计算。
PBCNV选择四个特征来识别CNV，即RD、GC含量、基因组中相邻位点之间的相关性和碱基质量。RD为短读映射后的密度，是NGS数据中的重要信号。GC含量是DNA碱基中鸟嘌呤和胞嘧啶的比值。由于测序深度的原因，GC含量容易出现偏差，影响CNV的检测。几种方法预先校正了GC含量偏差，但没有消除偏差对CNV检测的影响。因此，PBCNV将GC含量作为分析和检测的一个特征。基因组中相邻位点之间的相关性也反映了CNV区域。碱基质量反映了测序的准确性，也是影响CNV检测的因素之一。特征的计算方法可参考MFCNV。由四个特征组成的特征向量如公式所示：
Bi＝(Ri,Qi,Ci,Gi)
其中Ri值代表RD值，Qi代表碱基质量，Ci代表基因组中相邻位点之间的相关性，Gi代表GC含量，Bi代表第i个滑动窗口的特征向量。
为了平衡特征值之间的数值差异，我们对特征值进行归一化。以RD值为例，归一化方法如公式所示：
Figure BDA0003326889280000101
其中，Rmax代表RD最大值，Rmin代表RD最小值。
3.BP神经网络的构建。
BP神经网络由输入层、隐含层和输出层组成。不同的层完全互连，并且可能存在一个或多个隐藏层。在PBCNV中，采用了三层神经网络。BP神经网络的结构图如图2所示，输入层节点数由特征值的个数决定。PBCNV选择四个特征值，因此网络有四个输入层节点。输出层的节点数由CNV的类别决定。CNV有三种类别(正常、扩增或缺失)，因此网络在输出层有三个节点。一般认为增加隐层数可以减少网络误差，提高精度，但也会使网络变得更加复杂，从而增加网络训练时间，也可能导致过拟合。因此，可以通过增加隐藏层中的节点数来获得较低的误差。在PBCNV中，隐含层节点的数量是15个，这是经过大量实验确定的。输入层和隐含层之间的激活函数采用Sigmoid函数，隐含层和输出层之间的函数采用线性函数。Sigmoid函数如公式所示：
Figure BDA0003326889280000111
其中x表示特征向量。学习率衡量权重更新的大小。如果学习率较低，则训练进展缓慢；如果学习率较高，则可能会给损失函数带来不良后果。在PBCNV中，学习率被设置为0.1。
4.使用APSO算法优化BP神经网络。
粒子群优化(PSO)算法是计算智能领域中的一种群智能优化算法。粒子群算法的灵感来源于鸟类的行为特点，用于求解优化问题。算法中的每个粒子代表问题的一个潜在解，每个粒子对应一个由适应度函数确定的适应值。粒子的速度决定粒子运动的方向和距离，并根据粒子自身和其他粒子的运动经验动态调整速度，从而实现可解空间中的个体优化。APSO优化BP神经网络的流程图如图3所示。
所提PBCNV算法由两部分组成：BP神经网络部分和自适应粒子群优化部分。
(1)自适应粒子群优化算法部分。APSO算法首先在解空间中初始化一组粒子，每个粒子代表一个潜在的最优解。粒子特征由位置、速度和适应值表示。适应值由适应度函数计算，其值的好坏表示粒子的质量。粒子在解空间中运动，通过跟踪个体极值和种群极值来更新个体位置。每次粒子更新其位置时，都会计算一个适应值。通过比较新粒子与个体极值和群体极值的适应值，更新个体极值和群体极值的位置。具体步骤如下：
步骤1，种群初始化
随机初始化粒子位置和粒子速度，根据适应度函数计算粒子适应值。取预期输出和预测输出之间的绝对误差作为适应值，如公式所示。
Figure BDA0003326889280000112
其中，m表示输出节点的总数，n表示系数，di和pi分别表示神经网络的第i个节点的期望输出和预测输出。
步骤2，寻找初始极值
根据初始粒子适应度值找到个体极值和群体极值。
步骤3，更新粒子
根据公式更新粒子的位置和速度，然后根据新粒子的适应度值更新个体极值和群体极值。
Figure BDA0003326889280000121
Figure BDA0003326889280000122
其中ω是权重，D是空间维度，d＝1,2,...,D；I是粒子数，i＝1,2,...,I；k是当前迭代次数；Vid是粒子速度，c1和c2是非负常数，称为加速因子，r1和r2是分布在[0，1]之间的随机数。
步骤4，变异操作
为了提高算法找到更优值的可能性，使用变异算子以一定的概率初始化粒子。突变概率为0.1。
(2)BP神经网络部分。将APSO算法优化后的初始权值和阈值代入BP神经网络进行训练。算法二给出了BP神经网络的实现过程，并采用粒子群优化算法对个体进行适应度测试。对不符合适应度标准的个体进行变异，并对BP神经网络训练进行重新训练，直到个体达到适应度标准或精度要求。具体步骤如下：
步骤1，神经网络初始化
初始化输入层、隐含层和输出层的节点数。输入层、隐含层和输出层的节点数分别为n,l,m。输入层与隐含层之间的权值为ωij，隐含层与输出层之间的权值为ωjk，输入层与隐含层之间的偏差为aj，隐含层与输出层之间的偏差为bk，学习率为η。
步骤2，计算隐含层的输出
根据输入变量X，输入层和隐含层之间的连接权值ωij以及隐含层阈值a，计算隐含层输出H。
Figure BDA0003326889280000131
其中，l为隐含层节点数，s为隐含层激励函数。
步骤3，计算输出层的输出
根据隐含层输出H，隐含层与输出层之间的权值为ωjk和阈值b计算神经网络预测输出O。
Figure BDA0003326889280000132
步骤4，计算误差
根据网络预测输出O和期望输出Y，计算网络预测误差e。
Figure BDA0003326889280000133
步骤5，更新权值
根据网络预测误差e更新网络连接权值ωij，ωjk。
Figure BDA0003326889280000134
ωjk＝ωjk+ηHjek
步骤6，更新阈值
根据网络预测误差e更新网络节点阈值a，b。
Figure BDA0003326889280000135
bk＝bk+ηek
5.使用PBCNV预测CNV
将训练好的PBCNV算法应用于每个测试数据，对CNV进行预测。每个输入层节点对应一个特征，每个输出层节点对应一个CNV状态。对于每个基因组滑动窗口，PBCNV算法输出三个反映CNV状态的概率值。对应于三个值中最大值的CNV状态被认为是当前基因组框的CNV状态。三种概率值的预测结果提高了算法的预测精度。
6.实验分析
6.1仿真实验
在临床实验前，通过仿真实验验证了算法的合理性和可靠性。为了防止细胞污染等外界因素造成的误差，使用不同纯度和覆盖深度的样本作为训练数据。PBCNV的参数如表1所示。
表1 PBCNV方法参数设置
Figure BDA0003326889280000141
针对样品纯度低、肿瘤覆盖率低的特点，采用IntSIM软件生成模拟数据[50]。肿瘤纯度设置为0.2、0.3和0.4，覆盖深度设置为4x和6x。每个配置生成50个样本，总共300个样本。然后将PBCNV与现有的5种方法(FREEC、CNVnator、iCopyDAV、CNAnorm、MFCNV)进行比较。以灵敏度、精密度和F1得分作为参考指标。每项指标平均抽取50份样本。六种方法的比较结果如图4所示。
从图4可以看出，随着样本纯度和覆盖深度的降低，算法的检测能力也随之降低。当样品纯度为0.4，覆盖深度为6倍时，CNVnator的灵敏度为0.58；当样品纯度为0.2，覆盖深度为6倍时，CNVnator的灵敏度为0.31。当样本纯度为0.4，覆盖深度为6x时，FREEC的F1-score为0.72；当样本纯度为0.2，覆盖深度为6x时，FREEC的F1-score为0.57。与MFCNV相比，PBCNV的灵敏度、精度和F1-score均有不同程度的提高。在六种类型的样本中，PBCNV的F1-Score达到了最高值，并且PBCNV的性能优于其他五种等价算法。PBCNV的精度在5个样本中达到最高值，CNVnator在1个样本中达到最高值。PBCNV的灵敏度在3个样本中达到最高值，CNAnorm在3个样本中达到最高值。
边界偏差的检测也是衡量算法性能的指标之一。准确检测CNV的边界有利于准确定位CNV的突变位点，便于靶向药物的选择和使用。边界偏差越大，算法的检测能力越差。六种方法的边界偏差比较结果如图5所示。
从图5可以看出，随着样本纯度和肿瘤覆盖率的降低，六种方法的边界偏差增大。在这六类样本中，PBCNV算法的边界偏差小于MFCNV算法和其他四种算法。这表明PBCNV在CNV检测的准确性方面具有明显的优势。
另外，ROC曲线可以反映模型选择不同阈值时的假阳性率(FPR)和真阳性率(TPR)的变化趋势。ROC曲线接近左上角，表明算法的检测效果较好。六种方法(CNVnator、FREEC、GROM_RD、iCopyDAV、MFCNV、PBCNV)的ROC曲线如图6所示，随着覆盖率的降低，ROC曲线趋于右下角。PBCNV的ROC曲线优于其他五种算法，说明在选取不同的阈值时，PBCNV的检测能力仍然较好。
PBCNV比其他算法具有更好的检测能力，主要表现在：(1)训练样本数量大，需要对大量反映拷贝数状态的特征进行训练和分析；(2)单个特征的分析有限，因此PBCNV提取四个特征提高CNV状态预测的准确性；(3)PBCNV将BP神经网络与PSO算法相结合，增强BP神经网络的寻优能力；(4)利用变异算子解决PSO算法搜索精度低的问题，进一步提高了BP神经网络的寻优能力。
6.2真实样本引用
为了将PBCNV应用于临床诊断，用实际数据样本验证了该方法的有效性。实际样本数据来自1000基因组计划(http://www.internationalgenome.org/)，名为NA19238、NA19239和NA19240。PBCNV和其他五种方法(CNV、FREEC、iCopyDAV、ReadDepth、MFCNV)检测到的CNV数量如表2所示。基因组变异数据库(http://dgv.tcag.ca/)确认了每个真实样本中的CNV数量。
表2 PBCNV和五种对等方法检测到的CNV数量
Figure BDA0003326889280000161
从表2可以看出，PBCNV在NA19239和NA19240中检测到的CNV数量最多。在NA19238中，CNV的检出数量仅次于CNVnator。这表明PBCNV在实际样品应用中具有较强的检测能力，已达到临床试验水平。图7显示了这五种方法在实际数据上的性能比较结果。
从图7可以看出，在NA19238中，PBCNV的灵敏度略低于CNVnator，但PBCNV在其他两个样本中的灵敏度最高。PBCNV的F1-score高于其他五种算法，平均值为0.84，MFCNV的F1-score为0.79。结果表明，考虑了PBCNV的灵敏度和精度。
6.3结果分析
针对现有算法在样本纯度低、肿瘤覆盖率低的情况下检测能力低的问题，提出了一种基于BP神经网络和自适应粒子群优化算法的拷贝数变异检测方法。首先，PBCNV在不校正GC含量偏差的情况下提取和计算四个特征值，以避免误差。然后，用PBCNV构造了一个三层BP神经网络，以特征值为输入，CNV类型为输出。然后，利用粒子群优化算法优化BP神经网络的权值和阈值，并引入变异算子进一步提高算法的检测能力。最后，PBCNV利用大量的训练数据对模型进行训练，实现了CNV的预测。实验结果表明，PBCNV在灵敏度、精度、F1-score和边界偏差等方面具有明显的优势。

Claims (7)

1.一种改进粒子群优化的BP神经网络拷贝数变异检测方法与系统，其特征在于：所述方法包括：
S1：研究了基于下一代测序数据的拷贝数变异检测问题；
S2：对数据进行预处理；
S3：定义与CNV相关的四个特征值，计算并归一化；
S4：利用四个特征值构建BP神经网络；
S5：使用自适应粒子群算法优化BP神经网络的权值和阈值，训练BP神经网络，根据训练完成的BP神经网络模型预测每个窗的CNV状态(正常、扩增或缺失),
S6：所提出的算法对研究的问题进行了有效的验证。
2.根据权力要求1所述的一种基于自适应粒子群算法优化的BP神经网络拷贝数变异检测方法，其特征在于：
所述S1拷贝数变异检测问题中，采用下一代测序数据(即高通量测序数据)，与第一代测序技术合成终止后再测序不同的是，高通量测序技术加入了可逆终止末端，因此可边合成边测序，极大的降低了测序成本，缩短了测序用时，提高了测序通量。
3.根据权力要求1所述的一种基于自适应粒子群算法优化的BP神经网络拷贝数变异检测方法，其特征在于：所述S2对数据进行预处理时，采用经典的比对方法BWA处理输入的两个文件(即FASTA格式的参考文件和FASTQ格式的测序样本文件)，产生一个BAM格式的对齐文件,然后使用SAMtools软件从BAM文件中获取读取计数(RC)配置文件；这个RC配置文件的模板是参考基因组，其中大部分位置是常规碱基(“A”、“T”、“C”和“G”)，而其中一小部分位置无法确定；未确定的位置通常用字母“N”填充，因此无法与这些区域匹配测序读数；为了获得合理的RC文件，舍弃了参考基因组中包含“N”的区域。
4.根据权力要求1所述的一种基于自适应粒子群算法优化的BP神经网络拷贝数变异检测方法，其特征在于：所述S3定义了与CNV相关的四个特征值，即RD、GC含量、基因组中相邻位点之间的相关性和碱基质量；RD为短读映射后的密度，是NGS数据中的重要信号；GC含量是鸟嘌呤和胞嘧啶占DNA碱基总数的比值；由于测序深度的原因，GC含量容易出现偏差，影响CNV的检测；虽然很多方法预先校正了GC含量偏差，但没有消除偏差对CNV检测的影响；因此，PBCNV将GC含量作为分析和检测的一个特征；基因组中相邻位点之间的相关性反映了CNV区域与相邻区域的关联度；碱基质量反映了测序的准确性，也是影响CNV检测的因素之一；由四个特征组成的特征向量如公式(1)所示：
Bi＝(Ri,Qi,Ci,Gi) (1)
其中Ri值代表RD值，Qi代表碱基质量，Ci代表基因组中相邻位点之间的相关性，Gi代表GC含量，Bi代表第i个滑动窗口的特征向量；为了平衡特征值之间的数值差异，我们对特征值进行归一化,以RD值为例，归一化方法如公式(2)所示：
Figure FDA0003326889270000021
其中，Rmax代表RD最大值，Rmin代表RD最小值。
5.根据权力要求1所述的一种基于自适应粒子群算法优化的BP神经网络拷贝数变异检测方法，其特征在于：所述S4利用四个特征值构建BP神经网络，BP神经网络由输入层、隐含层和输出层组成；不同的层完全互连，并且可能存在一个或多个隐藏层；在PBCNV中，采用了三层神经网络；BP神经网络的结构图如图2所示，输入层节点数由特征值的个数决定；PBCNV选择四个特征值，因此网络有四个输入层节点；输出层的节点数由CNV的类别决定；CNV有三种类别(正常、扩增或缺失)，因此网络在输出层有三个节点；一般认为增加隐层数可以减少网络误差，提高精度，但也会使网络变得更加复杂，从而增加网络训练时间，也可能导致过拟合；因此，可以通过增加隐藏层中的节点数来获得较低的误差；在PBCNV中，隐含层节点的数量是15个，这是经过大量实验确定的；输入层和隐含层之间的激活函数采用Sigmoid函数，隐含层和输出层之间的函数采用线性函数；Sigmoid函数如公式(3)所示：
Figure FDA0003326889270000031
其中x表示特征向量。学习率衡量权重更新的大小；如果学习率较低，则训练进展缓慢；如果学习率较高，则可能会给损失函数带来不良后果；在PBCNV中，学习率被设置为0.1。
6.根据权力要求1所述的一种基于自适应粒子群算法优化的BP神经网络拷贝数变异检测方法，其特征在于：所述S5使用自适应粒子群算法优化BP神经网络的权值和阈值，训练BP神经网络，根据训练完成的BP神经网络模型预测每个窗的CNV状态(正常、扩增或缺失)；所提PBCNV算法由两部分组成：BP神经网络部分和自适应粒子群优化部分；
(1)自适应粒子群优化算法部分，APSO算法首先在解空间中初始化一组粒子，每个粒子代表一个潜在的最优解。粒子特征由位置、速度和适应值表示；适应值由适应度函数计算，其值的好坏表示粒子的质量；粒子在解空间中运动，通过跟踪个体极值和种群极值来更新个体位置；每次粒子更新其位置时，都会计算一个适应值；通过比较新粒子与个体极值和群体极值的适应值，更新个体极值和群体极值的位置；具体步骤如下：
步骤1，种群初始化
步骤2，寻找初始极值
步骤3，更新粒子
步骤4，变异操作
所述步骤1是这样实现的：
随机初始化粒子位置和粒子速度，根据适应度函数计算粒子适应值；取预期输出和预测输出之间的绝对误差作为适应值，如公式(4)所示：
Figure FDA0003326889270000041
其中，m表示输出节点的总数，n表示系数，di和pi分别表示神经网络的第i个节点的期望输出和预测输出；
所述步骤2是这样实现的：
根据初始粒子适应度值找到个体极值和群体极值；
所述步骤3是这样实现的：
根据公式(5)和公式(6)更新粒子的位置和速度，然后根据新粒子的适应度值更新个体极值和群体极值；
Figure FDA0003326889270000042
Figure FDA0003326889270000043
其中ω是权重，D是空间维度，d＝1,2,...,D；I是粒子数，i＝1,2,...,I；k是当前迭代次数；Vid是粒子速度，c1和c2是非负常数，称为加速因子，r1和r2是分布在[0，1]之间的随机数；
所述步骤4是这样实现的：
为了提高算法找到更优值的可能性，使用变异算子以一定的概率初始化粒子。突变概率为0.1；
(2)BP神经网络部分，将APSO算法优化后的初始权值和阈值代入BP神经网络进行训练；算法二给出了BP神经网络的实现过程，并采用粒子群优化算法对个体进行适应度测试；对不符合适应度标准的个体进行变异，并对BP神经网络训练进行重新训练，直到个体达到适应度标准或精度要求；具体步骤如下：
步骤1，神经网络初始化
步骤2，计算隐含层的输出
步骤3，计算输出层的输出
步骤4，更新权值
步骤5，更新阈值
所述步骤1是这样实现的：
初始化输入层、隐含层和输出层的节点数；输入层、隐含层和输出层的节点数分别为n,l,m；输入层与隐含层之间的权值为ωij，隐含层与输出层之间的权值为ωjk，输入层与隐含层之间的偏差为aj，隐含层与输出层之间的偏差为bk，学习率为η；
所述步骤2是这样实现的：
根据输入变量X，输入层和隐含层之间的连接权值ωij以及隐含层阈值a，计算隐含层输出H；
Figure FDA0003326889270000051
其中，l为隐含层节点数，s为隐含层激励函数；
所述步骤3是这样实现的：
根据隐含层输出H，隐含层与输出层之间的权值为ωjk和阈值b计算神经网络预测输出O；
Figure FDA0003326889270000052
所述步骤4是这样实现的：
根据网络预测输出O和期望输出Y，计算网络预测误差e；
Figure FDA0003326889270000053
所述步骤5是这样实现的：
根据网络预测误差e更新网络连接权值ωij，ωjk；
Figure FDA0003326889270000054
ωjk＝ωjk+ηHjek (11)
所述步骤6是这样实现的：
根据网络预测误差e更新网络节点阔值a，b；
Figure FDA0003326889270000061
bk＝bk+ηek (13)
7.根据权力要求1所述的一种基于自适应粒子群算法优化的BP神经网络拷贝数变异检测方法，其特征在于：所述S6使用所提出的算法对研究的问题进行了有效的验证；在临床实验前，通过仿真实验验证了算法的合理性和可靠性；为了将PBCNV应用于临床诊断，用实际数据样本验证了该方法的有效性；我们将HICA算法与其他当前流行的算法进行了比较，包括CNVnator，FREEC，iCopyDAV，ReadDepth，MFCNV，PBCNV，来验证提出的算法的优势。
CN202111265633.1A 2021-10-28 2021-10-28 一种改进粒子群优化的bp神经网络拷贝数变异检测方法与系统 Pending CN113903395A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111265633.1A CN113903395A (zh) 2021-10-28 2021-10-28 一种改进粒子群优化的bp神经网络拷贝数变异检测方法与系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111265633.1A CN113903395A (zh) 2021-10-28 2021-10-28 一种改进粒子群优化的bp神经网络拷贝数变异检测方法与系统

Publications (1)

Publication Number Publication Date
CN113903395A true CN113903395A (zh) 2022-01-07

Family

ID=79027634

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111265633.1A Pending CN113903395A (zh) 2021-10-28 2021-10-28 一种改进粒子群优化的bp神经网络拷贝数变异检测方法与系统

Country Status (1)

Country Link
CN (1) CN113903395A (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114863242A (zh) * 2022-04-26 2022-08-05 北京拙河科技有限公司 一种面向图像识别的深度学习网络的优化方法及系统
CN115630101A (zh) * 2022-10-24 2023-01-20 淮阴工学院 水文参数智能化监控与水资源大数据管理系统
CN117095744A (zh) * 2023-08-21 2023-11-21 上海信诺佰世医学检验有限公司 一种基于单样本高通量转录组测序数据的拷贝数变异检测方法
CN117724446A (zh) * 2023-12-14 2024-03-19 广州智业节能科技有限公司 一种预警监控系统、方法及装置

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114863242A (zh) * 2022-04-26 2022-08-05 北京拙河科技有限公司 一种面向图像识别的深度学习网络的优化方法及系统
CN114863242B (zh) * 2022-04-26 2022-11-29 北京拙河科技有限公司 一种面向图像识别的深度学习网络的优化方法及系统
CN115630101A (zh) * 2022-10-24 2023-01-20 淮阴工学院 水文参数智能化监控与水资源大数据管理系统
CN115630101B (zh) * 2022-10-24 2023-10-20 淮阴工学院 水文参数智能化监控与水资源大数据管理系统
CN117095744A (zh) * 2023-08-21 2023-11-21 上海信诺佰世医学检验有限公司 一种基于单样本高通量转录组测序数据的拷贝数变异检测方法
CN117724446A (zh) * 2023-12-14 2024-03-19 广州智业节能科技有限公司 一种预警监控系统、方法及装置

Similar Documents

Publication Publication Date Title
CN113903395A (zh) 一种改进粒子群优化的bp神经网络拷贝数变异检测方法与系统
CN110070141B (zh) 一种网络入侵检测方法
US7085690B2 (en) Unsupervised machine learning-based mathematical model selection
CN112581262A (zh) 一种基于鲸鱼算法优化lvq神经网络的欺诈行为检测方法
CN109637579B (zh) 一种基于张量随机游走的关键蛋白质识别方法
CN112215259B (zh) 基因选择方法和装置
CN110060738B (zh) 基于机器学习技术预测细菌保护性抗原蛋白的方法及系统
CN113488104B (zh) 基于局部和全局的网络中心性分析的癌症驱动基因预测方法及系统
CN114841257A (zh) 一种基于自监督对比约束下的小样本目标检测方法
CN114566216B (zh) 一种基于注意力机制的剪接位点预测及解释性方法
CN111145830A (zh) 基于网络传播的蛋白质功能预测方法
WO2019144798A1 (zh) 染色质相互作用差异的分析方法和系统
CN113555062A (zh) 一种用于基因组碱基变异检测的数据分析系统及分析方法
Zaman et al. Codon based back propagation neural network approach to classify hypertension gene sequences
CN115393632A (zh) 一种基于进化多目标神经网络架构构造的图像分类方法
CN117371511A (zh) 图像分类模型的训练方法、装置、设备及存储介质
Bai et al. A unified deep learning model for protein structure prediction
CN114974400B (zh) 一种全局生物网络比对方法
WO2021208993A1 (zh) 一种用于预测药物靶标的信息处理方法及装置
CN113192562B (zh) 融合多尺度模块结构信息的致病基因识别方法及系统
CN115249513A (zh) 一种基于Adaboost集成思想的神经网络拷贝数变异检测方法与系统
CN114139937A (zh) 一种室内热舒适数据生成方法、系统、设备及介质
CN114334168A (zh) 结合协同学习策略的粒子群混合优化的特征选择算法
CN114360641A (zh) 一种基于变分贝叶斯的基因调控网络结构辨识方法
CN109345274B (zh) 基于bp神经网络评分预测误差的近邻用户选择方法

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