DNA序列的统计分析

DNA序列的统计分析

【摘要】模型一统计了20个已知类别的DNA序列碱基的含量的概率分布,根据已知的类别就A,T,C,G的含量作为四个指标,采用判别分析对未知类别的序列给出了较满意的分类。模型二首先统计了已知类别的DNA序列的位置上各碱基出现的概率,发现A,B两类序列结构的不同,体现在密码子各位上的碱基概率分布有明显的差异,以嘌呤和嘧啶碱基为区别构造一个一维随机徘徊函数,从而给出A,B的分类法;接着,再从三个角度来划分碱基,对于每一种分类都构造一个一维随机徘徊函数,根据此函数得出拟和直线,用三条拟和直线的斜率作为分析的指标进行多元判别分析,由此给出A,B的分类法,较模型一分类的正确率明显提高。

一、问题简述与分析

人类基因组计划中DNA全序列图是由四个碱基A,T,G,C按照一定的顺序排成的长约30亿的序列,研究DNA全序列具有什么结构,探讨由这四个字符排成的看似随机的序列中到底隐藏着什么规律,是当代生物信息学最重要的课题之一。

DNA分子中唯一可变的部分是碱基(胸腺嘧啶T,鸟嘌呤G,胞嘧啶C,和腺嘌呤A)序列,人类发现在全序列中有一些是用于编码蛋白质的序列片段,即由这四个字符组成的64种不同的3字符串,其中大多数是用于编码构成蛋白质的20种氨基酸,研究表明,分析DNA序列的结构以及序列的某些片段之间具有的相关性对于理解DNA全序列有十分重要的意义,现提出给以下序列集合进行分类的问题:

1)由20个已知类别的序列中(序号1—10为A类,11—20为B类)提取特征,构造分类方法,并用这些已知类别的序列来衡量你的方法的好坏,然后对另外的20个未知类别的序列(标号21—40)进行分类。

2)对给出的182个DNA序列,用你的分类方法对他们进行分类,并给出分类结果。

研究表明,遗传密码所必要的碱基个数为3个,即密码子是由三个碱基组成,一串前后相依的密码子构成了氨基酸的排列次序,从而形成了具体的蛋白质,显然密码子使用的频率和数量,进一步,碱基出现的频率和数量,特别是排在一起的结构和序列片段的相关性都与研究DNA序列有十分紧密的联系,我们就是要挖掘这些统计特征,寻找出隐藏在这些序列中的规律。

首先,通过分析,我们可以看出给出的A,B两类的20 个样本数据中,四个碱基a, c, g, t 的含量有较明显的区别,因此我们可以通过其在含量方面的区别,以四种碱基的含量为四个指标利用SAS统计软件进行多元判别分析,以此来确定A,B的分类,并进而对其他的序列进行分类。(模型一)

其次,我们进一步判断,发现对a, c, g, t的含量完全相同的两个DNA序列来说,决定其分类的标准就不能再用碱基含量了,此时我们考虑用碱基的结构排列,即a, c, g, t出现在DNA序列中的每一位的顺序:我们先以嘌呤碱基与嘧啶碱基作为分类的标准,并构造一个一维随机徘徊函数,然后用据此得到的拟和直线的斜率来进行判断,但是我们进而发现仅从这一个角度来考虑是不完善的,因此经过研究我们从三个角度来分别构造一维随机徘徊函数,得到三条拟和直线,以这三条直线的斜率为指标再次用SAS统计软件进行多元判别分析,以此来判断A,B的分类。(模型二)

二模型假设与符号设定

1. 假定所给的DNA序列数据为起始密码子之后的第一个数据字符;

2. 每个碱基出现是随机的;

3. ha——一个序列中a的含量,hc——一个序列中c的含量;

4. hg——一个序列中g的含量,ht——一个序列中t的含量;

5. K1——按嘌呤与嘧啶碱基分类拟和的直线的斜率;

6. K2——按氨基与酮基碱基分类拟和的直线的斜率;

7. K3——按强氢键与弱氢键分类拟和的直线的斜率;

8. 其他的符号将在文中另外给出。

三模型一的建立和求解

一) 样本的统计分析

从含量的角度考虑,对于给出的20个已知类别的样本数据,我们利用MATLAB绘制出a, c, g, t 的分布图如下:(其中实线表示A类,虚线表示B类)

这里采用MATLAB的图形函数plot做图求解,其做图格式为:plot(x,a1,x,a2,'--')。其中X是横坐标,取1到10,a1与a2分别为A类与B类中的碱基含量,'--'代表线型是虚线。

DNA序列的统计分析

DNA序列的统计分析

a 的分布 c 的分布

DNA序列的统计分析

DNA序列的统计分析

g 的分布 t 的分布

图1 A 与B 两类a, c, g, t 的分布图

由上图可以看出,a, c, g , t 的含量明显不同,特别是g , t 的含量差别很大,因此我们可以根据a, c, g , t 的含量来区分A ,B 两类。于是我们将已知的20种序列和未知的20种序列的a, c, g , t 的含量计算出来并列表如下:

表1 A ,B 两类a ,c ,g ,t 的含量表

DNA序列的统计分析

DNA序列的统计分析

表2 未知分类的20 种序列a, c, g , t 的含量表

DNA序列的统计分析

在这里,衡量分类的标准有四个指标:ha , hc , hg , ht。这四种指标分别表示一个序列中的四种碱基a, c, g , t的含量。

接着我们用判别分析来进行处理。判别分析是多元统计分析中判别样品所属类型的一种重要方法,是根据多种因素(指标)对事物的影响,进行判别分类的统计方法,又称为鉴别分析。其具体的方法是:已知了研究对象分成若干个类型(或组别)并已取得各种类型的一批已知样品的观测数据,在此基础上总结出分类的规律性或判别信息,再根据某些准则建立判别式,然后对未知类型的样品进行判别分类。此处,我们正好已知了20个样本的类别,而需要对另外的20个样本进行分类,刚好符合判别分析的原理,因此我们采用判别分析来进

行处理。此处,我们用专门的统计软件SAS 来实现判别分析的过程。

二) 具体的程序操作

SAS 的具体操作程序见附录一。 程序运行结果如下:

DNA序列的统计分析

以上是一些基本情况

DNA序列的统计分析

上表为各组的基本情况,并列出了各组的先验概率值。因为指定了Prior Probability ,所以各组的先验概率按实际数据中各组比例计算。

DNA序列的统计分析

DNA序列的统计分析

上面为各组均值间广义距离平方的公式,即)()()(12j i j i i j X X S X X X D -'-=-(其

中S 为合并协方差阵)。

DNA序列的统计分析

DNA序列的统计分析

上面即线性判别函数的公式,给出了到第j类的线性评判别函数的常数项和自变量的系数向量的公式。下面具体给出了各类的线性判别函数的各常数项与系数值:

DNA序列的统计分析

下面为判别分析对训练数据集(Calibration Data)用线性判别函数进行回判的概况,先给出了广义平方距离函数的公式和每个已知类别的样本属于各类的后验概率的公式,然后

DNA序列的统计分析

DNA序列的统计分析

DNA序列的统计分析

DNA序列的统计分析

下面是对训练数据集进行交叉核实回判的情况。交叉核实的想法是:为了判断观测样本i 的判别正确与否,用删除第i 个观测样本后的训练数据集算出判别函数,然后用此判别函数来判别第i 观测样本。对每一个观测样本都进行这样的判别。在结果中,先给出了广义平方距离函数和后验概率公式,所以公式中用了j X X )(表示除了 X 所在观测样本后的第j 组的均值,用)(X COV 表示除了 X 所在观测样本后得到的合并协方差阵估计。

DNA序列的统计分析

DNA序列的统计分析

DNA序列的统计分析

DNA序列的统计分析

以下是用根据训练数据集得出的判别函数对检验数据集(Test Data)进行判别的结果。

DNA序列的统计分析

DNA序列的统计分析

DNA序列的统计分析

DNA序列的统计分析

三) 主要统计结果分析:

由表Linear Discriminant Function for type可以看出,线性判别函数为:

F1= -2398664+4797930X1+4789763X2+4805558X3+4797035X4;(A类)

F2= -2397786+4797053X1+4788916X2+4804612X3+4796190X4;(B类)

判别函数用于将未知的20种序列分类,即将未知的序列的四个指标值带入两个判别函数,哪一个判别函数的值最大,就估计属于哪一类。

由表Classification Results using Linear Discriminant Function Posterior Probability of Membership in type的结果可以看出对未知的20 种序列的分类如下:

A类:22 23 25 27 28 29 34 35 36 37;

B类:21 24 26 30 31 32 33 38 39 40;

由用线性判别函数进行回判的结果中的 Number of Observations and Percent Classified into type 表可以看出,其中有一种属于A类的序列被错判为B类,因此回代的符合率=19/20=95%,这个比率是比较高的,说明回代的效果比较理想;由交叉核实回判结果中的 Function Posterior Probability of Membership in type 表可以看出:第4序列由A类错判为B类,而第17序列由B类错判为A类,从而总的错判率=2/20=10% ,该错判率总的来看是比较低的,也说明利用此方法判断得出的结果是比较好的。但是我们认为仅仅如此分类是不完善的,因为就分类的标准而言,我们就可以看出仅从“不同的序列中碱基含量不同”这个角度来考虑,以碱基的含量作为判断的唯一标准太过于片面,因此我们有必要寻求更合理的分类方法来完善对于序列的分类以得到更准确的结果。

四模型二的建立和求解

一)碱基的概率统计分析

在模型一中,采用碱基A,C,G,T在DNA序列中的含量作为序列的特征来进行分类,这固然有一定的生物意义,并且也在DNA序列的分类中获得了较满意的结果,但是仅仅使用该特征作为分类的唯一标准并没有充分地体现碱基排列的信息量,仅仅只是考虑了碱基的含量并没有考虑到碱基在序列中的排列情况。例如,序列(AGCT)与序列(GCAT)有着相同的碱基含量,但我们不能就此认为他们是没有区别的。因此,直接从序列本身的碱基排列顺序来考察

序列就成为一种更合适的提取特征的方式。由生物学的知识可知,DNA链上开始密码子的后几十位和终止密码子的前几十位中,每一位出现碱基a ,c ,g ,t的概率呈现明显的特征。因此,我们有必要就a ,c ,g ,t对A,B两类的样本数据求出密码子各位点上的碱基出现的概率分布。对已知的20 种序列各位上的碱基出现频率,再次利用MATLAB的图形函数plot绘图如下(其中,实线表示A类,虚线表示B类):

DNA序列的统计分析

a

DNA序列的统计分析

c

DNA序列的统计分析

g

DNA序列的统计分析

t

图2 各点碱基出现概率分布

二)构造 A,B 的分类法

(1)构造一维随机徘徊函数

DNA序列是由a,g,c,t这四个字母组成的序列。1992年Peng等的工作已经揭示了DNA序列中存在长程相关性。由图2的分析可知,A,B两类每一位出现a,g,c,t的概率有明显的差异,正因为如此,我们可以再从“不同序列中碱基位置不同”这个角度来考虑。而A,G 同为嘌呤碱基,C,T 同为嘧啶碱基,我们可以利用序列先后出现嘌呤或嘧啶碱基概率的不同将DNA 序列转换表示为一个一维的随机徘徊函数。

方法:从第一个碱基(即第一个字母)算起,若是嘌呤碱基(即A或G)则负走一步,若是嘧啶碱基(即C或T)则正走一步。记n步后的净位移为f n,(n=1,2,3,…,L)L为序列的长度。在长度L的窗口里计算净位移f n。

下图是从前20 个样本数据中抽出序列号为3、7(A类),15、18(B类)的四个样本数据利用MATLAB的图形函数plot 绘图进行的定性的分析。

DNA序列的统计分析

DNA序列的统计分析

A 3 A 7

DNA序列的统计分析

DNA序列的统计分析

B 15 B 18

图3

由图3可以看出,经过一维随机徘徊函数转换以后,A类与B类存在明显的差异,A类序列呈下降趋势,而B类序列呈上升的趋势,经过定性分析,我们可以利用这L个离散点作出其拟和直线,把直线的斜率K作为衡量A,B两类的标准。

我们记n=0,f0 =0,把(0,0)作为基准点。则:

0, n=0;

f

= f n-1–1,序列的第n位为a或g;

n

f

+1,序列的第n位为 c或t;

n-1

(2)由直线回归结果进行分类

我们在XOY坐标系中以n 值作为x向量坐标,f n 值作为y向量坐标,设直线y i= Kx i+b 为序列号为i 的序列的拟和直线。在此,我们用SAS中的回归分析来拟和直线,进而得到斜率K值。以A类第3序列为例,在SAS 中做回归分析,结果如下:

DNA序列的统计分析

DNA序列的统计分析

DNA序列的统计分析

由以上的结果可知,第3序列的拟和直线为:y = -0.4432x + 2.05825,即斜率K = -0.4432。同理,可求出已知类别的20种序列中其他的序列的拟和直线的斜率以及未知类别的20种序列的拟和直线的斜率。经过运算,各序列拟和直线的斜率如下表所示:(精度:0.00001)

表3

DNA序列的统计分析

由上表可知,对需要判断的DNA序列样本,我们的步骤是:通过一维随机徘徊函数转换,采用线性回归分析方法计算出拟和直线的斜率K值,然后根据A,B两类的明显的K值特征来进行判断,即K值小于-0.11469的一定属于A类,而K值大于-0.01669的则一定属于B类。因此根据此方法我们可以判断出,分类结果如下:

A类:27 29 34 ; B类:21 24 26 28 30 31 32 33 36 38 39 40 ;

(3)对无法判断的序列进行分类

但是我们可以看到,对K值界于 -0.11469和 -0.01669之间的样本,我们就没有办法判断到底属于哪一类,此处就有22,23,25,35,37五个样本无法判断。于是我们利用系统聚类法将剩余的样本数据聚类成两大类,再利用该两大类的K值来判断A、B类。由于K值是一维随机的,因此可以利用最短距离法将界于 -0.11469和 -0.01669之间的样本聚类成两类。

最短聚类法的聚类步骤如下:

(1)选取聚类统计量。在这里我们选取了绝对距离。

(2)计算每个样本之间的两两距离,并记在分类距离表中,记为D(0),这是第0步的表,每

个样本为一类,D pq表示每两个类之间的距离(p、q为类号)。

(3)选择D(0)中的最短距离,设为D pq,则将p、q两类合并成一个新类,记为r类, r={p、q},

表示由p类和q类所组成。

(4)计算新类D r与其他类之间的距离,定义D rk = min{D pk, D qk}, (1)

(5)作D(1)表,将D(0)中的第p、q行p、q列删去,加上第r行、r列。第r行、r列与其他类的

距离按(1)式判断后记上,这样得到一个经过一次聚类后的新的分类距离对称表。需要注意的是D r类是由哪两类聚类得到的,应该在D(1)表下给予说明。

(6)对D(1)按(3)、(4)、(5)步重复类似D(0)的聚类工作,可以得到D(2),这是经过二步聚类

得到的一个新的分类距离对称表。

(7)重复以上的步骤,直到最后只剩下两类为止

这样我们就可以对剩下的三类按最短距离法判别是属于A类还是属于B类。此处,我们判别的结果是:23,25,35,37属于A类,22属于B类。

三)进一步的改善法

上面我们是利用一维随机徘徊函数来计算K值,只是从按嘌呤和嘧啶碱基分类这一个角度来考虑的,很显然这是不全面的。因此,为了更好的对序列进行分类,我们应该从更多的方面来充分的考虑。实际上,从纯化学的角度,我们就可以从两个角度来对碱基进行两类的划分:(1)如上的分类,按双环或单环结构,可分为:嘌呤碱基(A或G)与嘧啶碱基(C 或T);(2)按环中对应位置上是否存在氨基或酮基。可分为:氨基碱基(A或C)与酮基碱基(G或T):然后,从生物学的角度,在双螺旋结构中,按碱基对形成氢键的数目或强弱,碱基又可分为:强氢键碱基(G或C)与弱氢键碱基(A或T)。这样,我们就把化学和生物的信息都考虑进去了,对第二种和第三种分类方法如第一种分类方法一样,分别构造一个随机徘徊函数进行分析,我们发现,在每一种分类方法中,都存在一种线性关系,因此,我们考虑,对每一个序列分别按三种方法构造一维随机徘徊函数,从而得到三条拟和直线,以这三条直线的斜率K1,K2,K3作为三个指标再进行多元判别分析,以此得到未知序列的分类。(具体的程序过程见附录二。)

(1)程序运行结果:

注:程序运行结果中各表格所代表的具体含义在模型一中已有较详细的说明,此处就不再重复解释了。只将其列出如下:

DNA序列的统计分析

DNA序列的统计分析

Pooled Covariance Matrix Information Natural Log of the Covariance Determinant of the Matrix Rank Covariance Matrix

3 -11.10919

DNA序列的统计分析

DNA序列的统计分析

DNA序列的统计分析

DNA序列的统计分析

DNA序列的统计分析

DNA序列的统计分析

DNA序列的统计分析

DNA序列的统计分析

DNA序列的统计分析

DNA序列的统计分析

DNA序列的统计分析

DNA序列的统计分析

DNA序列的统计分析

DNA序列的统计分析

DNA序列的统计分析

DNA序列的统计分析

(2)主要结果分析:

由Linear Discriminant Function for type 表可以知道,判别函数如下:

F1= -3.80071 –20.00952K1+14.39466K2 – 2.86845K3;(A类)

F2= -5.01463 + 5.79269K1 + 2.66024K2 + 13.58387K3;(B类)

判别函数用于将未知的20种序列分类,即将未知的序列的三个K值分别带入两个判别函数,哪一个判别函数的值最大,就估计其属于哪一类。

由表 Classification Results using Linear Discriminant Function Posterior Probability of Membership in type的结果可以看出对未知的20 种序列的分类如下:

A类:22 23 25 27 29 34 35 36 37;

B类:21 24 26 28 30 31 32 33 38 39 40;

在SAS软件中,用判别分析进行数据处理时,一般只要正确率达到 85% 以上就说明判别的效果很好了。在此模型中,与模型一相同,正确率都已经在 85% 以上,达到了95%,由此可知判别效果已经很好了。而由交叉核实回判结果中的Posterior Probability of Membership in type表可以看出与模型一相比较,误判率又有了下降,而且由对未知类别的序列用线性判别函数进行分类的表 Posterior Probability of Membership in type 中也可以看出对于各种序列的判别分类,其后验概率也比模型一要好,这些都说明了分类的效果有了明显的提高。

四)对182个序列的分类结果

由上可知,用以上的方法对DNA序列进行分类,效果比较理想,于是我们以此作为判断分类的方法,对给出的另外182个序列进行分类,(具体程序见附录三),结果如下:

表4

DNA序列的统计分析

DNA序列的统计分析

五模型的优缺点

相关推荐
相关主题
热门推荐