CN111951889B - 一种rna序列中m5c位点的识别预测方法及系统 - Google Patents
一种rna序列中m5c位点的识别预测方法及系统 Download PDFInfo
- Publication number
- CN111951889B CN111951889B CN202010832292.0A CN202010832292A CN111951889B CN 111951889 B CN111951889 B CN 111951889B CN 202010832292 A CN202010832292 A CN 202010832292A CN 111951889 B CN111951889 B CN 111951889B
- Authority
- CN
- China
- Prior art keywords
- rna
- data set
- site
- prediction
- model
- 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
- 238000000034 method Methods 0.000 title claims abstract description 68
- 108091028043 Nucleic acid sequence Proteins 0.000 title claims abstract description 28
- 108091032973 (ribonucleotides)n+m Proteins 0.000 claims abstract description 70
- 239000012634 fragment Substances 0.000 claims abstract description 45
- 239000002773 nucleotide Substances 0.000 claims description 42
- 238000012549 training Methods 0.000 claims description 42
- 125000003729 nucleotide group Chemical group 0.000 claims description 39
- 239000000523 sample Substances 0.000 claims description 32
- OPTASPLRGRRNAP-UHFFFAOYSA-N cytosine Chemical compound NC=1C=CNC(=O)N=1 OPTASPLRGRRNAP-UHFFFAOYSA-N 0.000 claims description 30
- 238000012360 testing method Methods 0.000 claims description 29
- 239000013598 vector Substances 0.000 claims description 17
- 229940104302 cytosine Drugs 0.000 claims description 15
- ISAKRJDGNUQOIC-UHFFFAOYSA-N Uracil Chemical compound O=C1C=CNC(=O)N1 ISAKRJDGNUQOIC-UHFFFAOYSA-N 0.000 claims description 12
- UYTPUPDQBNUYGX-UHFFFAOYSA-N guanine Chemical compound O=C1NC(N)=NC2=C1N=CN2 UYTPUPDQBNUYGX-UHFFFAOYSA-N 0.000 claims description 12
- 238000004364 calculation method Methods 0.000 claims description 10
- 238000010276 construction Methods 0.000 claims description 8
- 230000035772 mutation Effects 0.000 claims description 8
- 230000008569 process Effects 0.000 claims description 8
- 229930024421 Adenine Natural products 0.000 claims description 6
- 229960000643 adenine Drugs 0.000 claims description 6
- GFFGJBXGBJISGV-UHFFFAOYSA-N adenyl group Chemical group N1=CN=C2N=CNC2=C1N GFFGJBXGBJISGV-UHFFFAOYSA-N 0.000 claims description 6
- 238000012512 characterization method Methods 0.000 claims description 6
- 229940035893 uracil Drugs 0.000 claims description 6
- 239000013614 RNA sample Substances 0.000 claims description 5
- 238000012545 processing Methods 0.000 claims description 5
- 238000011144 upstream manufacturing Methods 0.000 claims description 4
- 238000010353 genetic engineering Methods 0.000 claims description 3
- 238000005457 optimization Methods 0.000 claims description 3
- 230000015572 biosynthetic process Effects 0.000 claims description 2
- 239000000126 substance Substances 0.000 claims description 2
- LRSASMSXMSNRBT-UHFFFAOYSA-N 5-methylcytosine Chemical compound CC1=CNC(=O)N=C1N LRSASMSXMSNRBT-UHFFFAOYSA-N 0.000 abstract description 74
- 230000035945 sensitivity Effects 0.000 description 19
- 238000012706 support-vector machine Methods 0.000 description 17
- 238000002790 cross-validation Methods 0.000 description 16
- 230000006870 function Effects 0.000 description 8
- 230000004048 modification Effects 0.000 description 8
- 238000012986 modification Methods 0.000 description 8
- 230000002068 genetic effect Effects 0.000 description 6
- 238000007477 logistic regression Methods 0.000 description 6
- 238000002474 experimental method Methods 0.000 description 5
- 230000001124 posttranscriptional effect Effects 0.000 description 5
- 238000012163 sequencing technique Methods 0.000 description 5
- 238000004458 analytical method Methods 0.000 description 3
- 238000006243 chemical reaction Methods 0.000 description 3
- 238000002487 chromatin immunoprecipitation Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 238000010801 machine learning Methods 0.000 description 3
- 239000002777 nucleoside Substances 0.000 description 3
- 125000003835 nucleoside group Chemical group 0.000 description 3
- 108090000623 proteins and genes Proteins 0.000 description 3
- 102000004169 proteins and genes Human genes 0.000 description 3
- 101710159080 Aconitate hydratase A Proteins 0.000 description 2
- 101710159078 Aconitate hydratase B Proteins 0.000 description 2
- 108020004414 DNA Proteins 0.000 description 2
- 102000044126 RNA-Binding Proteins Human genes 0.000 description 2
- 101710105008 RNA-binding protein Proteins 0.000 description 2
- 230000008901 benefit Effects 0.000 description 2
- 230000027455 binding Effects 0.000 description 2
- 238000001369 bisulfite sequencing Methods 0.000 description 2
- 239000003795 chemical substances by application Substances 0.000 description 2
- 238000003066 decision tree Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 230000004807 localization Effects 0.000 description 2
- 239000011159 matrix material Substances 0.000 description 2
- 150000007523 nucleic acids Chemical group 0.000 description 2
- 238000000730 protein immunoprecipitation Methods 0.000 description 2
- 230000001105 regulatory effect Effects 0.000 description 2
- 238000012795 verification Methods 0.000 description 2
- LSNNMFCWUKXFEE-UHFFFAOYSA-M Bisulfite Chemical compound OS([O-])=O LSNNMFCWUKXFEE-UHFFFAOYSA-M 0.000 description 1
- 108020004705 Codon Proteins 0.000 description 1
- 102000016397 Methyltransferase Human genes 0.000 description 1
- 108060004795 Methyltransferase Proteins 0.000 description 1
- 230000004570 RNA-binding Effects 0.000 description 1
- 102000004389 Ribonucleoproteins Human genes 0.000 description 1
- 108010081734 Ribonucleoproteins Proteins 0.000 description 1
- 108091081021 Sense strand Proteins 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 230000010933 acylation Effects 0.000 description 1
- 238000005917 acylation reaction Methods 0.000 description 1
- 150000001413 amino acids Chemical class 0.000 description 1
- 230000000692 anti-sense effect Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000031018 biological processes and functions Effects 0.000 description 1
- 230000033228 biological regulation Effects 0.000 description 1
- 229910052799 carbon Inorganic materials 0.000 description 1
- 125000004432 carbon atom Chemical group C* 0.000 description 1
- 210000004027 cell Anatomy 0.000 description 1
- 239000003153 chemical reaction reagent Substances 0.000 description 1
- 230000002759 chromosomal effect Effects 0.000 description 1
- 210000000349 chromosome Anatomy 0.000 description 1
- 238000013145 classification model Methods 0.000 description 1
- 238000003776 cleavage reaction Methods 0.000 description 1
- 238000010367 cloning Methods 0.000 description 1
- 238000000205 computational method Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000012850 discrimination method Methods 0.000 description 1
- 239000002360 explosive Substances 0.000 description 1
- 230000012010 growth Effects 0.000 description 1
- 230000006058 immune tolerance Effects 0.000 description 1
- 230000002434 immunopotentiative effect Effects 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 230000007102 metabolic function Effects 0.000 description 1
- 230000004060 metabolic process Effects 0.000 description 1
- 108091070501 miRNA Proteins 0.000 description 1
- 239000002679 microRNA Substances 0.000 description 1
- 108020004707 nucleic acids Proteins 0.000 description 1
- 102000039446 nucleic acids Human genes 0.000 description 1
- 230000004481 post-translational protein modification Effects 0.000 description 1
- 238000000746 purification Methods 0.000 description 1
- 230000006798 recombination Effects 0.000 description 1
- 238000005215 recombination Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 230000007017 scission Effects 0.000 description 1
- 238000010187 selection method Methods 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 238000000638 solvent extraction Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 238000000844 transformation Methods 0.000 description 1
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/30—Detection of binding sites or motifs
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
- G06F18/214—Generating training patterns; Bootstrap methods, e.g. bagging or boosting
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N20/00—Machine learning
- G06N20/10—Machine learning using kernel methods, e.g. support vector machines [SVM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q10/00—Administration; Management
- G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
-
- 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)
- Data Mining & Analysis (AREA)
- Life Sciences & Earth Sciences (AREA)
- Business, Economics & Management (AREA)
- Bioinformatics & Cheminformatics (AREA)
- General Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Evolutionary Computation (AREA)
- Evolutionary Biology (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Bioinformatics & Computational Biology (AREA)
- Artificial Intelligence (AREA)
- Health & Medical Sciences (AREA)
- Economics (AREA)
- Software Systems (AREA)
- Human Resources & Organizations (AREA)
- Strategic Management (AREA)
- Medical Informatics (AREA)
- Marketing (AREA)
- Operations Research (AREA)
- Quality & Reliability (AREA)
- Tourism & Hospitality (AREA)
- General Business, Economics & Management (AREA)
- Entrepreneurship & Innovation (AREA)
- Game Theory and Decision Science (AREA)
- Computing Systems (AREA)
- Mathematical Physics (AREA)
- Development Economics (AREA)
- Chemical & Material Sciences (AREA)
- Analytical Chemistry (AREA)
- Biophysics (AREA)
- Genetics & Genomics (AREA)
- Molecular Biology (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Biotechnology (AREA)
- General Health & Medical Sciences (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
Abstract
本发明公开了一种RNA序列中M5C位点的识别预测方法及系统,属于M5C位点预测技术领域,包括以下步骤:S1:构建基准数据集;S2:利用特征表示RNA片段;S3:对特征进行优化选择;S4:构建预测模型;S5:利用模型进行预测。本发明构建了一个平衡的数据集Mat372,在此数据集的基础上,利用KNF、KSNPF和pseDNC三个特征对RNA片段进行编码,对特征进行优化选择后,利用SVM建立了M5C‑NSGAII预测模型,能够对RNA序列中的5‑甲基胞嘧啶(M5C)位点进行准确识别,值得被推广使用。
Description
技术领域
本发明涉及M5C位点预测技术领域,具体涉及一种RNA序列中M5C位点的识别预测方法及系统。
背景技术
自然界中,DNA的转录后加工与修饰是一种常见的现象,RNA的代谢过程和功能受特定RNA序列基序、RNA形成二级结构并组装成核糖核蛋白复合物的能力的影响。基于上述效应,记录这些分子的转录后修饰过程和程度,如剪切输出、免疫耐受等,具有重要意义。因此,人们越来越重视研究转录后修饰在RNA中的发生率和生物学相关性。5-甲基胞嘧啶(5-methylcytosine,M5C)是RNA中的一类高丰度的转录后修饰(PTCM),已经在多种生物中被发现。通过RNA甲基转移酶的作用,RNA序列中的胞嘧啶的第5位碳原子发生修饰,被转换成M5C。M5C在tRNA二级结构的稳定性、氨基酸酰化和密码子识别、应激反应调控等生物学过程中发挥着重要的作用,从而得到了广泛的研究。精确识别RNA序列中M5C位点,对于理解M5C的功能作用和机制具有重要意义和帮助。
一些实验方法,比如亚硫酸氢盐测序、M5C-rip、Aza-IP或miCLIP等已被开发用来识别M5C位点。亚硫酸氢盐测序使用亚硫酸氢盐来处理核酸序列,没有甲基化的胞嘧啶会转化为尿嘧啶,而甲基化的胞嘧啶不发生改变,可以由此识别位点。M5C-rip技术来源于RNA结合蛋白免疫沉淀(RNABindingProteinImmunoprecipitation),RNA BindingProteinImmunoprecipitation是研究细胞内RNA和蛋白结合情况的一种技术,是了解转录后调控网络动态过程的强有力工具,并且可以帮助找到miRNA的调节靶点。RIP是一种新兴的技术,运用针对目标蛋白的抗体沉淀出相应的RNA-蛋白复合物,然后经过分离纯化,结合在复合物上的RNA就可以被分析。RIP可以被视为一个类似普遍使用的染色质免疫沉淀ChIP技术的应用,但由于研究对象是RNA-蛋白复合物而不是DNA-蛋白复合物,RIP的优化条件实验不同于ChIP实验(如复合物不需要固定,RIP反应体系中的试剂和抗体必须不包含RNA酶,抗体需要通过RIP实验进行验证等等)。通过将PCR产物克隆成载体并加以测序,可以提升测序的成功概率,这种方法称为BSP-克隆测序法。miCLIP可以交联RNA-m6A抗体结合位点,当抗体结合的RNA被逆转录时,这些特异性位点会发生突变。这种独有的突变特征(例如,C-T转换或截短)的测序可以精确定位m6a,然后在此基础上发展一个可定位M5C的方法。然而,这些技术既费时又昂贵,此外,测序技术的快速发展导致RNA序列的爆炸性增长,这就需要更快、更经济的分析方法。
计算预测方法由于速度快、成本低,为RNA序列中M5C位点的识别提供了另一类方法。据悉,已有若干研究小组发展了一些计算方法用于预测RNA的M5C位点。Feng等人开发了一个模型来预测人类RNA的M5C位点。该模型建立在120个阳性样本和120个阴性样本的平衡数据集的基础上,以包含三种RNA理化性质的伪二核苷酸组合(PseDNC)作为特征对RNA序列进行编码,并采用支持向量机(SVM)作为模型的分类器。Qiu等人提出的另一个模型iRNAM5C-PseDNC,是基于一个包含475个正样本和1425个负样本的不平衡数据集建立的。该模型应用SVM作为分类器,但采用了改进的PseDNC对RNA序列进行编码。最近,Zhang等人提出了一种名为M5C-HPCR的模型。在该方法中,引入了一种启发式算法来选择部分PseDNC特征,然后采用集成方法建立了模型。M5C-HPCR分别在Feng等人和Qiu等人使用的平衡和不平衡数据集上进行了验证。然而当前此类方法都只使用PseDNC对RNA序列进行编码,其它一些序列编码方法,如K-核苷酸频率(K-nucleotidefrequency,KNF),K-间隔核苷酸对频率(K-spacednucleotidepairfrequency,KSNPF)等,都还没有被用于预测M5C位点。
虽然已经有人发展了一些计算方法用于预测RNA M5C位点,但是这些方法在使用中存在预测精度不够高等问题,因此,迫切需要开发更精确的计算方法来有效识别预测M5C位点。为解决上述问题,我们提出一种RNA序列中M5C位点的识别预测方法及系统。
发明内容
本发明所要解决的技术问题在于:如何解决现有技术中M5C位点识别方法在使用时存在的预测精度不够高等问题,提供了一种RNA序列中M5C位点的识别预测方法。
如图1所示,本发明是通过以下技术方案解决上述技术问题的,本发明包括以下步骤:
S1:构建基准数据集
建立基准数据集,将基准数据集分割成两个派生数据集分别为训练集与测试集,各数据集中的正子集包含的RNA片段拥有可以被修饰成M5C的中心胞嘧啶,负子集包含的RNA片段拥有不可以被修饰成M5C的中心胞嘧啶;
S2:利用特征表示RNA片段
利用KNF(K-核苷酸频率)、KSNPF(K-间隔核苷酸对频率)和pseDNC(伪二核苷酸组合)三个特征将基准数据集的RNA片段编码为特征向量;
S3:对特征及支持向量机的超参数进行优化选择
利用NSGAII方法对步骤S2中的特征向量及支持向量机的超参数进行优化选择;
S4:构建M5C-NSGAII预测模型
利用SVM(支持向量机)作为学习器在训练集上进行训练,建立M5C-NSGAII预测模型;
S5:利用模型对RNA序列中可能的M5C位点进行预测
利用步骤S4中的M5C-NSGAII预测模型对RNA片段样本上可能的M5C位点进行预测。
更进一步地,在所述步骤S1中,基准数据集的生成过程如下:
S11:通过GEO数据库中ID为GSE90963的记录获取高阈值的M5C位点信息;
S12:根据步骤S11中高阈值M5C位点在基因组中的位置信息,截取人类基因转录组中位于其两侧各20个碱基的RNA片段构成正样本,所有的正样本构成的数据集被命名为P1;
S13:排除GSE90963中记录的所有可能的M5C位点,根据基因转录组中其余的C位点及两侧各20个碱基的RNA片段构成负样本,将该组片段被命名为N1;
S14:去除P1和N1中的冗余序列,由于负样本的数量比较多因此采用下采样方法获得与正样本同样数量的负样本,即生成含有186例阳性样本的P2和含有186例阴性样本的N2;
S15:将各含有186例样本的P2和N2合并得到基准数据集Mat372。
更进一步地,在所述步骤S14中,利用CD-HIT去除P1和N1中的冗余序列,截断值分别为0.7。
更进一步地,在所述步骤S1中,训练集包括基准数据集的正样本和负样本中按比例各选取的149个RNA片段,样本剩余部分作为测试集。
更进一步地,在所述步骤S1中,所有数据集中RNA片段的长度均为41,将每个中心碱基处有一个潜在M5C位点的RNA样本(片段)表达如下:
Rξ(C)=N-ξN-(ξ-1)…N-1CN1…N+(ξ-1)Nξ
其中,N-ξ代表中心胞嘧啶上游的第ξ个核苷酸,而N+ξ代表中心胞嘧啶下游第ξ个核苷酸;
将上式简化如下:
R20(C)=N1N2…N20CN22…N40N41
其中,Ni(i=1,2,…20,21…41)表示RNA片段的第i位的核苷酸,为RNA中4个核苷酸碱基中的任意一个,即:
Ni∈{A(adenine),C(cytosine),G(guanine),U(uracil)}
其中,A表示腺嘌呤;C表示胞嘧啶;G表示鸟嘌呤;U表示尿嘧啶。
更进一步地,在所述步骤S2中,对于给定的K值,KNF表示RNA序列中每个K-mer核苷酸组分的发生频率,因此将RNA样本转化为一个4K维的特征向量,KNF的计算公式如下:
其中,n1n2…nK表示一个K-mer核苷酸组分,nK可以是RNA中4个核苷酸碱基中的任意一个,即:
nK∈{A(adenine),C(cytosine),G(guanine),U(uracil),N(n1n2…nK)表示样本序列中出现n1n2…nK的次数,L表示核苷酸段的长度(L=41)。
更进一步地,在所述步骤S2中,采用N1x{K}N2(N1,N2和x∈{A,C,G,U})表示一个K间隔的核苷酸对,KSNPF的计算公式如下:
其中,N(N1x{K}N2)表示样本序列中出现N1x{K}N2的次数,L表示核苷酸段的长度(L=41)。
更进一步地,在所述步骤S2中,pseDNC利用自由能、亲水性、堆积能三种物理化学性质来对RNA片段进行编码。
更进一步地,在所述步骤S3中,利用NSGAII方法对步骤S2中的特征向量进行优化选择的过程如下:
S31:总体初始化
随机产生固定数量规模的初始种群,对种群初始进行初始化;
S32:非支配排序
对初始化后种群总体按照非支配进行排序;
S33:计算拥挤距离
完成非支配排序后,对种群总体中的所有个体分配一个拥挤距离值;
S34:个体选择
在对各个体分配拥挤距离后,使用拥挤比较算子对种群中的个体进行选择;
S35:遗传操作
利用模拟二进制交叉、多项式突变操作得到第一代子代种群;
S36:形成新种群
从第二代开始,将父代种群与子代种群合并,重复步骤S33~S35,产生新的子代种群,直到满足预设条件。(种群中的个体是由选择出的特征及选择出的SVM的超参数g和c来构成的)
本发明还提供了一种RNA序列中M5C位点的识别预测系统,包括:
数据集构建模块,用于建立基准数据集,从基准数据集派生出两个子集分别为训练集与测试集;
特征表示模块,用于利用KNF、KSNPF和pseDNC三个特征将基准数据集的RNA片段编码为特征向量;
特征选择模块,用于利用NSGAII方法对步骤S2中的特征向量进行优化选择;
模型构建模块,用于利用SVM作为学习方法对训练集进行训练,建立M5C-NSGAII预测模型;
预测模块,用于利用M5C-NSGAII预测模型对RNA片段样本上M5C位点进行预测;
控制处理模块,用于向其他模块发出指令,完成相关动作;
所述数据集构建模块、特征表示模块、特征选择模块、模型构建模块、预测模块与控制处理模块电连接。
本发明相比现有技术具有以下优点:该RNA序列中M5C位点的识别预测方法及系统,构建了一个平衡的数据集Mat372,在此数据集的基础上,利用KNF、KSNPF和pseDNC三个特征对RNA片段进行编码,对特征进行优化选择后,利用SVM建立了M5C-NSGAII预测模型,能够对RNA序列中的5-甲基胞嘧啶(M5C)位点进行准确识别,值得被推广使用。
附图说明
图1是本发明识别方法的总体流程示意图;
图2是本发明实施例中生成基准数据集Mat372的流程图;
图3是本发明实施例中利用交叉验证结果绘制的ROC图;
图4是本发明实施例中M5C-Bayes模型在训练集上的ROC图;
图5是本发明实施例中M5C-knn模型在训练集上的ROC图;
图6是本发明实施例中利用独立测试集Test74预测结果绘制的ROC图。
具体实施方式
下面对本发明的实施例作详细说明,本实施例在以本发明技术方案为前提下进行实施,给出了详细的实施方式和具体的操作过程,但本发明的保护范围不限于下述的实施例。
如图2~6所示,本实施例提供一种技术方案:一种RNA序列中M5C位点的识别预测方法,主要技术方案包括以下内容:
1、选取基准数据集
在本实施例中,使用了一个基准数据集(即Mat372)来训练和验证我们的模型,将该数据集划分为两个派生数据集Train298和Test74。这些数据集皆由一个正子集和一个负子集组成。正子集包含的为RNA序列拥有可以被修饰成M5C的中心胞嘧啶,而负子集包含的为RNA序列拥有不可以被修饰成M5C的中心胞嘧啶。
如图2所示,Mat372的生成流程如下:
(1)从“GSE90963_Table_S1-M5C_candidate_sites”中获取M5C位点的高阈值信息。
下载自GEO(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc= GSE90963)
(2)根据(1)中的“高阈值”M5C位点在基因组中的位置信息,具体信息为该位点所处的是DNA正义链还是反义链,该位点所在的染色体是哪一个以及该位点在染色体序列中的具体位置。根据这些信息我们就可以在转录组中找到相应的位点及核苷,然后在其上下游各取20个核苷,从而获得了包含41个核苷的RNA片段,所有这些片段构成的集合被命名为P(正数据集)1。
(3)排除GSE90963_Table_S1-M5C_candidate_sites.xlsx文件中记录的可能存在的M5C位点,从而生成所有其他的以C为中心的41-tupleRNA片段,这些片段的集合被命名为N(负数据集)1。
(4)使用CD-HIT(CD-HIT是一款广泛用于对蛋白或者核酸序列进行聚类和比较的软件。)去除P1和N1中的冗余序列,截断值分别为0.7。对于P1即生成了含有186例阳性样本的P2,对于N1生成的序列数量远多于P2,因此使用下采样方法从中随机选取186例阴性样本作为N2。
(5)将各含有186例样本的P2和N2合并得到基准数据集Mat372。
Train298则是通过在Mat372的正样本和负样本中按比例各选取149个RNA片段组成,其余部分作为测试集Test74。
所有数据集中RNA片段的长度均为41,可以将中心碱基处存在一个潜在M5C位点的RNA样本(片段)表达如下:
Rξ(C)=N-ξN-(ξ-1)…N-1CN1…N+(ξ-1)Nξ (1)
其中,N-ξ代表中心胞嘧啶上游的第ξ个核苷酸,而N+ξ代表中心胞嘧啶下游第ξ个核苷酸。
为了进一步简化描述,将式(1)改写为:
R20(C)=N1N2…N20CN22…N40N41 (2)
其中,Ni(i=1,2,…20,21…41)表示RNA片段的第i位的核苷酸,可以是RNA中4个核苷酸碱基中的任意一个,即:
Ni∈{A(adenine),C(cytosine),G(guanine),U(uracil)} (3)
如下表1所示,表1显示了所用数据集的信息:
表1、三个数据集的信息
数据集 | 长度(bp) | 正子集 | 负子集 | 总计 |
Mat372 | 41 | 186 | 186 | 372 |
Train298 | 41 | 149 | 149 | 298 |
Test74 | 41 | 37 | 37 | 74 |
2、RNA片段的特征表示
将RNA片段编码为具有高识别率信息的特征向量,对建立预测M5C位点的机器学习模型具有重要意义。KNF(K-核苷酸频率)、KSNPF(K-间隔核苷酸对频率)和pseDNC(伪二核苷酸组合)常被用来预测RNA的翻译后修饰位点。在本实施例中,探讨了这三个特征预测M5C位置的能力,每个特征的详细信息如下:
2.1、KNF特征
KNF是表示核苷酸序列特征的经典方法。对于给定的K值,KNF表示RNA序列中每个K-mer核苷酸组分的发生频率。因此,RNA样本可以转化为一个4K维的特征向量。
采用下式计算KNF:
其中,n1n2…nK表示一个K-mer核苷酸组分,nK可以是RNA中4个核苷酸碱基中的任意一个,即
nK∈{A(adenine),C(cytosine),G(guanine),U(uracil) (5)。
其中,N(n1n2…nK)表示样本序列中出现n1n2…nK的次数,L表示核苷酸段的长度(L=41)。例如,当K=2时,有16种二核苷酸,RNA片段可以编码为:
R(2NF)=[fAAfACfAGfAUfCAfCCfCGfCUfGAfGCfGGfGUfUAfUCfUGfUU](42=16d) (6)
显然,随着K值的增加,特征向量的维数以指数式速度增长,为了避免“维数灾难”,本实施例中使用的最大K值设置为4。
2.2、KSNPF特征
KSNPF也被用来编码RNA序列。K-间隔的核苷酸对是由K个任意核苷酸分开的核苷酸对,例如UxxxG是一个3间隔的核苷酸对,在核苷酸U和G之间有3个任意核苷酸。为了说明这一点,采用N1x{K}N2(N1,N2和x∈{A,C,G,U})来表示一个K间隔的核苷酸对。显然,对于由N1和N2组成的核苷酸对,这里有(4×4=)16种可能的组合。
换句话说,对于一个固定的K值,有16种K-间隔的核苷酸对。例如,当K=3时,N1x{K}N2可能是AxxxA、AxxxC、……或UxxxU。
和KNF一样,可以用下面的公式来计算KSNPF:
N(N1x{K}N2)表示样本序列中出现N1x{K}N2的次数,L表示核苷酸段的长度(L=41)。在本实施例中,尝试了不同的K值(即1、2、3、4、5),例如当K=1时,RNA片段可以编码为:
R(1SNPF)=[f(AxA)f(AxC)f(AxG)……f(UxC)f(UxC)f(UxU)](16d) (7)
2.3、pseDNC特征
pseDNC是一个可以结合RNA片段的局部和全局序列模式信息的特征。经过一系列的自协方差和互协方差变换后,由物理化学矩阵导出pseDNC的分量。据悉,至少有23种物理化学性质可以被pseDNC用来编码RNA片段。在本实施例中,从这23种理化性质中选择了3种,即、自由能、亲水性和堆积能。表2列出了详细信息。
表2、RNA中二核苷酸的理化性质表
二核苷酸 | 自由能 | 亲水性 | 堆积能 |
GG | -3.260 | 0.170 | -11.100 |
GA | -2.350 | 0.100 | -14.200 |
GC | -3.420 | 0.260 | -16.900 |
GU | -2.240 | 0.270 | -13.800 |
AG | -2.080 | 0.080 | -14.000 |
AA | -0.930 | 0.040 | -13.700 |
AC | -2.240 | 0.140 | -13.800 |
AU | -1.100 | 0.140 | -15.400 |
CG | -2.360 | 0.350 | -15.600 |
CA | -2.110 | 0.210 | -14.400 |
CC | -3.260 | 0.490 | -11.100 |
CU | -2.080 | 0.520 | -14.000 |
UG | -2.110 | 0.340 | -14.400 |
UA | -1.330 | 0.210 | -16.000 |
UC | -2.350 | 0.480 | -14.200 |
UU | -0.930 | 0.440 | -13.700 |
3、NSGAII作为特征选择方法
3.1、NSGAII提出背景
首先第一代(NSGA)为一种流行的基于非控制的多目标优化遗传算法。它是一种非常有效的算法,但是由于计算复杂性、缺乏精英性以及为共享参数选择最优参数值而受到普遍的批评。其后Srinivas和Deb提出了一个改进版本,NSGAII,它具有更好的排序算法,结合了精英策略,不需要先验地选择共享参数。
3.2、NSGAII主要思想
首先,随机产生规模为N的初始种群,非支配排序后通过遗传算法的选择、交叉、变异三个基本操作得到第一代子代种群;其次,从第二代开始,将父代种群与子代种群合并,进行快速非支配排序,同时对每个非支配层中的个体进行拥挤度计算,根据非支配关系以及个体的拥挤度选取合适的个体组成新的父代种群;最后,通过遗传算法的基本操作产生新的子代种群。依此类推,直到满足预先设定的条件。
3.3、NSGAII
3.3.1、总体初始化
总体是根据问题范围和约束进行初始化的。
3.3.2、非支配排序
初始化后的总体是按照非支配排序的。
快速排序算法具体步骤如下:
(1)针对总体P中的每个个体p执行以下操作:
①初始化这个集合包含所有由p支配的个体;
②初始化np=0。这是支配p的个体数;
③对于P中的每个个体q,如果q被p支配,将q放入集合Sp中,也就是说Sp=Sp∪{q},反之如果q支配p,则增加p的支配计数值,也就是说np=np+1。
④如果np=0,意味着没有个体支配p,则p属于第一个前沿集;将个体p的等级设为1,即prank=1。通过将p加入第一个前沿集的第一个元素上更新第一个前沿集即为F1=F1∪{p}。
(2)对总体P中的所有个体执行上述操作;
(3)将前沿计数值初始化为1,i=1。
(4)当第i个前沿不为空时执行以下操作
①这个集合用于存储第(i+1)个前沿的个体
②对于前沿Fi的每个个体p:
对于每个Sp中的个体q(Sp是p所支配的个体的集合):
nq=nq-1,减少个体q的支配计数值;
如果nq=0则在随后的前沿中没有个体支配q。因此设qrank=i+1。加入q来更新集合Q(Q=Q∪q)。
③将前沿计数值增加1。
④现在集合Q是下一个前沿且Fi=Q。
3.3.3、计算拥挤距离
一旦非支配排序完成,即赋予了拥挤距离。由于个体是根据等级和拥挤距离选择的,因此种群中的所有个体都被赋予一个拥挤距离值。
拥挤距离的计算如下:
对于每个前沿Fi,n是个体的数量:
(1)对于所有的个体,初始化距离为0即Fi(dj)=0,其中j意味着前沿Fi中第j个个体。
(2)对于每个目标函数m:
①基于目标m为前沿Fi中的个体排序,即为I=sort(Fi,m);
②对于Fi中的每个个体将无限距离赋给边界值,即I(d1)=∞和I(dn)=∞;
③k取值2到(n-1):
其中,I(k).m是I中第k个个体的第m个目标函数的值。
拥挤距离的基本思想是基于m维超空间中每个个体的m目标,找出前沿中每个个体之间的欧氏距离。边界上的个体总是被选中,因为它们被分配了无限距离。
3.3.4、个体选择
一旦根据非支配性对个体进行排序,在分配了拥挤距离的情况下,使用拥挤比较算子(<n)进行选择。比较的基础如下:
(1)非支配等级为prank,前沿Fi中的个体有各自的等级为prank=i;
(2)拥挤距离Fi(dj):对于p<nq,如果prank<qrank或是p和q属于相同的前沿Fi,则Fi(dp)>Fi(dq),也就是说拥挤距离应该更大。
个体选择通过一个具有拥挤比较运算符的二元竞争选择进行。
3.3.5、遗传操作
GA(目标优化遗传算法)的实码使用模拟二进制交叉(SBX)运算符,用于交叉和多项式突变。
3.3.5.1、模拟二进制交叉
模拟二进制交叉模拟了自然界中观测到的二元交叉,如下所示:
其中,ci,k是第k分支的第i个子代,pi,k是被选择的父代且βk(≥0)是从具有以下密度函数的随机数中生成的样本,其分布如下所示:
此分布可以从(0,1)之间的均匀采样随机数μ获得。ηc是用来进行交叉的分布指数,也就是:
其中,μ是从(0,1)之间均匀采样得到的一个随机数。
3.3.5.2、多项式突变
利用如下公式进行多项式突变:
其中,ck是子代,pk是父代,是父代部分的上界,/>是下界,δk是由如下多项式分布计算出的微小变化:
rk是介于(0,1)之间的均匀采样随机数,ηm是变异分布指数。
3.3.6、重组和选择
后代群体和当代群体相结合,进行选择以确定下一代的个体。由于所有前代和当代最好的个体都加入了种群,因此精英主义得到了保证。种群现在是根据非支配性进行排序的。新一代被各个前沿填充,直到种群规模超过当前种群规模。如果将前沿Fj中的所有个体累计,种群大小超过N,则前沿Fj中的个体依据拥挤距离从大到小进行选择直到种群大小达到N。重复这个过程来产生下一代。
4、利用SVM(支持向量机)作为预测引擎进行预测
SVM是目前最常用的适用于小样本数据集和高维度特征集的一种分类方法。对于低维空间内的线性不可分问题,SVM通过使用某种非线性变换将输入的特征值投射到高维空间,变成高维空间中的线性可分问题,进而在高维空间进行求解。SVM方法最重要的一环就是引进了核函数,用以解决高维空间的内积问题。核函数包括很多种,其中最常用的是高斯核函数,它可以将数据映射到无穷维,也称作径向基函数。在本实施例中使用SVM作为分类器,并采用MATLAB实现算法。
5、性能评估
本实施例中运用交叉验证来评估模型的泛化性能。在统计学习领域,经常使用以下三种验证方法来评估预测器的性能:二折交叉验证测试、K折交叉验证测试和刀切(jackknife)测试。二折交叉验证把数据集分成两个大小相同的子集,选择一个子集作为训练集,另一个子集作为测试集,然后交换两个集合的职能。总误差为对两次运行求和。K折交叉验证(k-foldcross-validation)是对二折的推广,首先将所有数据分割成K个子样本,不重复的选取其中一个子样本作为测试集,其他K-1个样本用来训练。总共重复K次,平均K次的结果或者使用其它指标,最终得到一个单一估测。这个方法的优点是确保每个子样本都参与训练并进行测试,以降低泛化误差。其中,10折交叉验证是最常用的。在jackknife测试中,基准数据集中的所有样本将被逐个挑选出来,并通过对剩余样本进行训练的预测器进行测试。因此,jackknife测试被认为是最不随意的验证方法,对于给定的基准数据集总是可以产生唯一的结果。由于交叉验证中使用的训练集仅比原始数据集少一个样本,因此其结果与真实模型最为接近。因此,使用jackknife测试来评估我们的模型。对于预测模型,其泛化性能非常重要的。在这里,我们使用独立的测试集来进一步证明我们的模型具有良好的泛化性能。
在实验中,我们采用了敏感性(Sen)、特异性(Spe)、精确率(Pre)、准确率(Acc)和马修斯相关系数(MCC)五个常用的评价指标来检验本发明方法的性能。
它们的定义如下:
其中,TP、TN、FP和FN分别代表真阳性、真阴性、假阳性和假阴性的预测结果计数。
这个公式也可以写成:
其中,N+是指阳性样本或真正的5-甲基胞嘧啶修饰位点的总数,是错误地预测为非5-甲基胞嘧啶位点的真实5-甲基胞嘧啶修饰位点的数目,N-是阴性样本的总数或非5-甲基胞嘧啶位点的数目,而/>是错误地预测为5-甲基胞嘧啶位点的非5-甲基胞嘧啶位点的数目
根据以上公式,可以清楚地看到:当意味着没有真正的5-甲基胞嘧啶位点被错误地预测为非5-甲基胞嘧啶位点时,具有Sen=1的灵敏度。当/>意味着所有真正的5-甲基胞嘧啶位点被错误地预测为非5-甲基胞嘧啶位点时,具有Sen=0的灵敏度。同样,当/>意味着没有非5-甲基胞嘧啶位点被错误地预测为5-甲基胞嘧啶位点时,具有特异性Spe=1;而/>意味着错误地预测了所有非5-甲基胞嘧啶位点为真正的5-甲基胞嘧啶位点,有特异性Spe=0。/>的意义为真正的5-甲基胞嘧啶位点被正确地预测为5-甲基胞嘧啶位点的数量与非5-甲基胞嘧啶位点被错误地预测为5-甲基胞嘧啶位点的数量相同,此时具有Pre=0.5的精确率,当/>意味着没有非5-甲基胞嘧啶位点被错误地预测为5-甲基胞嘧啶位点时,具有Pre=1的精确率。当/>时,意味着没有错误的预测正数据集中的所有5-甲基胞嘧啶位点,也没有错误预测负数据集中的非5-甲基胞嘧啶位点,总体准确度Acc=1,MCC=1;当/>和/>时,意味着正数据集中所有真实的5-甲基胞嘧啶位点和负数据集中所有的非5-甲基胞嘧啶位点都被错误地预测,整体精度Acc=0且MCC=-1;而当/>和/>时,Acc=0.5和MCC=0,意味着与随机猜测一致。因此,第二个方程使敏感度,特异性,精确率,整体准确性和稳定性的含义更加直观和易于理解。
此外,为了进一步评价模型,我们采用了受试者工作特性曲线(ROC)和曲线下面积(AUC)。ROC曲线(receiveroperatingcharacteristiccurve),又称为感受性曲线(sensitivitycurve)。之所以叫这个名字,是因为曲线上的每一个点都反映着相同的感受性,它们都是对相同的信号刺激的反应,只不过是在两个不同的判定标准下所得的结果而已。ROC曲线是根据一系列不同的二分类方式(分界值或决定阈),同时以假阳性率(1-特异度)为横坐标,真阳性率(灵敏度)为纵坐标绘制而成的。ROC曲线是反映连续变量灵敏度和特异性的综合指标,该曲线通过构图法揭示了灵敏度和特异性之间的关系,设定多个连续变量的临界值,计算出一系列灵敏度和特异性的值,再以特异性为横坐标,灵敏度为纵坐标绘制成曲线。在ROC曲线上,最靠近坐标图左上方的点为特异性和灵敏度均较高的阈值点。ROC曲线下方的面积(AreaundertheCurveofROC(AUCROC)),其意义是:
(1)因为是在1x1的方格里求面积,AUC必在0到1之间;
(2)假设阈值以上是阳性,以下是阴性;
(3)随机选择一个阳性样本和一个阴性样本,分类器正确判断阳性样本的值高于阴性样本的值的概率为AUC。
ROC(receiveroperatingcharacter)曲线是能直观展现预测器性能的一个重要指标。研究表明,曲线下侧包含的面积(areaundercurve,AUC)越接近1,模型对样本的预测性能就越好,根据AUC判断分类器(预测模型)优劣的标准如下:
(1)AUC=1时,为完美分类器。使用该预测模型时,至少存在一个能得出完美预测的阈值。但在大多数情况下,没有完美分类器。
(2)0.5<AUC<1时,优于随机猜测。如果阈值设置正确,则该分类器(预测模型)可以具有预测价值。
(3)AUC=0.5时,跟随机猜测一样(例:扔硬币),模型没有预测价值。
(4)AUC<0.5时,比随机猜测还差;模型不可用。
简单说:AUC值的大小决定了分类器正确率的高低。
6、特征优化选择
在本实施例中选择1SNPF+pseDNC+4NF为起始特征组合,以这种特征组合得出训练集的特征和测试集的特征。
使用NSGAII遗传算法对特征进行优化选择,参数pop=50,gen=50,得出295×50的矩阵,其中前289列(数据集特征的列数为289)为是否选择相应特征的随机概率值,其中随机概率值大于0.5的对应特征将会被选择,第290、291列为进行模型建立时重要的参数c和g的指数(以2为底),第292列为遗传算法选择特征占所有特征的百分比,第293列为1/F1-score的值,也是评价模型性能的重要数值,第294、295列分别为rank值和crowdingdistance(拥挤距离)值。
表3、NSGAII特征选择运行结果(部分)
/>
由表3可知1/F1-score的最小数值为1.3595(此行F1-score值最大),即这一行将作为特征选择的标准,用这一行前289列来对特征进行优化,选出153列特征。由参数c和g(以2为底)的指数(四舍五入)可得参数c和g分别为0.25和32,这将用于后续的建模。
7、交叉验证及结果
在本实施例中,将数据集分成了训练集Train298和测试集Test74,我们先在Train298上做训练和分析,而Test74必须在模型训练完成后才被用来评估模型误差。交叉验证是在训练集上的进一步的分析,它可以从有限的数据中获得尽可能多的有效信息,从而可以从多角度学习样本,避免受局部极值的影响。交叉验证是评价预测模型拟合性能的有效方法。其对模型的各方面性能都能提供有参考意义的数值,以便于我们进行准确的分析。
表4、交叉验证结果
Index | Sen(%) | Spe(%) | Pre(%) | Acc(%) | Mcc |
Value | 78.52 | 57.05 | 64.64 | 67.79 | 0.3642 |
预测结果中敏感度Sen=0.7852,表示此模型能够正确判断正样本的概率达到0.7852,大于0.5,能够正确预测出绝大部分的正样本;同样特异性Spe=0.5705,表示此模型能够正确预测出较大部分的负样本;精确率Pre=0.6464,表示模型预测正样本的准确率较高。准确度Acc=0.6779,表示模型判断真例的正确性较高。Mcc=0.3642,处于0-1之间,表示模型的稳定性能较好。从整体上看,该模型在训练集上的预测结果较为准确,效果较好。
由图3可以看出AUC>0.5,matlab得出的AUC=0.7018,说明本方法的模型在训练集上的分类结果正确率较高,交叉验证结果较为理想。
7.1、基于训练集的不同模型比较
为了更好的体现本发明方法建立的M5C-NSGAII模型性能,将分别运用朴素贝叶斯(bayes)算法、knn算法和逻辑回归(logistics)算法建模,并与M5C-NSGAII模型进行比较。
(1)朴素贝叶斯算法
朴素贝叶斯法为基于贝叶斯定理与特征条件独立假设的分类方法。早期机器学习方法中应用最为广泛的两种分类模型是决策树模型(DecisionTreeModel)和朴素贝叶斯模型(BayesianModel,NBM)。
由于和决策树模型相比,朴素贝叶斯分类器(BayesClassifier,或NBC)有古典数学的理论基础,以及坚实的数学基础,其分类效率也较为稳定,因此本实施例要建立的便是朴素贝叶斯模型。朴素贝叶斯模型不需要估计太多参数,缺失数据对其没有过多影响,算法相对简单。按理说,朴素贝叶斯算法是最优的,但在实际操作中,该方法容易受到其限定假设的影响,即模型的属性相互独立。
我们将运用此种方法建立的模型命名为M5C-Bayes,在建立模型时将其在训练集上的表现结果与M5C-NSGAII在训练集上的表现结果(交叉验证结果)进行对比,如表5。
表5、M5C-NSGAII、M5C-Bayes在训练集上表现结果
模型 | Sen(%) | Spe(%) | Pre(%) | Acc(%) | Mcc |
M5C-NSGAII | 78.52 | 57.05 | 64.64 | 67.79 | 0.3642 |
M5C-Bayes | 69.80 | 63.76 | 65.82 | 66.78 | 0.3362 |
由表5可以看出,M5C-Bayes的敏感度Sen、特异性Spe、精确率Pre、准确度Acc和稳定性Mcc指标数值整体与M5C-NSGAII的相近,但仍有一定差距,M5C-NSGAII的性能高于M5C-Bayes的性能。
我们绘制出了M5C-Bayes模型在训练集上的ROC图,由图4可看出AUC明显高于0.5,而运算结果AUC=0.7251也验证了这一看法。与M5C-NSGAI I的AUC=0.7018相比较,可知M5C-NSGAII模型在训练集上分类的表现与M5C-Bayes相近。
综合所有参数考虑,M5C-NSGAII模型在训练集上的性能略好于M5C-Bayes。
(2)knn算法
knn算法的核心思想是,如果特征空间中一个样本的k个相邻样本中多数都属于某个类,则该样本也属于该类,并且具有该类上样本的特征。在分类决策中,该方法只依据最近的一个或多个样本的类别来确定要分类的样本类别。knn方法的类别决策过程中只涉及到非常少的邻近样本。由于knn方法主要依赖于周围有限的相邻样本,而不是依靠类的判别方法来确定类别,所以knn方法比其他方法更适合于用更多重叠或重叠类域来划分样本集。因此我们利用knn算法建立一个模型M5C-knn,得出模型在训练集上敏感度、特异性、准确度和稳定性的值,将其在训练集上的表现结果与M5C-NSGAII在训练集上的表现结果(交叉验证结果)进行对比。
表6、M5C-NSGAII、M5C-knn在训练集上表现结果
模型 | Sen(%) | Spe(%) | Pre(%) | Acc(%) | Mcc |
M5C-NSGAII | 78.52 | 57.05 | 64.64 | 67.79 | 0.3642 |
M5C-knn | 73.15 | 51.68 | 60.22 | 62.42 | 0.2543 |
由表6可以看出,M5C-NSGAII的灵敏度Sen、特异性Spe、精确率Pre、准确度Acc和稳定性Mcc指标数值都比M5C-knn的高,综合来看M5C-NSGAII比M5C-knn的性能要好。
我们绘制出了M5C-knn模型在训练集上的ROC图,由图5可看出AUC高于0.5,实际运算结果AUC=0.5941<0.7018(M5C-NSGAII的AUC)由此可判断M5C-NSGAII在训练集上的表现比M5C-knn更好。
(3)逻辑回归算法
逻辑回归(LogisticsRegression)是机器学习中的一种二分类模型,该算法简单、高效,在实际中有着非常广泛的应用。Logistics回归通过使用其固有的logistics函数估计概率,来衡量因变量(我们想要预测的标签)与一个或多个自变量(特征)之间的关系。Logistics回归通过线性边界将您的输入划分为两个区域,每个区域对应一个类别。因此,在建立逻辑回归模型时,数据应当是线性可分的,我们用逻辑回归算法构建一个模型M5C-Logistics,该模型在训练集上的结果如表7所示。
表7、M5C-NSGAII、M5C-knn在训练集上表现结果
模型 | Sen(%) | Spe(%) | Pre(%) | Acc(%) | Mcc |
M5C-NSGAII | 78.52 | 57.05 | 64.64 | 67.79 | 0.3642 |
M5C-Logistics | 53.02 | 51.68 | 52.32 | 52.35 | 0.0470 |
由表7可明显看出M5C-Logistics模型的结果较差,与M5C-NSGAII的差距较大,因此无需绘制ROC图便可判断M5C-NSGAII模型的性能优于M5C-Logistics模型。
综合(1)、(2)、(3)三点所述,M5C-NSGAII的性能最优。
7.2、独立测试集上预测的结果
在模型符合预期的前提下,我们用M5C-NSGAII对测试集上的数据进行预测,表8表明了我们的模型M5C-NSGAII在独立测试集Test74上的性能。
表8、基于独立测试集Test74预测的结果
Sen(%) | Spe(%) | Pre(%) | Acc(%) | Mcc |
75.68 | 67.57 | 70.00 | 71.62 | 0.4339 |
与表4中在训练集上得出的预测结果相比,预测的敏感度Sen=0.7568<0.7852,表明在训练集上正确预测出正样本的概率较测试集的高,特异性Spe=0.6757>0.5705,表明在测试集上正确预测出负样本的概率较高,精确率Pre=0.7000>0.6464,准确度Acc=0.7162>0.6779,表明在测试集上判断真例的正确性较高,MCC=0.4339>0.3642,表明模型的稳定能较好。
除此之外,我们还利用matlab软件根据预测的结果绘制出了模型的ROC曲线(如图6所示)。
此处AUC=0.7531>0.7,说明此模型对于预测M5C修饰位点性能良好。
综合考虑,本发明的预测器在测试集上的表现相对较为理想,达到了预期,模型能在测试集上实现最优性能。
综上所述,本实施例的RNA序列中M5C位点的识别预测方法,构建了一个平衡的数据集Mat372,在此数据集的基础上,利用KNF、KSNPF和pseDNC三个特征对RNA片段进行编码,对特征进行优化选择后,利用SVM建立了M5C-NSGAII预测模型,能够对RNA序列中的5-甲基胞嘧啶(M5C)位点进行准确识别,值得被推广使用。
尽管上面已经示出和描述了本发明的实施例,可以理解的是,上述实施例是示例性的,不能理解为对本发明的限制,本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。
Claims (6)
1.一种RNA序列中M5C位点的识别预测方法,其特征在于,包括以下步骤:
S1:构建基准数据集
建立基准数据集,将基准数据集划分为两个派生数据集,分别为训练集与测试集;
在所述步骤S1中,所有数据集中RNA片段的长度均为41个碱基,将每个中心碱基处有一个潜在M5C位点的RNA片段表达如下:
Rξ(C)=N-ξN-(ξ-1)…N-1CN1…N+(ξ-1)Nξ
其中,N-ξ代表中心胞嘧啶上游的第ξ个核苷酸,而N+ξ代表中心胞嘧啶下游第ξ个核苷酸;
将上式简化如下:
R20(C)=N1N2…N20CN22…N40N41
其中,Ni(i=1,2,…20,21…41)表示RNA片段的第i位的核苷酸,为RNA中4个核苷酸碱基中的任意一个,即:
Ni∈{A,C,G,U}
其中,A表示腺嘌呤;C表示胞嘧啶;G表示鸟嘌呤;U表示尿嘧啶;
S2:利用特征表示RNA片段
利用KNF、KSNPF和pseDNC三个特征将基准数据集的RNA片段编码为特征向量;
在所述步骤S2中,对于给定的K值,KNF表示RNA序列中每个K-mer核苷酸组分的发生频率,因此将RNA样本转化为一个4K维的特征向量,KNF的计算公式如下:
其中,n1n2…nK表示一个K-mer核苷酸组分,nK可以是RNA中4个核苷酸碱基中的任意一个,即:
nK∈{A,C,G,U},N(n1n2…nK)表示样本序列中出现n1n2…nK的次数,L表示核苷酸段的长度;
在所述步骤S2中,采用N1x{K}N2(N1,N2,和x∈{A,C,G,U})表示一个K间隔的核苷酸对,KSNPF的计算公式如下:
其中,N(N1x{K}N2)表示样本序列中出现N1x{K}N2的次数,L表示核苷酸段的长度;
在所述步骤S2中,pseDNC利用自由能、亲水性、堆积能三种物理化学性质对RNA片段进行编码;
S3:对特征进行优化选择
利用NSGAII方法对步骤S2中的特征向量进行优化选择;
S4:构建预测模型
利用SVM作为学习器基于特征向量进行训练,建立预测模型;
S5:利用模型进行预测
利用步骤S4中的预测模型对基准数据集中的RNA片段样本上M5C位点进行预测。
2.根据权利要求1所述的一种RNA序列中M5C位点的识别预测方法,其特征在于:在所述步骤S1中,基准数据集的生成过程如下:
S11:通过GEO数据库中ID为GSE90963的记录获取高阈值的M5C位点信息;
S12:根据步骤S11中高阈值M5C位点在基因组中的位置信息,截取人类基因转录组中位于其两侧各20个碱基的RNA片段构成正样本,所有的正样本构成的数据集被命名为P1;
S13:排除GSE90963中记录的所有可能的M5C位点,根据基因转录组中其余的C位点及两侧各20个碱基的RNA片段构成负样本,将该组片段被命名为N1;
S14:使用CD-HIT去除P1中的冗余序列,即生成含有186例阳性样本的P2;使用CD-HIT去除N1中的冗余序列并从中随机选择186例阴性样本得到N2;
S15:将各含有186例样本的P2和N2合并得到基准数据集。
3.根据权利要求2所述的一种RNA序列中M5C位点的识别预测方法,其特征在于:在所述步骤S14中,利用CD-HIT去除P1和N1中的冗余序列,截断值分别为0.7。
4.根据权利要求1所述的一种RNA序列中M5C位点的识别预测方法,其特征在于:在所述步骤S1中,训练集包括基准数据集的正样本和负样本中按比例各选取的149个RNA片段,样本剩余部分作为测试集。
5.根据权利要求1所述的一种RNA序列中M5C位点的识别预测方法,其特征在于:在所述步骤S3中,利用NSGAII方法对步骤S2中的特征向量进行优化选择的过程如下:
S31:总体初始化
随机产生固定数量规模的初始种群,对种群初始进行初始化;
S32:非支配排序
对初始化后种群总体按照非支配进行排序;
S33:计算拥挤距离
完成非支配排序后,对种群总体中的所有个体分配一个拥挤距离值;
S34:个体选择
在对各个体分配拥挤距离后,使用拥挤比较算子对种群中的个体进行选择;
S35:遗传操作
利用模拟二进制交叉、多项式突变操作得到第一代子代种群;
S36:形成新种群
从第二代开始,将父代种群与子代种群合并,重复步骤S33~S35,产生新的子代种群,直到满足预设条件。
6.一种RNA序列中M5C位点的识别预测系统,其特征在于,根据如权利要求1~5任一项所述识别预测方法对RNA序列中M5C位点进行识别预测,包括:
数据集构建模块,用于建立基准数据集,从基准数据集派生出两个数据集分别为训练集与测试集;
特征表示模块,用于利用KNF、KSNPF和pseDNC三个特征将基准数据集的RNA片段编码为特征向量;
特征选择模块,用于利用NSGAII方法对步骤S2中的特征向量进行优化选择;
模型构建模块,用于利用SVM作为学习方法对训练集进行训练,建立预测模型;
预测模块,用于利用预测模型对新的RNA片段样本上M5C位点进行预测;
控制处理模块,用于向其他模块发出指令,完成相关动作;
所述数据集构建模块、特征表示模块、特征选择模块、模型构建模块、预测模块与控制处理模块电连接。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010832292.0A CN111951889B (zh) | 2020-08-18 | 2020-08-18 | 一种rna序列中m5c位点的识别预测方法及系统 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010832292.0A CN111951889B (zh) | 2020-08-18 | 2020-08-18 | 一种rna序列中m5c位点的识别预测方法及系统 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN111951889A CN111951889A (zh) | 2020-11-17 |
CN111951889B true CN111951889B (zh) | 2023-12-22 |
Family
ID=73343095
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010832292.0A Active CN111951889B (zh) | 2020-08-18 | 2020-08-18 | 一种rna序列中m5c位点的识别预测方法及系统 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN111951889B (zh) |
Families Citing this family (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN115831224B (zh) * | 2022-11-09 | 2024-05-03 | 内蒙古大学 | 一种预测微生物益生潜力的方法及其装置 |
Citations (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101710362A (zh) * | 2009-12-10 | 2010-05-19 | 浙江大学 | 一种基于支持向量机的microRNA靶位点预测的方法 |
KR101593045B1 (ko) * | 2014-11-12 | 2016-02-12 | 인하대학교 산학협력단 | 결합 상대방을 고려하여 dna 서열에서 단백질과 결합하는 부위를 예측하는 방법 |
KR20180017827A (ko) * | 2016-08-11 | 2018-02-21 | 인하대학교 산학협력단 | 염기 프로파일과 조성을 이용하여 단백질과 결합하는 rna 서열 영역을 예측하는 방법 및 시스템 |
DE102017002092A1 (de) * | 2017-03-04 | 2018-09-06 | Johannes-Gutenberg-Universität Mainz | Verfahren zur Detektion von bekannten Nukleotid-Modifikationen in einer RNA |
CN109859798A (zh) * | 2019-01-21 | 2019-06-07 | 桂林电子科技大学 | 一种细菌中sRNA与其靶标mRNA相互作用的预测方法 |
CN110379464A (zh) * | 2019-07-29 | 2019-10-25 | 桂林电子科技大学 | 一种细菌中dna转录终止子的预测方法 |
CA3107649A1 (en) * | 2018-08-08 | 2020-02-13 | Deep Genomics Incorporated | Systems and methods for determining effects of therapies and genetic variation on polyadenylation site selection |
CN111161793A (zh) * | 2020-01-09 | 2020-05-15 | 青岛科技大学 | 基于stacking集成的RNA中N6-甲基腺苷修饰位点预测方法 |
-
2020
- 2020-08-18 CN CN202010832292.0A patent/CN111951889B/zh active Active
Patent Citations (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101710362A (zh) * | 2009-12-10 | 2010-05-19 | 浙江大学 | 一种基于支持向量机的microRNA靶位点预测的方法 |
KR101593045B1 (ko) * | 2014-11-12 | 2016-02-12 | 인하대학교 산학협력단 | 결합 상대방을 고려하여 dna 서열에서 단백질과 결합하는 부위를 예측하는 방법 |
KR20180017827A (ko) * | 2016-08-11 | 2018-02-21 | 인하대학교 산학협력단 | 염기 프로파일과 조성을 이용하여 단백질과 결합하는 rna 서열 영역을 예측하는 방법 및 시스템 |
DE102017002092A1 (de) * | 2017-03-04 | 2018-09-06 | Johannes-Gutenberg-Universität Mainz | Verfahren zur Detektion von bekannten Nukleotid-Modifikationen in einer RNA |
CA3107649A1 (en) * | 2018-08-08 | 2020-02-13 | Deep Genomics Incorporated | Systems and methods for determining effects of therapies and genetic variation on polyadenylation site selection |
CN109859798A (zh) * | 2019-01-21 | 2019-06-07 | 桂林电子科技大学 | 一种细菌中sRNA与其靶标mRNA相互作用的预测方法 |
CN110379464A (zh) * | 2019-07-29 | 2019-10-25 | 桂林电子科技大学 | 一种细菌中dna转录终止子的预测方法 |
CN111161793A (zh) * | 2020-01-09 | 2020-05-15 | 青岛科技大学 | 基于stacking集成的RNA中N6-甲基腺苷修饰位点预测方法 |
Non-Patent Citations (4)
Title |
---|
一种新的融合统计特征的DNA甲基化位点识别方法;孙佳伟;张明;王长宝;徐维艳;程科;段先华;;江苏科技大学学报(自然科学版)(02);全文 * |
基于SVM和优化特征集的MicroRNA靶标预测;王宝文;齐晓阳;王常武;刘文远;司亚利;;生物医学工程学杂志(06);全文 * |
基于多特征融合的基因调控网络构建方法研究;孟军;周广博;黄楚冰;;小型微型计算机系统(04);全文 * |
基于机器学习的microRNA预测方法研究进展;王颖;李金;王磊;徐成振;才忠喜;;计算机科学(02);全文 * |
Also Published As
Publication number | Publication date |
---|---|
CN111951889A (zh) | 2020-11-17 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN111161793B (zh) | 基于stacking集成的RNA中N6-甲基腺苷修饰位点预测方法 | |
Ranawana et al. | A neural network based multi-classifier system for gene identification in DNA sequences | |
Rostami et al. | A clustering based genetic algorithm for feature selection | |
CN107463799B (zh) | 交互融合特征表示与选择性集成的dna结合蛋白识别方法 | |
CN111951889B (zh) | 一种rna序列中m5c位点的识别预测方法及系统 | |
Hussein et al. | Deep learning and machine learning via a genetic algorithm to classify breast cancer DNA data | |
CN117252114B (zh) | 一种基于遗传算法的电缆耐扭转实验方法 | |
Baten et al. | Fast splice site detection using information content and feature reduction | |
Golenko et al. | IMPLEMENTATION OF MACHINE LEARNING MODELS TO DETERMINE THE APPROPRIATE MODEL FOR PROTEIN FUNCTION PREDICTION. | |
Platon et al. | IRSOM, a reliable identifier of ncRNAs based on supervised self-organizing maps with rejection | |
Lee et al. | Survival prediction and variable selection with simultaneous shrinkage and grouping priors | |
CN116153396A (zh) | 一种基于迁移学习的非编码变异预测方法 | |
Arshak et al. | A new dimensional reduction based on cuttlefish algorithm for human cancer gene expression | |
Iqbal et al. | Data mining of protein sequences with amino acid position-based feature encoding technique | |
CN111755074B (zh) | 一种酿酒酵母菌中dna复制起点的预测方法 | |
CN111599412B (zh) | 基于词向量与卷积神经网络的dna复制起始区域识别方法 | |
Akey Sungheetha | An efficient clustering-classification method in an information gain NRGA-KNN algorithm for feature selection of micro array data | |
Sengupta et al. | A scoring scheme for online feature selection: Simulating model performance without retraining | |
CN112801197A (zh) | 一种基于用户数据分布的K-means方法 | |
Bonidia et al. | Selecting the most relevant features for the identification of long non-coding RNAs in plants | |
CN114155910B (zh) | 一种癌症体细胞突变功能影响预测方法 | |
Alzubaidi et al. | A new hybrid global optimization approach for selecting clinical and biological features that are relevant to the effective diagnosis of ovarian cancer | |
Bhat et al. | OTU clustering: A window to analyse uncultured microbial world | |
Mo et al. | Applications of Machine Learning in Phylogenetics | |
CN103646159B (zh) | 一种基于约束性布尔网络的最大评分预测方法 |
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 |