论文解读：《ILGBMSH：基于集成学习算法用于shRNA目标预测的可解释分类模型》

文章地址：https://academic.oup.com/bib/article-abstract/23/6/bbac429/6731717?redirectedFrom=fulltext&login=false
DOI：https://doi.org/10.1093/bib/bbac429
期刊：BRIEFINGS IN BIOINFORMATICS
2022年影响因子/JCR分区：13.994/Q2
发布时间：2022年9月30日
Github：https://github.com/ChengkuiZhao/ILGBMSH

1.文章概述

短发夹RNA (Short hairpin RNA, shRNA)介导的基因沉默是实现RNA干扰的重要技术，设计有效可靠的shRNA分子起着至关重要的作用。然而，通过生物技术高效地选择shRNA靶标是昂贵且耗时的。因此，开发一种更精确和有效的计算方法来设计有效和可靠的shRNA分子是至关重要的。在这项工作中，作者提出了一个可解释的分类模型，用于shRNA目标预测，使用（Light Gradient Boosting Machine，ILGBMSH）。作者并没有仅仅利用shRNA序列特征，而是提取了554个生物学和深度学习特征，这些特征在之前的shRNA预测研究中没有考虑到。与最先进的shRNA目标预测模型进行了比较，评估了模型的性能。此外，还从模型的参数和可解释方法Shapley Additive explanation（SHAP）来研究特征解释，这为从模型中获得生物学见解提供了依据。作者使用来自其他资源的独立shRNA实验数据来证明所提出模型的预测能力和鲁棒性。最后，利用模型设计了miR30-shRNA序列，并进行了基因敲低实验。实验结果与预期完全一致，Pearson相关系数为0.985。综上所述，ILGBMSH模型可以实现最先进的shRNA预测性能，并从机器学习模型参数中提供生物学见解。

2.关键点

虽然有许多工具可以预测有效的siRNA靶点，但由于生物成因的差异，它们不能用于选择有效的shRNA靶点。
作者提出了一个计算框架ILGBMSH，以预测具有miR30主干的有效shRNA靶点。
通过结合LightGBM模型训练的生物和深度学习特征，实现了shRNA预测的最佳性能。
基于 SHAP 的预测结果和可解释方法，研究 shRNA 函数中的特征重要性。
作者进行了一个生物实验来验证模型的预测能力和稳健性，这使得它成为指导shRNA实验设计的有用工具。

3.背景

设计用于基因沉默的人工分子是抑制RNA表达的一种有前途的方法。具体来说，小干扰RNA (siRNA)和短发夹RNA (shRNA)是导致有效基因敲低的两种主要人工分子。SiRNA是合成的短干扰双核苷酸，可直接转染到细胞细胞质中形成RNA诱导沉默复合体;然后，这个复合物被进一步用于使目标RNA沉默。shRNA从DNA载体转录后进入内源性RNA干扰(endogenous RNA interference, RNAi)通路，形成有效的siRNAs。虽然siRNA可以在实验室中方便地获得，并被广泛用于基因敲低，但它在细胞中只能持续几天。这限制了它的科学研究和治疗。与合成siRNA相比，shRNA可以提供同样有效的性能，同时允许更长的稳定敲低，这导致了更大范围的潜在应用。根据以往的研究，shRNA可能引起较少的脱靶效应，有助于避免毒性细胞反应。shRNA分子也可以嵌入到microRNA的主干中，以一种更自然的生物方式发挥作用。此前，由于我们对RNAi处理的了解有限，缺乏支持算法，很难精确地得到基于microrna的效能高、脱靶效果低的最佳shRNA。近年来，随着实验研究人员生成越来越多的shRNA数据，基于microrna的shRNA的计算设计是一种很有前途的方法，可以潜在地提高基因沉默性能并减少不必要的副作用。
许多计算工具已经开发用于siRNA靶效预测。由于siRNA和shRNA在生物学上的相似性，一些研究人员利用siRNA预测工具选择有效的shRNA靶点来实现RNA沉默。但由于这两种分子在细胞内的生物学加工方式不同，siRNA预测工具所选择的成熟序列在shRNA功能中可能没有意义。因此，有必要开发一种预测有效shRNA的计算算法。虽然shRNA预测工具的开发受到实验数据缺乏的限制，但Fellmann等人通过构建和评估约20000个基于microrna的shRNA来填补这一空白，提供了足够大的数据集来训练shRNA机器学习模型。基于该数据集和回归算法，Matveeva等人开发了“miR_Scan”软件，以选择最佳的基于mi30的shRNA。另一项工作使用支持向量机(SVM)分类算法预测了另一个基于mir-30的大规模数据集，该数据集包含以前的大约20 000个shRNA和大约10 000个新的有效shRNA。结果表明，与siRNA模型相比，shRNA预测模型可以显著提高shRNA预测精度。
而最佳shRNA回归算法“miR_Scan”的Pearson相关值为0.46，最佳shRNA分类算法“SplashRNA”的精度-召回曲线(AUPRC)下面积为0.696。因此，目前shRNA的预测精度并不令人满意，还有很大的提高空间。为此，作者使用最先进的算法深入研究shRNA有意义的生物学和深度学习特征，并改进了有效的shRNA预测。在本研究中，基于Light Gradient Boosting Machine (LightGBM)的先进集成算法开发了一个预测模型。LightGBM已被证明在数千个样本上的基因组选择预测上优于其他算法(ANN, KNN, SVR等)，Pearson相关性约为0.9，并且与其他Gradient Boosting算法相比，它消耗的CPU内存和时间最少。使用LightGBM预测具有生物学和深度学习特征的有效shRNA，包括shRNA序列特征、热力学特征、目标f层区域特征等。作者使用可解释方法Shapley Additive Explanations(SHAP)来评估这些特征的重要性，该方法可以通过计算每个特征对预测的贡献来解释任何模型。然后，将模型的性能与之前的工作和其他主流深度学习算法进行了比较，如CNN和LSTM单元。最后，使用独立的shRNA数据集来验证模型的鉴别能力。

4.数据

shRNA数据集检索自两个已发表的基于mir30的shRNA实验。来自Functional identification of optimized RNAi triggers using a massively parallel sensor assay的第一个RNAi传感器检测(称为“TILE”)数据集包含9个哺乳动物基因的18720个基于mir30的shrna，其中5736个有效序列和12685个弱有效序列。来自Prediction of potent shRNAs with a sequential classification algorithm的第二个数据集(称为“M1”)包含来自333个新基因(332个小鼠基因和1个大鼠基因)的9602个有效shRNA序列。对于“M1”数据集，只向训练数据集中添加了阳性数据，因为这样会得到一个更平衡的数据集，而平衡的数据集可以给我们最好的模型性能。去除无相关标签或相应基因特征的序列，得到一个大型平衡数据集，包含24032条shRNA序列，其中有效shRNA序列11 529条，阴性shRNA序列12 503条。该数据集用于与之前的si/shRNA算法进行性能比较。对于最终版本的数据集，作者进一步删除了306个重复，得到了23 726个shRNA序列，这有助于我们训练一个更广义的预测模型。
对于外部独立的shRNA数据集，作者从Multi‑target inhibition by four tandem shRNAs embedded in homo‑ or hetero‑miRNA backbones收集了5个具有miR30主干的有效shRNA，从自己实验室收集了6个具有一般主干的有效shRNA。利用训练好的模型对PD1基因的所有shRNA序列效率进行了预测。最后，设计了5个mir30-shrna，其中4个预测得分最高(0.95、0.87、0.86、0.85)，最后一个预测得分中等(0.195)，用于进行基因敲低实验。每次体外shRNA实验都有相应的阴性对照序列。

5.方法

作者开发了一个用于shRNA预测的分类模型。训练数据来自两个独立的实验，根据实验得分分布分为敲低效率高组和敲低效率低组。作者使用并比较了两种策略来开发机器学习模型。第一个是提取生物和深度学习特征，可以帮助模型更好更快地理解数据。第二个是让算法自动从原始序列中学习特征，比如深度学习算法。作者选择LightGBM进行生物和深度学习特性的训练，并选择LSTM、BiLSTM、CNN和CNN-LSTM网络作为深度学习算法。采用保留一个基因的交叉验证来评价模型的性能。

5.1 特征提取

根据shRNA特征和mRNA靶区特征以及CNN潜在节点三个主要方面总结了8种生物和深度学习特征。从shRNA序列中提取了单核苷酸、二核苷酸和三核苷酸的特征，通过滑动窗口算法分别得到22个、21个和20个特征。从22个核苷酸中共提取了63个shRNA序列特征。

还总结了单核苷酸、二核苷酸和三核苷酸的84个shRNA序列频率特征，根据四种核苷酸的组合，分别得到4个、16个和64个特征。
使用来自Using an RNA secondary structure partition function to determine confidence in base pairs predicted by free energy minimization的算法计算shRNA二级结构，该算法计算shRNA核苷酸与序列中所有其他核苷酸之间的配对概率。最后，获得了24个结构特征，包括22个shRNA核苷酸的配对概率以及它们的和和最大值。
使用了shRNA双相的两个重要热力学特征：(1)shRNA双相吉布斯自由能之和(ΣΔG);(2) shRNA正义和反义链之间的末端五核苷酸ΔΔG。
mRNA靶区定义为包含22-nt shRNA靶标的序列，靶标两侧各有20个核苷酸。与22-nt靶序列相比，f区靶mRNA能提供更多影响结合和敲低的信息。采用与shRNA序列特征相同的表示方法提取目标区域序列特征。得到了183个序列表示特征和84个频率特征。此外，还总结了64个靶区结构特征。
在经过良好训练的CNN模型的最后一层中获得50个节点作为深度学习特征。

5.2 解释模型预测

1. 单独使用每一类特征来训练模型，进行了保留一个基因的交叉验证；
2. 对于训练好的模型，使用LightGBM参数来评估每个单一特征的贡献程度；
3. 采用SHAP方法对模型预测进行解释。该方法不仅可以直观地显示每个特征的重要性，而且可以显示该特征对模型预测的消极或积极影响。

5.3 实验

作者设计了两个载体：PD1载体和shRNA载体。荧光素酶通过2A肽与靶基因PDCD1结合，显示靶基因是否被敲低，载体上也有蓝色的f荧光蛋白，显示载体是否被转导到靶细胞。使用PerCP/ cyanine5.5偶联的抗人CD279 (PDCD-1) (BioLegend, San Diego, CA, USA)验证了稳定Jurkat细胞中PDCD1的残留表达。shRNA载体包含shRNA和绿色荧光蛋白，也用于识别转导。将载体转入293T细胞系，根据两种载体的比例进行三组实验。作者还准备了空白组，只有PD1载体被转导，对照组是不相关的基因。采用荧光素酶生物发光和f低细胞仪检测靶基因敲低的效率。

6.结果

6.1 分类模型性能比较

使用“TILE”数据集和“M1”正数据集对模型进行训练并验证其预测精度。进行了保留一个基因的交叉验证，这样每次都将来自“TILE”数据集中9个基因中的1个的shrna用作验证数据集，其余的用作训练数据集。将来自所有9个基因的预测shrna放在一起，并计算AUPRC作为最终结果。

6.2 回归模型性能比较

将=的模型的性能与之前的si/shRNA预测方法进行了比较，这些方法都是回归模型，在同一数据集上。为了便于比较，作者用提取的特征训练了一个LightGBM回归模型。在Mysara和Matveeva总结的“TILE”数据集中收集了基于mir30的shrna的预测结果，并将4 /5的数据用于训练，1/5的数据用于测试。

6.3 通过模型解释特征的重要性

前20个重要特性如图4A所示。shRNA序列和RNA f链区域的三核苷酸特征(标记为“shNNN”和“trNNN”)是最具预测性的特征。“shNNN1”和“shNNN19”特征分别代表5和3端shRNA的两端是两个最重要的特征，说明shRNA的两端对敲低效率至关重要。双链shRNA的热力学特征“sigma_g”在特征重要性排名中也排名靠前。
使用SHAP方法来解释模型的预测。训练良好的模型的SHAP结果如图4B所示。

6.4 预测不同骨链的shRNA

作者继续检验模型是否能够准确地预测具有miR30和一般主干的shrna。选择PD1基因进一步验证模型性能，所有shrna预测评分见补充文件1。作者的模型是基于22-nt shRNA序列数据训练的。根据之前的研究，错配在3端shRNA末端耐受较好。因此，使用22-nt shRNA的预测作为对应21-nt shRNA的敲低得分近似值。
从之前的研究中收集了5个基于mir30的有效的PD1基因21-nt shRNAs。通过量化对打乱shRNA对照的抑制率，这5种shRNA都具有很高的敲低功效(>75%)。模型预测了PD1基因的所有2076个shRNA序列。预测分数(范围从0到1)显示强大的击倒概率。然后根据预测得分对所有shRNA序列进行降序排序，并将所有序列中前25%的数字[上分位数(Q3)以上]设为高秩序列。所获得的5个shrna在我们的预测中都是高秩序列，这意味着模型对其他人的实验数据有很好的预测能力(数据见补充文件2)。
然后，将该模型应用于预测具有一般主干的shRna。将6个有效的PD1 shRNA序列组装在总主干中，在之前的实验中敲低率为>80%(数据见补充文件2)。在ILGBMSH预测中，6个有效shRNA中有5个作为有效shRNA在上四分位数中排名较高。这些结果表明，模型具有预测具有其他主干的shrna的能力，并且可以作为选择具有一般主干的有效shrna的可选模型。

6.5 实验验证

为特定基因设计了22-nt shRNAs，并进行了实验进行直接验证。对于PD1基因，共选择2076个22-nt shRNA序列进行模型预测。我们选择了5个shrna，其中4个是排名靠前的序列，最后一个是排名靠后的序列(图5A)。方法部分所描述的实验载体示意图如图5B所示。
在荧光素酶分析中，通过shRNA组与空白组相比发光强度下降的百分比来计算效率。总的来说，在这5个shRNA中，’ shRNA 732 ‘的敲低效率最高(平均= 89.3%)，’ shRNA 736 '的敲低效率最低(平均= 52.3%)，而其他3个shRNA也表现出强敲低效率(平均分别为84.7%、82%、86%)(图5C)。
他的值被设定为参考值。然后，我们将所选shRNA组的PD1 f荧光强度均值与空白组的PD1 f荧光强度均值进行比较，计算敲down效率。
在低细胞分析中，细胞被载体转导24小时。首先，检测了没有shRNA载体的空白组中的平均PD1f荧光强度，并将该值设置为参考值。然后，将所选shRNA组的PD1 f荧光强度均值与空白组的PD1 f荧光强度均值进行比较，计算敲down效率。因此，可以得到每个shRNA的敲低效率，包括不相关基因有shRNA的对照组(补充图S2A)。最后，比较了所选shRNA组与对照组的PD1 f荧光强度，计算出shRNA的相对敲低效率。结果显示了与荧光素酶分析相似的趋势，shRNA 732是最高的(平均= 73.3%)，shRNA 736是最低的(平均= 19.7%)，而其他三个shRNA也表现出强有力的敲低(平均= 67.6%，58%，62.3%)(图5D)。

7.结论

综上所述，作者开发了一种结合多种生物学上有意义的特征来预测基于mir30的shRNA敲除效力的方法，并获得了最佳性能。利用可解释方法研究特征的重要性，为shRNA的生物发生提供了更多的生物学见解。与siRNA预测方法相比，证明了shRNA预测方法在shRNA设计上的优越性。该方法可以为研究人员提供更好的shRNA预测和最优的shRNA设计。