Donor剪接位点识别算法与新基因的寻找

Donor剪接位点识别算法与新基因的寻找

-WAM模型等等

华中科技大学生命科学与技术学院生信息基地1101班2014/6/1

摘要：Donor剪接位点的识别，在国际上已经有很多相关的算法，如WMM、WAM等模型。本文主要是应用C语言，将WMM算法独立编写出来，以进一步掌握对Donor剪接位点的认识，并学会利用现有知识初步开发生物信息学相关算法。而对于我们要寻找的新基因序列，则是指在数据库中已经存在，但在蛋白质水平上还没有完全匹配的基因序列，或者是在蛋白质水平上也有完全匹配的但却来自于另一个物种的基因序列。本文同时也介绍寻找新基因的一般流程，并给出实例。

关键词：Donor剪接位点、WAM模型、新基因、BLAST

1 引言

当人类基因组研究进入一个系统测序阶段时，基因组的研究热点开始转向揭示基因信息结构的复杂性与遗传语言的根本规律。其中，基因预测算法的研究也成为对基因组序列进行统计分析的重要目标。所谓基因预测，一般是指预测DNA序列中编码蛋白质的部分，它是在对DNA序列编码潜能（coding potentials）提出某种模式（pattern）描述的基础上，对一未知的DNA序列上完整的基因结构进行注释。基因预测的最终目标是预测完整的基因结构，正确地识别出一个基因的所有外显子及其边界。

预测DNA序列中蛋白质编码区域的方法主要是基于特征信号的识别。真核基因外显子（编码区域）具有一些特别的序列信号，如内部的外显子被剪切接受体位点（Accepted sites）和供体位点（Donor sites）所界定。根据这些序列特征信号确定外显子的边界，从而达到识别编码区域的目的。最初基因分析方法是进行简单的核苷酸统计，而后加上剪切保守位点的检测。以后采用了人工神经网络（ANN）、隐马尔柯夫模型（HMM）等先进的信息处理和分析技术，并与同源序列搜索结合起来，通过与已知基因序列或者EST序列的比较，提高基因识别的准确率。常见的编码区分析工具通常将多种技术组合起来，给出对编码区的综合判别。但鉴于这样的算法过于复杂，本文选择的是最为简单的模型，即WMM模型。

而新基因的寻找，事实上是为了发现基因的功能，而非单纯的寻找某个基因。它是基于现有的公共数据库，如华盛顿特区的NCBI、英国Hinxton的EBI和日本的DDJB，将感兴趣的基因与这些数据库中的每个序列进行比较，鉴别出相似的序列。搜索结果显示出与最佳匹配序列的的对位排列及匹配计分。如果一个查询序列很容易和一个已知结构、功能或生物活性的数据库序列进行对位排列，则这个查询序列被认为具有相同的结构、功能或生物活性。这些预测的力度依赖于序列间对位排列的质量。作为一条粗略规则，如果待查询的序列和数据库序列中多于一半的碱基或氨基酸在对位排列中相同，这个预测就非常强。随着相似程度的下降，预测的可信度也下降。用于数据库搜索的程序给出了统计评价，以评价对位排列计分。

2 Donor剪接位点识别算法

必须清楚，要想设计一个100%识别编码区域的程序几乎是不可能的。问题是如何提高一个识别算法的敏感性Sn和特异性Sp。Sn 和Sp都应该比较高，若一个算法的测试结果仅仅一个很高，而另一个很低，则该算法是不成功的。例如，假设有一个识别编码区域的算法，它将所有介于AG和GT之间的序列片段都找出来作为识别结果，那么该算法的敏感性Sn将达到100%，但其特异性Sp却近似于0%。因此，对于一个识别算法，往往用敏感性和特异性的平均值作为衡量其准确率的指数，即(Sn+Sp)/2。在一般情况下，调整程序的参数，使得Sn?Sp。

2.1.Donor剪接位点附近碱基信息的提取

在Training Set提供的30个文件中，已经注释出了外显子的存在位置，而Donor剪接位点正好存在于外显子与内含子的分界线处，由此我们可以统计处Donor剪接位点处前3个位置和后六个位置四种碱基出现的概率。具体算法思想如下：

⑴．读取Training Set中的30个*.txt文件(文件读取方法是将所有文件名采取通配符的办法，然后加上文件类型.txt，在Training Set所在的路径中搜索并依次打开该文件)，利用一个string 类的字符串数组Buffer[30]。每一个字符串Buffer[i](i<30)存入一个文件中的所有字符，包括基因序列前的注释性说明。

⑵．判断Buffer[i](i<30)中的基因序列字符。

for (i=0;i<30;i++)

While (Buffer*i+*j+! =’)’) j++;

sequence[i] =Buffer[i].substr (j+1);

于是我们得到一个只存入了基因序列的字符串数组sequence[30]，其中每一个字符串sequence[i](i<30)存入一个文件的基因序列。

⑶．利用Training Set中文件事实上已经提供给我们了Donor剪接位点的确切位置，获取Donor 剪接位点附近的碱基信息。现定义一个char型的二维数据donor[][9],行表示共有9个位置，列则表示这9个位置中每个位置的碱基a、c、g或者t。同时还定义一个二维数组num[150][1000]，行表示文件序号，列存取每个文件中Donor剪接位点的位置。读取Buffer[i]（i<30，表示文件的序号）中的字符，每当遇到一个“,”时，便存下它前面的数字，该数字就是一个已知的Donor剪接位点的确切位置。由此，我们就能得到已知的Donor剪接位点附近9个碱基的信息。其算法思想是：

if (num[i][a]!=0)

for (j=0;j<9;j++)

{c=num[i][a]-3;

donor[b][j]=sequence[i][c+j];

b++

}

其中，b表示存储到第b个Donor剪接位点，i表示第i个文件，j表示第j个位置。

⑷．统计出每个位置a、c、g、t出现的次数。继续定义4个二维数组num_a[2][9]，num_c[2][9]，num_g[2][9]，num_t[2][9]。行表示共有位置，列则表示每个位置该碱基出现的总的次数。其中第一行输出Donor剪接位点的9个位置的编号，第二行则输出每个位置碱基出现的次数。现假设总共有N个Donor剪接位点，故其算法应为：

for(j=0;j<9;j++)

{

for(i=0;i

{ if(donor[i][j]=='a') num_a[1][j]=num_a[1][j]+1;

if(donor[i][j]=='c') num_c[1][j]=num_c[1][j]+1;

if(donor[i][j]=='g') num_g[1][j]=num_g[1][j]+1;

if(donor[i][j]=='t') num_t[1][j]=num_t[1][j]+1;

}

}

⑸．计算各个位置a、c、g、t出现的概率。定义4个一位数组p_a[9]，p_c[9]，p_g[9]，p_t[9]，分别用来存放4个碱基在9个位置出现的概率。基本算法如下：

for(j=0;j<9;j++)

p_x[j]=(double)num_a[1][j]/N;（x∈{a、c、g、t}，N表示Donor剪接位点的总数