前列腺癌细胞核酸适配体支持向量机分类模型

禹新良

(1.湖南工程学院环境催化与废弃物再生化湖南省重点实验室湖南湘潭 411104；2.湖南大学化学生物传感与计量学国家重点实验室湖南长沙 410082)

核酸适配体(Aptamer)是单链寡核苷酸，它能结合特定的生物靶标，如蛋白质分子、细胞、病毒[1，2]。与抗体相比，核酸适配体具有许多独特的优势，包括合成简单、体外修饰容易、免疫原性小、易于穿透肿瘤组织等[3,4]。核酸适配体的出现，为探索重大疾病早期诊断及靶向治疗方法开辟了崭新的途径[5，6]。因此，核酸适配体相关技术的研发已经成为化学、生物医学等领域的热点课题。

核酸适配体可采用以细胞为靶标的指数富集配体系统进化(Cell-SELEX)技术获得，从人工合成的DNA/RNA文库中筛选得到，即通过重复的寡核苷酸序列与标靶的吸附分配、寡核苷酸序列与标靶复合物洗脱、回收，以及聚合酶链扩增3个基本步骤，从单链寡核苷酸文库中筛选出能与靶标特异性结合的DNA/RNA适配体[7，8]。

成功的适配体SELEX筛选实验通常循环次数多，且每轮循环操作耗时、耗力、耗材[9]。基于分子构效关系(QSAR/QSPR)模型方法可以用于分子性质的估算，然后将最有希望满足某种性能要求的分子用于实验室的合成与测试[10]。因此，核酸适配体分子构-效关系研究可以节约资源、促进新型适配体开发。每年关于分子构效关系的研究报道超过5 000篇[11]，通常这些研究基于的分子参数是直接从研究的分子体系计算而得。由于核酸适配体实验数据的缺乏，尤其是核酸适配体分子结构参数提取的困难，目前对核酸适配体开展构-效关系研究的报道极少。Musafia等人[12]采用简单的拓扑参数，对识别流感病毒的核酸适配体的亲和性进行构-效关系研究，模型训练集、测试集相关系数R2分别为0.702、0.66。另外，我们课题组对肝癌SMMC-7721细胞核酸适配体中心序列形成的二级结构(即环结构)开展计算，将得到参数建立亲和性、专一性的分类模型，其模型的精度高、预测能力强[9]。

前列腺癌(Prostate Carncer，PCa)是最常见的、严重威胁男性健康的恶性肿瘤之一[2,6,13，14]，而目前还缺乏对前列腺癌有效的早期诊断方法。因此，对前列腺癌细胞的核酸适配体开展构-效关系研究，有利于开发适配体分子识别探针，促进适配体应用于前列腺癌早期诊断与治疗[5,6]。本研究拟采用间接分子参数提取模型计算参数，建立以前列腺癌PC-3M-1E8为靶细胞的核酸适配体构-效关系模型。

1 实验部分

1.1 实验数据

本研究采用以前列腺癌PC-3M-1E8为靶细胞，通过Cell-SELEX方法得到候选核酸适配体[6]。单链DNA文库(5′-AGAAGGAAGGAGAGCGACAC-N40-TATCAGTGGTCGGTCGTCAT-3′)制备时，荧光标记的正向引物为：FP-FAM，5′-FAM-AGAAGGAAGGAGAGCGACAC-3′；生物素标记的反向引物为：RP-Biotin，5′-Biotin-ATGACGACCGACCACTGATA-3′。而未做任何标记的、用于克隆测序的正向引物为：FP，5′-AGAAGGAAGGAGAGCGACAC-3′；反向引物为：RP，5′-ATGACGACCGACCACTGATA-3′。所用到的全部序列由上海生工生物工程技术有限公司合成纯化；前列腺高转移癌细胞PC-3M-1E8及前列腺低转移癌细胞PC-3M-2B4均由中国科学院上海科学研究院提供。以前列腺高转移癌细胞PC-3M-1E8进行正筛选，从第4轮开始引入反筛选，选择前列腺低转移癌细胞PC-3M-2B4作为对照细胞，以此减少非特异性吸附，提高筛选效率。

经12轮SELEX筛选，DNA文库与PC-3M-1E8细胞的结合力达到饱和，且文库中核酸序列的富集度高。SELEX实验的第3、5、7、9、11轮产品序列采用高通量测序技术测序。

分别选取第3、5、7、9、11轮产品中，中心序列碱基为40±1，且富集水平最高的各100条序列作为研究对象。其中第3轮中与11轮重复的序列予以剔除，保证第3轮的100条序列与第11轮的100条序列不出现重复。由于对第12轮产品测序的序列少，不具有代表性，且与第11轮序列重复，因此本文没有对12轮序列单独进行分析。

随机地将第3、11轮筛选所得的序列各自的3/4(75条)作为训练集，建立分类模型；各自剩下的1/4(25条)作为测试集，对分类模型进行检验。所得模型同时对第5、7、9轮所得序列进行预测。将第3轮的100条序列类标签值(或目标值)设置为“1”，代表低亲和性、低特异性的候选核酸适配序列；将第11轮筛选所得核酸序列类标签值设置为“2”，代表高亲和性、高特异性候选核酸适配序列。

1.2 分子参数计算

通常，QSAR/QSPR模型所用的分子参数直接从分子结构计算得到。基于DNA(或RNA)核酸适配体碱基序列的一维分子参数在本质上是拓扑参数，与基于分子二维、三维结构参数相比，缺乏对适配体分子的电子、几何等特征的精确反映[9]。由于核酸适配体序列分子较大，直接计算分子二维、三维结构参数存在困难，包括分子结构优化存在困难，以及参数提取存在困难。如最常用的分子参数提取软件Dragon[15]，不能用于核酸适配体(15～45碱基)的分子参数计算。

分子参数的实质是将分子的结构特征转换成数字特征。由于基因碱基序列隐含着氨基酸序列信息，即基因片段结构与相应的氨基酸序列结构相关联。而QSAR/QSPR研究的实质就是建立分子性质与其结构的关系式；核酸适配体性质与其序列结构有关联，而后者又与相应的氨基酸序列结构相关联。我们课题组已经成功采用DNAMAN 6.0软件将抗流感病毒核酸适配体序列翻译成氨基酸序列，计算氨基酸序列的分子参数，然后用这些分子参数建立抗流感病毒核酸适配体亲和性构效关系模型[16]。这种方法用于本文的分子参数计算。

在建立QSAR/QSPR模型时，理论上讲分子参数越精确越好。由于本文参数计算模型是从核酸序列翻译而来，因此这种计算模型本身带有近似性。采用半经验或分子力学等方法对分子进行优化、计算得到分子参数，用于构-效关系研究是可行的。因为采用回归分析时，分子参数中的系统误差可以降低或消除。本文的二元分类模型与回归模型在本质上是一样的，只是分类模型将回归模型的输出离散化。因此，本文在计算氨基酸序列分子参数时，采用ChemBioDraw 11.0构建分子结构，接着使用ChemBio3D 11.0的分子力学(MM2)方法优化，收敛标准为均方根误差小于0.42 kJ/molÅ。然后采用Dragon 6.0软件[15]，为每条氨基酸序列计算4 885个分子参数。将分子参数为常数，或分子参数之间相关性大于0.90的参数予以剔除，最后从每条氨基酸序列得到747个分子参数。

2 结果和讨论

2.1 分类模型建立

采用IBM SPSS Statistics 19的向前逐步(Forward Wald)方法对第3、11轮的200条候选核酸适配体序列类标与相应氨基酸序列分子参数进行二元Logical回归分析，得到分类回归表和变量特征表，分别见表1和表2。表1显示，模型(Model) 7含有7个分子参数(定义见表3)，模型类标(Class)1、类标2的分类准确率分别是79%、81%，平均准确率为80%，模型预测理想。表2显示，模型7中所有变量Sig.值小于0.05，表明7个分子参数均有统计意义。基于上述7个分子参数的二元Logical回归模型，对训练集进行预测。结果显示，对第311轮序列预测总的准确率均为80.0%。

由于本文用“1”、“2”代表核酸适配体亲和性高低，属于典型的二分法研究。而支持向量机分类(SVC)算法被认为是针对小样本统计估计和预测学习的最佳理论，在分类建模中取得了理想的结果[9]，因此本文进一步采用SVC用来建立更为精确分类模型。将从第3、11轮得到的各100条核酸序列随机分成测试集与训练集，每轮中3/4的序列作为训练集，余下的样本作为测试集。用高斯径向基函数(RBF)作为SVC的核函数，模型最佳正则化常数C与Gaussian函数参数γ采用粒子群优化(PSO)算法搜寻。

实验中使用的参数值设置如下：认知学习因子c1初始值和社会学习因子c2初始值设置分别为1.5、1.7；种群最大数量为20；最大迭代次数为200；SVM参数搜索范围C是100～1 000；SVM参数γ搜索范围是0.01～1.0。建模时采用留一法(LOO)进行交叉验证。采用LIBSVM在MATLAB R2012a平台上建立SVC模型。结果显示，相对最佳的SVC模型参数C和γ分别是C=357.7和γ=0.011；SVC模型对训练集的预测准确率为87.3%，其中第3轮核酸序列预测精度为86.7%，第11轮核酸序列的预测准确率为88.0%；对测试集的预测准确率为86.0%，其中第3轮序列预测准确率为88.0%，第11轮序列的预测准确率为84.0%。对第3轮与第11轮序列预测总的准确率均为87.0%。基于相同的数据集，用二元Logical回归分析时，对第3轮与第11轮序列预测总的准确率均为80.0%。很明显，SVC分类模型优于二元Logical回归分类模型。

2.2 模型参数讨论

表3中第1个分子参数LEB属于Burden特征值(Burden Eigenvalues)参数组，该组分子参数基于Burden提出的一种用改进的邻接矩阵计算参数方法。此矩阵中，对角元素的原子性质的计算要折合成C元素的性质，非对角元素与成键原子对的键级有关。Dragon能提供基于Burden矩阵的、8个具有最大特征值及8个最小特征值的分子参数。已经证实Burden特征值参数可以反映分子的有关结构特征，并用于分子结构相似性或差异性的区分[15]。

第2个分子参数SME属于边邻接指数(Edge Adjacency Indices)。边邻接指数本质上为拓扑参数，是基于隐氢分子图的邻接矩阵计算而来，反映分子图形边相连信息。矩阵的维数与取决于非氢原子形成的化学键数目，当化学键相邻或不相邻时，矩阵中对应的元素为1或0。当邻接矩阵的零对角元素用具体的化学键性质取代时，形成的增广边邻接矩阵可描述具体的化学结构信息。用于加权处理的化学键性质包括与共振积分有关的参数(Parameters Related to the Resonance Integral)、邻边数目(Edge Degree)、键级(Conventional Bond Order)、偶极矩(Dipole Moment)。分子参数SME能反映分子的几何特征等结构信息。

第3个分子参数MVV为采用van der Waals体积加权处理的3D-MoRSE参数[15]。这类参数表示在各个方向的电子散射强度。分子参数MVV能反映分子三维结构等信息，也具有区分同系列分子结构的相似性或差异的能力[9,15]。

第4个分子参数NSA为脂肪族化合物中含仲酰胺基团数目；第5个分子参数FNN为键的拓扑距离为5出现N-N键的频率；第6个分子参数FCO为键的拓扑距离为6出现C-O键的频率；最后1个分子参数FCS为键的拓扑距离为7出现C-S键的频率。

表1 基于二元Logical回归分类表Table 1 Classification table with binary logical regression analysis

表2 模型所用分子参数特征Table 2 Characteristics of descriptors appearing in model

(续表2)

ModelDescriptorBhS.E.iWalsjdfkSig.lExp(B)mModel6fLEB126.42754.8105.32110.0218.067E54SME1.8900.58510.43110.0016.621MVV-0.4950.13513.49210.0000.609NSA0.7790.2658.68110.0032.180FCO-0.3350.05044.65210.0000.715FCS-0.3430.1198.31310.0040.709Constant-949.942407.2765.44010.0200.000Model7gLEB124.84957.9594.64010.0311.664E54SME2.3710.63214.06710.00010.704MVV-0.7060.15919.71910.0000.494NSA0.8340.2759.19510.0022.302FNN-0.4470.1489.07910.0030.640FCO-0.3740.05644.82610.0000.688FCS-0.4170.12610.91710.0010.659Constant-940.762430.7594.77010.0290.000

aVariable(s) entered on Model 1:FCO.bVariable(s) entered on Model 2:MVV.cVariable(s) entered on Model 3:FCS.dariable(s) entered on Model 4:LEB.eVariable(s) entered on Model 5:SME.fVariable(s) entered on Model 6:NSA.gVariable(s) entered on Model 7:FNN.hB is Logistic regression coefficients;iS.E.is the standard error;jWals is the statistic;kdf is degree of freedom;lsig.is p-value test;andmExp(B) is the odds ration.

表3 分子参数符号、定义与分类Table 3 The symbol,definition and class of molecular descriptors selected

2.3 分类模型讨论

在SELEX筛选中，开始时核酸适配体亲和性较小；随着筛选轮数的增加，亲和性开始迅速增大；而随着SELEX筛选轮数继续增加时，亲和性增大变得缓慢。亲和性达到一个饱和点之后，随着筛选轮数的进一步增加，亲和性增加并不明显。因此，本研究将从第3轮筛选获得的候选核酸序列类标签值设置为“1”(代表低亲和性和低特异性候选核酸适配)，而把从第11轮筛选所得序列类标签值设置为“2”(表示高亲和性和高特异性候选核酸适配)是合理的。

本文模型LOO交叉验证准确率高(76%)；测试集与训练集的准确率接近；且SVC模型最佳正则化常数C(=357.7)及Gaussian函数参数γ(=0.011)相对较小，模型结构简单。因此本文所建模型不存在过拟合。

将SVC模型对第5、7、9轮的序列进行了预测，含类标“2”的样本分别是41、61、64条。通常在SELEX实验过程中，当靶标分子数目远小于核酸的序列数目时，候选核酸适配体序列与靶标结合的亲和性随筛选轮数的增加而增大。因此本文预测的第3、5、7、9、11轮含类标值为“2”(高亲和性、高专一性)的样本分率逐渐增大，分别是0.32、0.41、0.61、0.64、0.87，这与核酸适配体进行规律一致。

基于模型预测，表4中的4条序列应该为高亲和性、高特异性的核酸适配体。而实验结果显示，表4显示核酸适配体序列Xq1、Xq2、Xq3、Xq4与标靶PC-3M-1E8及PC-3M-2B4结合解离常数为2～47 mmol[6]，预测结果非常满意。

表4 基于分类模型验证的4条候选序列Table 4 Four candidate aptamers validated with classification model

目前流行的核酸序列分类(如序列同源性比对)是基于序列的相似性。我们应用DNAStar软件包中MegAlign软件对上述数据集第3、11轮核酸序列进行进化树分析，比对方法为ClustalW。根据SELEX筛选特点，随着筛选轮数的增加，产品中核酸序列的相似性也会增加，即进化树中相似性最大的序列应该主要由第11轮的序列组成。进化树分析结果显示，相似性最大的84条序列中，仅有50条序列来自第11轮筛选所得的产品，分类准确率接近60%。此结果低于本文SVC模型预测总的准确率(87.0%)。因此，目前基于序列相似性进行比对的算法，不能实现本文核酸序列亲和性高、低两类的区分。其原因是基于序列相似性的比对分析，所利用的信息是序列中碱基的排列顺序与位置，该方法的实质是利用了序列的一维拓扑信息；而本文的分类方法是基于序列的三维结构信息。因此本文分类方法与目前流行的分类方法在基于的序列结构信息方面有着本质区别。

QSAR/QSPR研究的主要目的是建立起分子结构与性能之间的关系式。由于核酸适配体序列结构与其翻译后的氨基酸序列相对应，即适配体序列化学结构与相应的氨基酸序列化学结构相关联。因此从氨基酸计算得到的7个分子参数隐含着核酸适配体的结构信息，由这些分子参数建立的二元回归模型可以用于核酸适配体序列亲和性、专一性的预测。

3 结论

本文成功建立了以前列腺癌PC-3M-1E8细胞为标靶的核酸适配体亲和性的SVC模型，模型对训练集(150条核酸序列)、测试集(50条核酸序列)的预测准确度分别为87.3%和86%。SVC模型对SELEX筛选实验的第3、5、7、9、11轮的核酸适配体进行了预测，其中高亲和性与高特异性核酸适配体的分率分别是0.23、0.41、0.61、0.64、0.87，预测结果符合SELEX筛选的适配体进行规律。因此，间接从氨基酸序列计算分子参数、建立模型，预测以前列腺癌PC-3M-1E8细胞为标靶的核酸适配体亲和性是可行的。