DNA序列的统计分析.docx
《DNA序列的统计分析.docx》由会员分享，可在线阅读，更多相关《DNA序列的统计分析.docx（43页珍藏版）》请在冰豆网上搜索。
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类中的碱基含量，'--'代表线型是虚线。
a的分布c的分布
g的分布t的分布
图1A与B两类a,c,g,t的分布图
由上图可以看出，a,c,g,t的含量明显不同，特别是g,t的含量差别很大，因此我们可以根据a,c,g,t的含量来区分A，B两类。
于是我们将已知的20种序列和未知的20种序列的a,c,g,t的含量计算出来并列表如下：
表1A，B两类a，c，g，t的含量表
数据号
含量
a%
c%
g%
t%
1
0.297
0.171
0.396
0.135
2
0.272
0.155
0.418
0.155
3
0.275
0.220
0.440
0.064
4
0.426
0.111
0.176
0.287
5
0.243
0.224
0.421
0.112
6
0.349
0.132
0.396
0.123
7
0.352
0.105
0.352
0.190
8
0.288
0.154
0.365
0.192
9
0.204
0.214
0.427
0.155
10
0.198
0.267
0.396
0.139
11
0.360
0.050
0.090
0.500
12
0.333
0.020
0.152
0.495
13
0.245
0.091
0.112
0.551
14
0.330
0.082
0.113
0.474
15
0.302
0
0.063
0.635
16
0.379
0.095
0.094
0.432
17
0.383
0.255
0.128
0.234
18
0.301
0.118
0.097
0.484
19
0.207
0.163
0.076
0.554
20
0.198
0.165
0.066
0.571
表2未知分类的20种序列a,c,g,t的含量表
数据号
含量
a%
c%
g%
t%
21
0．274
0．195
0．168
0．363
22
0．289
0．240
0．250
0．221
23
0．176
0．255
0．382
0．186
24
0．209
0．191
0．191
0．409
25
0．248
0．229
0．305
0．219
26
0．219
0．211
0．184
0．386
27
0．231
0．202
0．337
0．231
28
0．256
0．145
0．154
0．445
29
0．149
0．218
0．446
0．188
30
0．290
0．243
0．215
0．252
31
0．241
0．179
0．223
0．357
32
0．174
0．229
0．266
0．330
33
0．270
0．189
0．207
0．333
34
0．235
0．235
0．363
0．167
35
0．243
0．214
0．340
0．204
36
0．229
0．305
0．257
0．210
37
0．214
0．252
0．330
0．204
38
0．222
0．171
0．171
0．436
39
0．274
0．283
0．208
0．236
40
0．198
0．198
0．172
0．431
在这里，衡量分类的标准有四个指标：
ha,hc,hg,ht。
这四种指标分别表示一个序列中的四种碱基a,c,g,t的含量。
接着我们用判别分析来进行处理。
判别分析是多元统计分析中判别样品所属类型的一种重要方法，是根据多种因素（指标）对事物的影响，进行判别分类的统计方法，又称为鉴别分析。
其具体的方法是：
已知了研究对象分成若干个类型（或组别）并已取得各种类型的一批已知样品的观测数据，在此基础上总结出分类的规律性或判别信息，再根据某些准则建立判别式，然后对未知类型的样品进行判别分类。
此处，我们正好已知了20个样本的类别，而需要对另外的20个样本进行分类，刚好符合判别分析的原理，因此我们采用判别分析来进行处理。
此处，我们用专门的统计软件SAS来实现判别分析的过程。
二）具体的程序操作
SAS的具体操作程序见附录一。
程序运行结果如下：
TheDISCRIMProcedure
Observations20DFTotal19
Variables4DFWithinClasses18
Classes2DFBetweenClasses1
以上是一些基本情况
ClassLevelInformation
VariablePrior
typeNameFrequencyWeightProportionProbability
1_11010.00000.5000000.500000
2_21010.00000.5000000.500000
上表为各组的基本情况，并列出了各组的先验概率值。
因为指定了PriorProbability，所以各组的先验概率按实际数据中各组比例计算。
PooledCovarianceMatrixInformation
NaturalLogofthe
CovarianceDeterminantofthe
MatrixRankCovarianceMatrix
4-32.32840
TheDISCRIMProcedure
PairwiseGeneralizedSquaredDistancesBetweenGroups
____
D2（i|j）=（Xi–Xj）'COV-1（Xi–Xj）
上面为各组均值间广义距离平方的公式，即
（其中S为合并协方差阵）。
GeneralizedSquaredDistancetotype
Fromtype12
1027.76546
227.765460
LinearDiscriminantFunction
（线性判别函数）
___
Constant=-.5X'jCOV-1XjCoefficientVector=COV-1Xj
上面即线性判别函数的公式，给出了到第j类的线性评判别函数的常数项和自变量的系数向量的公式。
下面具体给出了各类的线性判别函数的各常数项与系数值：
LinearDiscriminantFunctionfortype
VariableLabel12
（变量）（标签）
Constant-2398664-2397786
（常数）
x1a的含量47979304797053
x2c的含量47897634788916
x3g的含量48055584804612
x4t的含量47970354796190
下面为判别分析对训练数据集（CalibrationData）用线性判别函数进行回判的概况，先给出了广义平方距离函数的公式和每个已知类别的样本属于各类的后验概率的公式，然后是每一类被判入各类的个数和百分比：
TheDISCRIMProcedure
ClassificationSummaryforCalibrationData:
WORK.DNA1
ResubstitutionSummaryusingLinearDiscriminantFunction
（用线性判别函数回代结果摘要）
GeneralizedSquaredDistanceFunction
（广义平方距离函数）
__
D2j（X）=（X-Xj）'COV-1（X-Xj）
PosteriorProbabilityofMembershipinEachtype
（各类的后验概率公式）
Pr（j|X）=exp（-.5D2j（X））/SUMexp（-.5D2k（X））
k
NumberofObservationsandPercentClassifiedintotype
（样本个数与分类比率）
Fromtype12Total
19110
90.0010.00100.00
201010
0.00100.00100.00
Total91120
45.0055.00100.00
Priors0.50.5
下面为各类的错判率：
ErrorCountEstimatesfortype
12Total
Rate0.10000.00000.0500
Priors0.50000.5000
下面是对训练数据集进行交叉核实回判的情况。
交叉核实的想法是：
为了判断观测样本i的判别正确与否，用删除第i个观测样本后的训练数据集算出判别函数，然后用此判别函数来判别第i观测样本。
对每一个观测样本都进行这样的判别。
在结果中，先给出了广义平方距离函数和后验概率公式，所以公式中用了
表示除了X所在观测样本后的第j组的均值，用
表示除了X所在观测样本后得到的合并协方差阵估计。
TheDISCRIMProcedure
ClassificationResultsforCalibrationData:
WORK.DNA1
Cross-validationResultsusingLinearDiscriminantFunction
GeneralizedSquaredDistanceFunction
__
D2j（X）=（X-X（x）j）'COV-1（x）（X-X（x）j）
PosteriorProbabilityofMembershipinEachtype
Pr（j|X）=exp（-.5D2j（X））/SUMexp（-.5D2k（X））
k
下面是进行交叉核实回判的具体结果：
TheDISCRIMProcedure
ClassificationResultsforCalibrationData:
WORK.DNA1
Cross-validationResultsusingLinearDiscriminantFunction
（用线性判别函数交互证实回代结果）
PosteriorProbabilityofMembershipintype
FromClassified
Obstypeintotype12
412*0.00001.0000
1721*0.99870.0013
*Misclassifiedobservation
接着是对各类进行交叉核实回判的概况：
NumberofObservationsandPercentClassifiedintotype
Fromtype12Total
19110
90.0010.00100.00
21910
10.0090.00100.00
Total101020
50.0050.00100.00
Priors0.50.5
下面是用交叉核实计算的各类的错判率：
ErrorCountEstimatesfortype
12Total
Rate0.10000.10000.1000
Priors0.50000.50
以下是用根据训练数据集得出的判别函数对检验数据集（TestData）进行判别的结果。
先给出了广义平方距离函数的公式和每个未知分类的样本属于各类的后验概率的公式：
TheDISCRIMProcedure
ClassificationResultsforTestData:
WORK.DNA2
ClassificationResultsusingLinearDiscriminantFunction
GeneralizedSquaredDistanceFunction
__
D2j（X）=（X-Xj）'COV-1（X-Xj）
PosteriorProbabilityofMembershipinEachtype
Pr（j|X）=exp（-.5D2j（X））/SUMexp（-.5D2k（X））
k
接着下表就是具体的分类结果：
TheDISCRIMProcedure
ClassificationResultsforTestData:
WORK.DNA2
ClassificationResultsusingLinearDiscriminantFunction
（用线性判别函数分类的结果）
PosteriorProbabilityofMembershipintype
Classified
Obsintotype12
120.00050.9995
210.78650.2135
311.00000.0000
420.00060.9994
510.99820.0018
620.00050.9995
710.99990.0001
811.00000.0000
911.00000.0000
1020.10050.8995
1120.04280.9572
1220.15470.8453
1320.00990.9901
1411.00000.0000
1510.99990.0001
1610.73560.2644
1710.99910.0009
1820.00010.9999
1920.07590.9241
2020.00001.0000
然后是各待判样本被判入各类的个数和百分比：
NumberofObservationsandPercentClassifiedintotype
12Total
Total101020
50.0050.00100.00
Priors0.50.5
三）主要统计结果分析：
由表LinearDiscriminantFunctionfortype可以看出，线性判别函数为：
F1=-2398664+4797930X1+4789763X2+4805558X3+4797035X4；（A类）
F2=-2397786+4797053X1+4788916X2+4804612X3+4796190X4；（B类）
判别函数用于将未知的20种序列分类，即将未知的序列的四个指标值带入两个判别函数，哪一个判别函数的值最大，就估计属于哪一类。
由表ClassificationResultsusingLinearDiscriminantFunctionPosteriorProbabilityofMembershipintype的结果可以看出对未知的20种序列的分类如下：
A类：
22232527282934353637；
B类：
21242630313233383940；
由用线性判别函数进行回判的结果中的NumberofObservationsandPercentClassifiedintotype表可以看出，其中有一种属于A类的序列被错判为B类，因此回代的符合率=19/20=95%，这个比率是比较高的，说明回代的效果比较理想；由交叉核实回判结果中的FunctionPosteriorProbabilityofMembershipintype表可以看出:
第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类）：
a
c
g
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步后的净位移为fn,（n=1,2,3,…，L）L为序列的长度。
在长度L的窗口里计算净位移fn。
下图是从前20个样本数据中抽出序列号为3、7（A类），15、18（B类）的四个样本数据利用MATLAB的图形函数plot绘图进行的定性的分析。
升级会员 