目录
论文简介
作者单位:
第一作者:大连海事大学理学院
通讯作者:莫纳什大学生物医学发现研究所 (BDI)癌症与感染与免疫项目的副教授和组长
论文开发WEB: http://flagship.erc.monash.edu/XG-m7G/
发表期刊:《Molecular Therapy-Nucleic Acids》
期刊分区及影响因子如下表:
1. 背景
7-甲基鸟苷(m7G)理解:一种被修饰的嘌呤核碱基。它是一种被甲基化的鸟苷(甲基取代氢原子),带正电的mRNA修饰,对基因的高效表达和细胞活力起着至关重要的作用(包括RNA剪接、聚腺苷酸化和mRNA稳定性),并且当它存在于人类尿液之中时,它可能被作为某种癌症的生物标记,分子结构如下图所示:
2. 数据
阳性样本:对张等人HeLa和HepG2细胞中801个碱基分辨率的m7G位点的测序,获得了801个m7G位点的序列,将m7G位点位于序列中心的序列提取为41bp,上游为20bp,下游为20bp(窗口的大小为41)。然后用CD-HIT18去除冗余序列19,阈值为80%,之后保留741个含有m7G位点的阳性样本。
阴性样本:为了进一步避免数据不平衡带来的Sn偏低的潜在问题,选择了741个序列相似度小于80%的非m7G位点的序列组成最终的阴性样本数据集。
3. 方法
3.1 概述
XG-m7G的流程图如下图所示,XG-m7G的开发五个主要步骤:
- 搜集m7G位点基准数据集,正负样本各741个(RNA序列 )
- 利用六种编码方式对RNA序列进行特征编码( binary
encoding, CKSNAP, ENAC, NCP, ND, and SCPseDNC)。 - 利用SHAP计算特征的重要性,进行特征选择。
- 利用XGBoost计算相关的评价指标。
- 开发Web服务器。
3.3 特征编码方法
3.3.1 Binary Encoding
编码原理:
A:(1,0, 0, 0)
C:(0, 1, 0, 0)
G:(0, 0, 1, 0)
U:(0, 0, 0, 1)
结果:每个大小为41窗口的RNA编码为164维的特征向量
3.3.2 CKSNAP
**编码原理:**由k个残基分隔的核苷酸对的频率,CKSNAP功能包含与核酸对相对应的16个值{AA, AC,AG, …, UG, UU},当K=1时得到如下:
:间隔的一个碱基
N
X
∗
Y
N_{X*Y}
NX∗Y:碱基对
N
T
o
t
a
l
N_{Total}
NTotal:序列中一间隔核酸对的总数
结果: k = 0, 1, 2最终每个大小为41窗口的RNA编码得到48维的特征向量(163)
3.3.3 ENAC
编码原理:基于固定长度的窗口计算核酸组成,该窗口从每个核苷酸序列的5到3个末端连续滑动(可以理解为窗口大小为41里面的子窗口),ENAC可按如下方式计算:
S:滑动窗口的大小
N
t
,
w
i
n
N_{t,win}
Nt,win:窗口内四种碱基{ A; C; G; U}的数量
结果:每个大小为41窗口的RNA编码为160维的特征向量((41-1)*4)
3.3.4 NCP
编码原理:A、C、G和U的每个核苷酸都有不同的化学结构和化学结合特征。NCP编码方案将这四种类型的核苷酸的化学性质合并到表示中。核苷酸可以分为三个不同的组,如下图所示:
结果:
A:(1,1,1)
C:(0,1,0)
G:(1,0,0)
U:(0,0,1)
3.3.5 ND
编码原理:
ND也被称为累积核苷酸频率(ANF),它综合了核苷酸频率信息和每个核苷酸在RNA序列中的分布。RNA序列中第i位的任何核苷酸的密度Di可由下式给出:
||
S
i
S_i
Si||:滑动子串的长度
L:而l为相应定位器的序列位置
将 NCP和ND组合在一起理解:
例子:ACGCGGAUUA编码为: { ( 1 , 1 , 1 , 1 ) , ( 0 , 1 , 0 , 0 . 5 ) , ( 1 , 0 , 0 , 0 . 3 3 ) , ( 0 ,1, 0, 0.5), (1, 0, 0, 0.4), (1, 0, 0, 0.5), (1, 1, 1, 0.29), (0, 0, 1, 0.125),(0, 0, 1, 0.22), (1, 1, 1, 0.3)}
结果:
每个大小为41窗口的RNA编码为164维的特征向量(41*4)
3.3.6 SCPseDNC
f
u
f_u
fu:u=1;2;…;16,是序列中第i个二核苷酸的归一化出现频率
w
w
w:W为权重因子,取值范围为0 ~ 1
⋀
\bigwedge
⋀:理化指标的数量,六个索引 (i.e., Rise (RNA), Roll (RNA), Shift
(RNA), Slide (RNA), Tilt (RNA), Twist (RNA)) 。RNA序列。qj (j = 1;2;;L)是j级相关因子,定义为:
公式中相关函数表示为:
λ
\lambda
λ:沿核苷酸序列的最高计数的相关等级(或层)
u:1; 2; .;
⋀
\bigwedge
⋀
m:1; 2; .;
λ
\lambda
λ
i:1;2; .; L- m- 2
P
u
(
R
i
R
i
+
1
)
P_u(R_iR_{i+1})
Pu(RiRi+1):二核苷酸RiRi+1at位置i的第u (u = 1,2, ., L)理化指标的数值
结果: λ \lambda λ = 20和w = 0.9,我们生成了一个136 维向量
3.4 XGBoost算法
XGBoost算法的基本思想与GBDT类似,不断地地进行特征分裂来生长一棵树,每一轮学习一棵树,其实就是去拟合上一轮模型的预测值与实际值之间的残差。当我们训练完成得到k棵树时,我们要预测一个样本的分数,其实就是根据这个样本的特征,在每棵树中会落到对应的一个叶子节点,每个叶子节点就对应一个分数,最后只需将每棵树对应的分数加起来就是该样本的预测值。
3.5 SHAP算法
SHAP是一个用于解释预测的统一框架,是2017年提出的唯一一种基于预期的一致且局部准确的特征归因方法。这种技术可以从复杂的训练模型中解释特征的重要性分数,并为测试样本提供可解释的预测。SHAP值被提出作为特征重要性的统一度量,因为它们为每个特征分配了一个重要值(fi),表示在模型预测中包含该特征的效果。在合作博弈理论中,SHAP值可计算为:
F:所有特征的集合
S:移除第t个特征后从F获得特征的子集
f
S
U
(
i
)
f_{SU_(i)}
fSU(i):重新训练模型
f
s
f_s
fs:预测模型
X
s
X_s
Xs:集合S中的输入特征值
4 结果
4.1 与其他现今算法比较
比较方法(KNN,SVM,LR,RF)及最优参数设定:
4.1.1 交叉验证
根据上表的参数进行10折交叉验证运行100轮的评价指标如下所示:
结论:
XGBoost算法除sn低于KNN约1.44%外,其他算法评价指标均高于比较的算法。
比较方法10折交差验证的方法ROC曲线如下图所示:
结论:
XGBoost的ROC值均高于其他比较方法。
4.1.2 T检验
采用学生t检验分析了预测结果差异的显著性。从表2可以看出,所有的p值都远小于0.01,说明XGBoost与其他四种算法的差异有统计学意义,),“Win”表示XGBoost超过它的次数,“Draw”表示两者性能相当,“Lost”表示XGBoost更差。
4.2 特征编码对模型预测的影响
4.2.1 最好的20个特征的影响
本研究采用六种不同的特征编码方案来生成特征向量。表S1列出了每种特性的性能。然后,我们使用SHAP算法来刻画特征的重要性和评估特征的行为。具体做法包括一下几个步骤:
- 将所有特性命名为:binary-1, binary-2, ., binary-164;cksnap cksnap - 165 -
166,cksnap - 212。scpsdnc -537, scpsdnc -538,scpsednc - 672。 - 根据式15计算出SHAP值,将所有样本的前20个特征绘制在图3中。
- 在图3A中,每一行代表一个特征,每个点是实例的SHAP值。样本点越红表示特征值越大,样本点越蓝表示特征值越小;横轴代表SHAP值,feature NCP-466所在的区域在图3B中被放大了。
SHAP值的意义: - 如果SHAP值为正,这意味着该特性推动了对m7G位点的预测,并具有积极的影响;如果是负面的,该功能会将预测导向非m7g位点,例如当NCP-432和NCP-438等特征具有较高的SHAP值时,该模型被驱动向正向的m7G预测。
- SCPseDNC-546, SCPseDNC-543, SCPseDNC-567, 和 SCPseDNC-579。促进非m7g位点的预测,
- 而其他功能SCPseDNC-541, SCPseDNC-538, SCPseDNC-599,and SCPseDNC-549仅促进模型对m7G位点的预测 。
4.2.2 去除这最好20个特征后的影响
对每个特征的SHAP绝对值的平均值来对特征的重要性进行排序在补充表格2中(excel文件)。jackknife检验也是统计预测中常用的交叉验证方法。我们保留了前50名,根据对100、150、200、250、300、350、400、450、500、550、600个特征进行排序,通过执行折刀测试找到最佳特征组合,表S3列出了所有性能结果,当特征维数降至150时,模型的整体性能最好。
为了说明特征选择的有效性,我们提供了表S4和图4中使用原始特征和最优特征训练的模型的折刀测试性能对比结果。结果还表明,使用SHAP选择的最优特征训练的模型比使用所有原始特征训练的模型明显提高了预测性能。
4.3 特征选择评价指标的影响
使用F-score和最小冗余-最大相关性(mRMR)进行了相同的特征排序和选择过程,这在生物信息学和计算生物学领域被广泛采用来降低特征维数的变化。
结果:
按照f分数排序后,前400个功能的组合达到了最佳Acc的91.36%,相比之下,经mRMR排序后,前600个特征的组合达到了最佳Acc值的90.35%。这些结果表明了SHAP对m7G位点识别的有效性。
4.4 与iRNA-m7G比较
iRNA-m7G是目前唯一建立的用于在RNA序列中搜索m7G位点的模型。因此,我们通过10倍交叉验证试验比较了XG-m7G和iRNA-m7G的性能。表3列出了XG-m7G和iRNA-m7G在Sn、Sp、Acc、MCC和AUC值。
结果:
- 交叉验证:XG-m7G对Sn和Acc分别提高了2.82%和1.41%,XG-m7G的MCC和AUC分别为0.025和0.026
- jackknife 测试:XG- m7G的Sn为92.17%,Sp为91.77%,Acc为91.97%,MCC为0.839,AUC为0.972。
- 总之,这些结果证实了我们提出的模型XG-m7G在识别m7G位点方面优于iRNA-m7G
5 结论
- 通过10倍交叉验证和折刀试验,比较了XG-m7G与最新的m7G位点识别模型iRNA-m7G的性能。
- 基于模型解释方法SHAP,探讨了最重要的特征,并验证了它们对m7G站点识别的贡献。
- 构建了一个web服务器,它允许感兴趣的用户使用我们的模型来识别m7G站点,并根据自己的数据集方便地训练他们的特定模型.
6. 启发
- 论文特征可以采用shap进行特征选择并进行可视化
- 当没有独立测试集的时候可以采用交差验证和Jackknife来验证算法的效果