为了构建一个高质量的与剪接突变相关(Splicing-Associated Variant, SAV)的数据集, 首先整合人类基因突变数据库(Human Gene Mutation Database, HGMD专业版2022.4)^[34]上与剪接相关的突变数据, 共得到34 852条剪接突变数据, 其中包含突变类型为splice、exonic-splice、canonical-splice的数据.同时移除在Y染色体上的唯一一条剪接突变数据, 最终仅保留染色体1~22和X上的突变数据.然后, 根据人类参考基因组版本为GRCh37和致病性标签为DM(Disease Causing Mutation), 对剪接突变数据进行筛选, 最终获得26 550条致病剪接突变数据.另外, 还从ClinVar(https://www.ncbi.nlm.nih.gov/clinvar, 发布于2022年7月9日的VCF文件)^[35]中获取402 763条临床意义为Benign的突变数据.

由于本文的致病剪接突变数据仅包含单核苷酸突变, 因此进一步筛选突变类型为单核苷酸的良性突变数据(GRCh37版本).同时, 为了保证负样本数据的置信度较高, 选择ReviewStatus为criteria provi-ded, multiple submitters, no conflicts的突变数据, 并限定染色体为1~22和X.

经过上述筛选, 最终得到32 064条良性突变数据.根据文献[27], 将剪接区域定义为外显子及其周围50 bp, 并筛选位于剪接区的突变数据, 从而获得30 772条良性剪接突变数据.

致病剪接突变与良性剪接突变数据的来源不同, 导致存在少量数据标签矛盾.因此, 移除正负样本标签不一致的数据, 最终得到26 543条致病剪接突变数据和30 765条良性剪接突变数据.另外, 鉴于非经典剪接突变的临床意义较难解释, 根据剪接区域的不同(如表1所示), 进一步将剪接突变数据划分到经典剪接区(Canonical Splicing Region)与非经典剪接区(Non-Canonical Splicing Region).

表1

Table 1

表1(Table 1)

表1 数据集在不同剪接区的分布情况 Table 1 Dataset distribution in different splicing regions 样本类别 经典剪接区 非经典剪接区 总计 致病(正样本) 21965 4578 26543 良性(负样本) 25 30740 30765

表1 数据集在不同剪接区的分布情况 Table 1 Dataset distribution in different splicing regions

为了减少来自同一基因的不同突变造成的偏差, 采用Cytogenetic Band划分训练集和测试集.首先, 从总共721个染色体条带中随机选择20%(144个)的染色体条带, 将这些条带中包含的突变作为测试集.然后, 将剩下的80%(577个)染色体条带中包含的突变作为训练集.染色体条带是指突变位点所处的染色体区段, 每个突变都与特定的染色体区段关联.该方法的目的是减少位于同一基因中的不同突变同时被用于训练和测试的几率, 因为相邻的突变可能具有相似的特征, 可能导致结果偏差^[19].通过Cytogenetic Band, 可确保同一染色体区段内的突变要么用于训练, 要么用于测试, 而不会同时用于训练和测试.划分后得到包含45 999条数据的训练集和包含11 309条数据的测试集, 具体如表2所示.

表2

Table 2

表2(Table 2)

表2 数据集划分情况 Table 2 Division of dataset 数据集 致病突变 良性突变 染色体条带数 总计 训练集 21224 24775 577 45999 测试集 5319 5990 144 11309

表2 数据集划分情况 Table 2 Division of dataset

本文使用PyTorch框架构建BCsplice.模型以DNABERT为基础架构, 采用其默认参数, 详细参数见表3.对于输入的k-mer序列, 文献[35]研究表明, 当k=6时, 模型具有最优性能, 因此设置k=6, 即每6个碱基作为一个词.

表3

Table 3

表3(Table 3)

表3 BCsplice参数信息 Table 3 Parameters of BCsplice 模块 参数 数值 输入模块 序列长度 k 200 6 BERT模块 隐藏层数量 注意力头数 隐层单元数 12 12 768 CNN模块 卷积核大小 卷积核数量 2× 768 3 训练模块 批次大小 迭代次数 学习率 预热学习率 随机失活率 权重衰退 32 5 2e-4 0.1 0.1 0.01

表3 BCsplice参数信息 Table 3 Parameters of BCsplice

本文使用7个常用指标评估BCsplice分别在全基因组区域和非经典剪接区的性能表现.具体指标如下:Accuracy(ACC)、ROC(Receiver Operating Cha-racteristic Curve)曲线下面积(Area under the ROC Curve, AUC), 精确率-召回率曲线下面积(Area under the Precision-Recall Curve, AUPR)、F1 score(F1)、马修斯相关系数(Matthews Correlation Coeffi-cient, MCC)、精确率(Precision)、召回率(Recall).这些评价指标是根据混淆矩阵计算的, 参数包括真正例(True Positive, TP)、真反例(True Negative, TN)、假正例(False Positive, FP)和假反例(False Nega-tive, FN), 具体计算方式定义如下:

$\begin{array}{l} A C C=\frac{T P+T N}{T P+F P+T N+F N}, \\ F 1=2\left(\frac{\text { Precision } \cdot \text { Recall }}{\text { Precision }+ \text { Recall }}\right), \\ M C C= \\ \frac{T P \cdot T N-F P \cdot F N}{\sqrt{(T P+F P)(T P+F N)(T N+F P)(T N+F N)}}, \\ \text { Precision }=\frac{T P}{T P+F P}, \\ \text { Recall }=\frac{T P}{T P+F N} .\end{array}$

一般来说, 各指标值越高, 模型性能越优.在进行性能对比时, 主要考虑AUC与AUPR, 因为在这7个评价指标中, AUC与AUPR是最重要的2个综合评价指标.

经典剪接区的突变位于外显子与内含子边界± 2 bp位置, 具有明显的基序特征和位置特征, 通常容易解释.非经典剪接区的突变相对更隐蔽, 致病性较难确定.因此, 本文将突变分为2个区域:经典剪接区和非经典剪接区.外显子和内含子边界的GT与AG二核苷酸位置的突变是发生在经典剪接区的突变, 这个位置的突变绝大多数都是致病的^[16].在本文设计的数据集上, 经典剪接区的致病突变有21 965条, 而良性突变只有25条.由此可见, 经典剪接区上的良性突变极少, 仅占0.08%, 良性突变基本上都位于非经典剪接区.此外, 根据本文的分类结果, 经典剪接区的正负样本非常不平衡.尽管非经典剪接区的正负样本也存在不平衡, 但相比之下, 在非经典剪接区发生的突变并不像经典剪接区那样直观, 其致病性需要进一步研究.基于不同区域的致病性情况和突变数据的分布情况, 在接下来的性能对比中主要关注非经典剪接区和全基因组区域(包括经典剪接区和非经典剪接区).测试集在非经典剪接区和全基因组区域上的具体分布情况如表4所示.

表4

Table 4

表4(Table 4)

表4 测试集在不同剪接区的分布情况 Table 4 Distribution of test set in different splicing regions 剪接区 致病突变 良性突变 总计 非经典 916 5985 6901 全基因组 5319 5990 11309

表4 测试集在不同剪接区的分布情况 Table 4 Distribution of test set in different splicing regions

为了研究BCsplice中BERT的隐层向量经过TextCNN处理后对模型性能提升的重要性, 定义BERT-TextCNN和BERT, 进行针对TextCNN模块的消融实验.由于经典剪接突变通常会对剪接产生有害影响, 故预测非经典剪接突变的致病性更具挑战性.因此, 除了在全基因组区域进行分析以外, 还对非经典剪接区进行单独分析.在列出的评价指标中, AUC与AUPR是最重要的两个综合评价指标, 因此主要根据这两个指标进行性能对比.

在全基因组区域, 首先对模型进行五折交叉验证和测试数据集上的性能评估, 具体如表5所示.根据表中结果可见, BERT-TextCNN在各个评价指标上都有性能提升.而在非经典剪接区域, BERT-TextCNN和BERT的指标值对比如表6所示.由表可见, BERT-TextCNN性能提升更显著, BERT-Text-CNN的AUC值为0.947, AUPR值为0.828, 相比BERT分别提高5.1%和15.2%.

表5

Table 5

表5(Table 5)

表5 各模型在五折交叉验证和测试集上的性能 Table 5 Performance of different models on 5-fold cross-validation and test set 模型 数据集 ACC F1 MCC Precision Recall AUC AUPR BERT- TextCNN 训练集 测试集 0.953± 0.008 0.959 0.949± 0.010 0.957 0.905± 0.017 0.918 0.944± 0.010 0.954 0.954± 0.010 0.960 0.985± 0.004 0.987 0.984± 0.004 0.988 BERT 训练集 测试集 0.942± 0.008 0.935 0.938± 0.008 0.931 0.883± 0.016 0.870 0.932± 0.008 0.934 0.945± 0.009 0.928 0.979± 0.005 0.975 0.978± 0.005 0.977

表5 各模型在五折交叉验证和测试集上的性能 Table 5 Performance of different models on 5-fold cross-validation and test set

表6

Table 6

表6(Table 6)

表6 各模型在非经典剪接区测试集上的性能 Table 6 Performance of models on test set of non-canonical splicing region 模型 ACC F1 MCC Precision Recall AUC AUPR BERT- TextCNN 0.940 0.785 0.752 0.753 0.821 0.947 0.828 BERT 0.906 0.656 0.602 0.639 0.674 0.896 0.676

表6 各模型在非经典剪接区测试集上的性能 Table 6 Performance of models on test set of non-canonical splicing region

每个编码层的句子向量表示不同级别的语义信息.浅层的编码层能更好地捕捉低级别的语义信息, 而深层的编码层能更好地捕捉高级别的语义信息.通过TextCNN将多个编码层的语义信息进行组合提取, 既能获取原始的全局语义信息, 又能获得序列信息中关键的特征.这可能是模型性能提升的重要原因之一.此外, 非经典剪接区突变的影响更隐晦, 对于序列上下文的深层语义信息更依赖.

因此, 提取多级别的隐含信息对于预测非经典剪接突变的致病性具有重要意义.

为了进一步验证BCsplice的优越性, 选择如下剪接突变预测方法进行对比:SQUIRLS^[19]、Spli-ceAI^[22]、CADD-Splice^[25]、SPiP^[26]、S-CAP^[27]、MLC-splice(Machine Learning-Based Classification of Splice Sites Variants)^[36].

SpliceAI是一种可从pre-mRNA转录本序列中预测剪接连接(Splice Junctions)的方法, 通过对比突变前后序列的分数差异以预测突变的致病性.CADD-splice整合其它工具的剪接分数, 并优化对剪接突变效应的预测.MLCsplice是一个集成的机器学习方法, 用于预测突变对剪接的影响.SPiP也是一种机器学习工具, 用于检测外显子和内含子突变对mRNA剪接的影响.

为了对比方法的性能, 使用11 309条测试数据, 并在全基因组区域和非经典剪接区域进行实验, 结果如表7和表8所示.

表7

Table 7

表7(Table 7)

表7 各方法在全基因组区上的性能对比 Table 7 Performance comparison of different methods in whole genome region 类别 方法 ACC F1 MCC Precision Recall AUC AUPR 存在缺失值 CADD-Splice SpliceAI MLCsplice S-CAP 0.914 0.962 0.965 0.936 0.860 0.937 0.965 0.918 0.800 0.910 0.932 0.874 0.821 0.924 0.983 0.851 0.904 0.949 0.948 0.996 0.962 0.983 0.984 0.926 0.940 0.975 0.982 0.740 无缺失值 SPiP SQUIRLS BCsplice 0.951 0.925 0.959 0.945 0.921 0.957 0.903 0.850 0.918 0.985 0.906 0.954 0.909 0.938 0.960 0.963 0.974 0.987 0.974 0.979 0.988

表7 各方法在全基因组区上的性能对比 Table 7 Performance comparison of different methods in whole genome region

表8

Table 8

表8(Table 8)

表8 各方法在非经典剪接区上的性能对比 Table 8 Performance comparison of different methods in non-canonical splicing region 类别 方法 ACC F1 MCC Precision Recall AUC AUPR 存在缺失值 CADD-Splice SpliceAI MLCsplice S-CAP 0.897 0.956 0.958 0.907 0.441 0.709 0.936 0.520 0.406 0.689 0.906 0.558 0.349 0.645 0.965 0.353 0.599 0.787 0.908 0.981 0.835 0.930 0.977 0.914 0.446 0.771 0.965 0.227 无缺失值 SPiP SQUIRLS BCsplice 0.933 0.880 0.940 0.696 0.594 0.785 0.679 0.528 0.752 0.879 0.540 0.753 0.576 0.660 0.821 0.855 0.864 0.947 0.735 0.705 0.828

表8 各方法在非经典剪接区上的性能对比 Table 8 Performance comparison of different methods in non-canonical splicing region

在本文的测试数据中, 并不是每个工具都能对所有突变给出预测结果.因此将对比预测工具分为两类:1)不能对所有突变给出预测结果的工具, 即存在缺失突变(CADD-splice、SpliceAI、S-CAP、MLC-splice); 2)可对所有测试突变给出预测结果, 即没有缺失突变(SPiP、SQUIRLS).在表中使用“ 类别” 进行标注以区分这两类工具.需要特别注意的是, BCsplice也属于没有缺失突变的这一类工具.因此BCsplice与SPiP、SQUIRLS可在同一水平上进行对比.

在全基因区域, 使用由5 319个致病剪接突变和5 990个良性突变组成的测试集.根据表7的评估结果显示, BCsplice在AUC和AUPR指标上取得最高值(AUC=0.987, AUPR=0.988).相比存在缺失突变的工具, BC-splice表现出明显的优势.BCsplice能对所有突变数据进行预测, 而其余4种方法都存在一定的缺失突变.具体缺失突变的数量见表9.MLCsplice缺失的突变数量最多, 达到10 092条(10092÷ 11309≈ 89%).其次是S-CAP, 缺失5 076条突变(5076÷ 11309≈ 45%).相比没有缺失突变的工具, BCsplice在AUC和AUPR值上的提升约为1%~2%(见图2), 也显示出一定的优势.

表9

Table 9

表9(Table 9)

表9 各方法在全基因组区上的混淆矩阵值 Table 9 Confusion matrix values of different methods in whole genome region 方法 TN TP FN FP 未返回 预测结果 BCsplice 5741 5108 211 249 0 SPiP 5916 4834 485 74 0 SQUIRLS 5471 4987 332 519 0 CADD-splice 5482 2237 238 489 2863 SpliceAI 5725 2346 126 192 2920 S-CAP 3613 2223 9 388 5076 MLCsplice 592 583 32 10 10092

表9 各方法在全基因组区上的混淆矩阵值 Table 9 Confusion matrix values of different methods in whole genome region

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 各方法在全基因组区域上的综合性能对比Fig.2 Comprehensive performance comparison of different methods in whole genome region

综合来看, 相比对比方法, BCsplice性能更优或相当.这是因为对于使用者来说, 一个理想的方法不仅需要高精确度, 还需要能分析更广泛范围的突变.BCsplice可对全基因组上任何区域的突变进行预测, 而当前的预测工具大多存在剪接突变预测分数缺失的情况.那些能对所有剪接突变进行预测的工具, 在性能上并未高于BCsplice, 特别是在非经典剪接区域中.

另外需要说明的是, 由于获取MLCsplice和S-CAP的预测分数需要先计算其它多种工具的特征分数, 过程较耗时, 因此本文直接使用预先计算好的MLCsplice和S-CAP的人类基因组突变分数, 这可能是两个预测工具缺失值较多的原因.对于CADD-Splice, 本文使用线上预测功能, 下载包含PHRED标注突变致病性评分的tsv文件.为了评估SpliceAI, 根据其文献中提供的网站资源, 在服务器上下载SpliceAI工具包, SpliceAI为每个突变提供4个Δ 分数, 这4个分数中的最大值代表该突变的Δ Score.然而, 对于CADD-Splice和SpliceAI, 目前无法明确它们存在缺失值的具体原因.可能是由于它们依赖于某些剪接特征或受到基因组区域的限制, 导致无法预测所有的剪接突变.对于所有对比方法, 在本文的测试集上, 使用find_optimal_cutoff函数生成具有高AUC的最佳阈值, 确保性能对比的准确性与一致性.

非经典剪接突变的临床意义较难解释, 发生在非经典剪接区的突变也难以识别.大量的预测方法专注于经典剪接区, 对非经典剪接区突变的预测能力较弱.因此本文在非经典剪接区上也进行性能对比, 具体如表10所示.在本文的测试集上, 非经典剪接区的致病剪接突变有916条, 良性突变有5 985条.

表10

Table 10

表10(Table 10)

表10 各方法在非经典剪接区上的混淆矩阵值 Table 10 Confusion matrix values of different methods in non-canonical splicing region 方法 TN TP FN FP 未返回 预测结果 BCsplice 5738 752 164 247 0 SPiP 5912 528 388 73 0 SQUIRLS 5469 605 311 516 0 CADD-splice 5480 261 175 486 499 SpliceAI 5725 343 93 189 551 S-CAP 3610 212 4 388 2687 MLCsplice 592 277 28 10 5994

表10 各方法在非经典剪接区上的混淆矩阵值 Table 10 Confusion matrix values of different methods in non-canonical splicing region

由表10可见, 相比预测区域在全基因组区域, 所有方法在非经典剪接区上都显示出一定程度的性能下降.其中CADD-Splice、SPiP与SQUIRLS性能下降最明显, 这可能与非经典剪接区的样本数据不平衡有关.在非经典剪接区域, 本文依旧将对比方法分为两类.相比存在缺失突变的方法, MLCsplice在AUC和AUPR指标上具有最高值.然而, 在本文的测试集上, MLCsplice存在5 994个突变的缺失情况(占总数的87%), 具体缺失情况见表10.相比没有缺失突变的方法, BCsplice具有显著优势(AUC=0.947, AUPR=0.828).相比SPiP(AUC=0.855, AUPR=0.735)与SQUIRLS(AUC=0.864, AUPR=0.705), BCsplice的综合性能提升8%~12%, 具体如图3所示.

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 各方法在非经典剪接区上的综合性能对比Fig.3 Comprehensive performance comparison of different methods in non-canonical splicing region

因此, BCsplice在不同情况下都表现出独特的优势.非经典剪接突变通常更依赖序列上下文的深层语义信息.BCsplice通过结合CNN和BERT以捕捉多级别的语义信息, 这可能是其性能较优的原因之一.