DNA序列分类模型.doc

DNA序列分类模型重庆市数学建模竞赛一等奖王 勇, 莫志锋, 秦力顼(1999级自动化学院)摘要 本文根据题中所给两个已知类别的DNA序列进行结构特征分析，从中提取信息和构造分类模型，对未知类别的DNA序列进行分类。我们构造了三个分类模型，它们分别是：特征密码子概率分布判别模型、图论最小支撑树模型和向量空间直观判别模型。后两种分类结果几乎一致，判别率在90%左右，误判率控制在（0.05-0.1）范围。问题一结果为：模型一的结果：A类有7个：22，23，27，29，34，35，37；B类有10个：21，24，26，28，30，31，32，33，38，40；不能判断的有3个：25，36，39；模型三的结果：A类有10种：22，23，25，27，29，34，35，36，37，39；B类有10种：21，24，26，28，30，31，32，33，38，40；问题二结果为：模型二的结果：A类有108个，B类有74个。具体情况见文中答案。模型三的结果：A类有120个，B类有62个。具体情况见文中答案。 我们还对三种分类方法进行了类比，认为模型二、三方法新颖独特，结果稳定，它们是一种较好的分类方法。并且对各种计算结果进行误差分析和检验等工作。一、 问题的重述 本问题为一个DNA序列分类问题。假定已知两组人工已分类的DNA序列（20个已知类别的人工制造的序列），其中序列标号110 为A类，11-20为B类。要求我们从已经分类了的DNA序列片段中提取共同特征构造分类方法，并评价所用分类方法的好坏，从而构造或选择一种较好的分类方法。测试对象是20个未标明类别的人工序列（标号2140）和182个自然DNA序列。二、模型的假设及符号说明1、名词解释：碱基： 在生物学中，用A，T，C，G四个字符代表组成DNA序列的四种碱基；密码子：在遗传学中每三个碱基的组合被称为一个密码子，可以编码一个氨基酸，共有64个，还可以由密码子组成20个氨基酸。特征密码子：若类中某个密码子出现的概率pA与类中同一个密码子出现的概率pB之差的绝对值大于0.05，并且，pA pB，定义该密码子为类特征密码子。同理，可定义B类特征密码子。2、模型假设： 所给的DNA序列片段中没有断句和标点符号。对密码子的计数没有固定的起点 和终点。 假定各密码子近似呈现正态统计规律。 DNA序列的某些片段之间具有一定的相关性。 DNA序列不构成环状，划分密码子时，碱基剩余个数小于3时舍去不考虑。3、符号说明：- 表示第k种密码子在第i个序列片段中出现的频率；- 类特征密码子；- 类特征密码子；-最适A类密码子的概率；-最适B类密码子的概率；-误判率；-空间向量夹角余弦之差；三、模型的分析和建立1、A、B两类的DNA序列特征提取（1） 碱基含量分析对A，B两类中的碱基a，t，c，g出现的频率进行了统计，其分布状况如下： A类：t，c 的总含量在25 41%左右，a，g 的总含量在59 75%左右；B类：t，c 的总含量在50 74%左右，a，g 的总含量在26 49%左右；以上的碱基含量分布作为已知DNA序列的第一特征。由此，我们可以根据碱基含量对未知类型的DNA序列片段进行粗略的分类：问题（1）：A类:29, 34;B类:21, 24, 26, 28, 31, 32, 33, 36, 38, 39, 40;不能辨别类：22，23，25，27，30，35，37；在此基础上再对A、B两类进行第二、第三特征提取。针对不能判别类型利用其他特征和分类方法进行鉴别。（2）密码子的概率分布状况分析由遗传学中的知识可知,任何基因都是由64种密码子构成,不同的DNA序列片段中各密码子出现的概率不同。以密码子出现的概率不同为出发点，我们分别对A、B两类中关于DNA序列中的64种密码子的平均频率进行了统计，并且分析了A、B两类的特征密码子。如下统计表1：表1：64种密码子的平均概率分布表序号密码子A(pA)B(pB)|pA-pB|序号密码子A(pA)B(pB)|pA-pB|1AAA0.0450.0450.000433CAA0.01380.00650.00732AAT0.01560.04260.027034CAT0.00830.00830.00013AAC0.01930.01020.009135CAC0.0101 0.0037 0.00644AAG0.02570.00930.016436CAG0.00920.0037 0.00555ATA0.01560.01850.002937CTA0.0037 0.0111 0.00746ATT0.01650.07410.0576*38CTT0.0083 0.0241 0.01587ATC0.00730.01570.008439CTC0.0037 0.0083 0.00478ATG0.01560.01110.004540CTG0.0092 0.0046 0.00469ACA0.01470.00830.006441CCA0.0064 0.0019 0.004610ACT0.00640.01390.007542CCT0.0009 0.0083 0.007411ACC0.00550.00560.000143CCC0.0009 0.0037 0.002812ACG0.02850.00740.021144CCG0.0092 0.0019 0.007313AGA0.00550.01300.007545CGA0.0046 0.0028 0.001814AGT0.00280.00830.005646CGT0.00460.0056 0.001015AGC0.00370.0000.003747CGC0.0092 0.0009 0.008316AGG0.05970.00930.0504*48CGG0.07340.00560.0679*17TAA0.01010.04630.036249GAA0.03670.01200.024718TAT0.01280.05740.044650GAT0.01380.01300.000819TAC0.00830.01940.011251GAC0.01740.0019 0.015620TAG0.00730.01300.005652GAG0.0294 0.0046 0.024821TTA0.01380.09720.0835*53GTA0.0055 0.0102 0.004722TTT0.01560.17870.1631*54GTT0.0110 0.0241 0.013123TTC0.01290.01390.001055GTC0.0028 0.0056 0.002824TTG0.00920.00930.000156GTG0.0055 0.0009 0.004625TCA0.00180.00830.006557GCA0.0193 0.0037 0.015626TCT0.00370.02410.020458GCT0.0138 0.0019 0.011927TCC0.00460.00650.001959GCC0.0064 0.0009 0.005528TCG0.01650.00460.011960GCG0.0386 0.0009 0.037729TGA0.00460.00740.002861GGA0.0854 0.0083 0.0770*30TGT0.00640.00930.002862GGT0.0110 0.0083 0.002731TGC0.00180.00280.000963GGC0.0634 0.0046 0.0587*32TGG0.02750.00650.021164GGG0.02020.0009 0.0193* 表示该行所在的密码子在两类序列中的概率之差的绝对值大于0.05。（2） 特征密码子分析由特征密码子的定义知，表1 中的 *所在行的密码子就是特征密码子。 A类序列的特征密码子：GGA，CGG，GGC，AGG，其对应密码子序号为61，48，63，16。B类序列的特征密码子：TTA，TTT，ATT，TAT，其对应密码子序号为22，21，6，18。（二）模型的建立模型一：基于特征密码子分布概率构造分布辨析纸由特征密码子的分析可知，给定一个DNA序列，我们可以通过计算特征密码子出现概率的大小来进行判别分类。在A类中，提取八个密码子GGA，CGG，GGC，AGG，TTA，TTT，ATT，TAT，按A类中每一行DNA序列（1-10）T统计出它们的特征密码子频率，记为：（）其中。计算： 定义为最适A类密码子的概率； 定义为最适B类密码子的概率；建立平面直角坐标系，以为横坐标，为纵坐标，以（），其中,描绘出A类的样本点（见图一，用“*”表示）；同理，描绘出B类的样本点坐标，用“”号表示。建立平面直角坐标系，以为横坐标，为纵坐标，以（），其中,描绘出A类的样本点（见图一，用“*”表示）；同理，描绘出B类的样本点坐标，用“”号表示。再计算出未知DNA序列（21-40）的八个密码子（GGA，CGG，GGC，AGG，TTA，TTT，ATT，TAT）的统计概率，以平面坐标（），（如图一）。从图一中可清楚的看出，在由、构成的二维空间内，A类和B类呈现出非常好的聚类倾向。理想的情况是能够找到一条直线将A类和B类的点完全分开，那么自然的，与已知类同侧的点就属于该类。从图上可以看到这条直线位于过原点的角平分线附近，但是A类和B类都分别有一个点位于这条直线附近。所以为了准确起见，我们分别过原点和这两个点做了两条直线作为分界线，将二维空间分成了三个域，待判点位于右下角则属A类；位于左上域则属于B类；如果待判点位于两条分界线之间，则进入了本方案的辨析盲区，方案失效。将待判点描到图上，无须计算就可以将它们归类： 人工基因序列的分类结果：A类：34 35 23 27 29 37 22 B类：31 21 32 40 33 24 26 28 38 30不能判别的点：39 25 36这个分类方案简单，计算量小。如果将坐标和界线固定的画成一张坐标纸，就可以象正态概率纸一样方便的对DNA序列进行分类了。但它最主要的缺点就是存在判断盲区，而且，当我们把它应用在182个自然DNA序列的分类时（如图二），聚类就不那么明显了。所以该模型只能判断序列移位的人工模型，对自然DNA序列判断结果不理想。对自然DNA序列分类的结果：A类：点数143个B类：点数24个不能判别的点数：15个误判率在10%左右。模型二：图论最小生成树聚类法从图二中可以看出待判点聚类不明显，为了使以上结果进行精确地分类，为了构造出更好的分类模型，判别率更高，我们以模型一的分类结果（182个点的分类结果）为基础，研究了图论最小生成树的聚类法。基于模型一中对182个点的分类结果，采用贪婪算法先求出这182个点的最小生成树（minimum-weight spanning tree），再求出最小生成树的“颈”，按“颈”制定分类的方法和原则。对自然DNA序列分类的结果：A类：点数108个 B类：点数74个判别的主要步骤是：1、 在最小生成树上，找“主干”，最小生成树上边数最多的没有回路的链，称为“主干”。2、 求“子主干”，在最小生成树的主干上，“度”大于2的一点作为起点，肯定能够找到除主干本身外的没有回路的链，边树最多的就称为“子主干”，而它的边数就是该点的“深度”。3、 找“颈”，在的主干上找出“深度”大于2的点，如果这些点之间存在“深度”为零的点，则这些“深度”为零的点之间的边就称为“颈”。4、 根据最小生成树，去掉“颈”，就得到基因序列的分类。（如图三所示）通过以上步骤，找出了三个“颈”，分别是（72105）、（6499）、（28159），自然分成四类，再将相邻两类合并，得到A、B两类的分类结果，该方法也非常直观，易懂。该方法的缺点是误判概率较难估计，因为不能判别的点数以在“颈”附近相邻点数记，相对来讲，应该很少。模型三 多维向量空间的判别分析模型:图4 二维空间向量夹角示意图B（）0YXC（）A（）前面两个模型都有自身的局限性，判别分析就是要将与已知类别特征相近的物质归到一类，其关键就是要寻求待分样品和已知类别的相似特征。在这个模型中我们利用多维向量之间的夹角余弦差构造判别函数W(i)，对待判DNA序列片段进行归类。下面我们以二维向量空间为例介绍建模的原理：如图四所示在二维空间有两个类别已知的向量和向量，以及待判类别的向量，与、之间的夹角分别为和。若，则认为更加靠近，也就是在性质上与相近。若其中，向量、分别代表了A、B两类向量的重心位置，则可以断定属于A类；反之，则属于B类。当 时， 当 时， 根据平面几何的知识： 其中（i=1，2）分别为A、B、C各点的平面坐标值。令 则的判别式可写为： 当 时， ，属于A类；当 时，属于B类。若有i个待归类的向量则W可以表示为W（i），它是向量坐标的函数，称之判别函数。对于DNA序列片段的分类问题，虽然没有二维向量分类这么简单，但如果我们能用向量表示每个片段的结构，并且找出A、B两类片段的重心向量，就可以用上面同样的方法对DNA片段进行归类了。我们知道，DNA的排序是千变万化的，但是反映遗传信息的是由4种碱基所组成的64种密码子的有效排列，据此可以用一个64维向量表示一个DNA片段,向量中的每个元素对应一个密码子在这个片段中出现的频率，第i个片段的向量表示为： k=1,2,364表示第k种密码子在第i序列片段中出现的频率；那么所有属于A类和属于B类的DNA片段向量的集合就将整个64维空间分为两个互不重叠的域，我们姑且称之为A域和B域。若已知这两个域的重心向量，则可以分别求出待判别向量和A域重心向量间的夹角以及它和B域重心向量间的夹角，夹角越小表明它离所对应的重心向量越近。如果则表明这个向量在空间上更加接近A域，也就是这个向量所代表的DNA片段在特征上更类似于A类DNA片段，则可以认为这个片段是属于A类的。反之，则是属于B类的。值得注意的是，在这里“夹角”已经失去了它在二维和三维空间里的具体含义，而仅仅是表示两个向量间关系的“关系系数”，因此用夹角余弦所表达的意义和用夹角表达的意义是完全一致的，夹角越小，夹角余弦值越大，二者的关系也就越密切。所以我们用关系函数来判断待判向量与哪一类序列的关系更加密切。现在我们就以给出的A、B两类片段来估计这两个域的中心向量。对已给出的两类序列中的各个密码子的频率求平均,由这些平均频率所组成的向量就是这两个域的中心向量估计。两个域的重心向量和的估计为：第j(j=21,22,23,40)个待判序列的向量和A域重心向量的夹角余弦: 第j(j=21,22,23,40)个待判序列的向量和B域重心向量的夹角余弦: 第j个待判序列的判别函数为： (j=21,22,40)当 时, 向量与的夹角较小 ,第j个DNA片段属于A类;当 时, , 向量与的夹角较小 ,第j个DNA片段属于B类。四、结果分析模型一的求解及结果分析：问题（1）的分类结果： A类：22 23 27 29 34 35 37 B类：21 24 26 28 30 31 32 33 38 40不能判断的点： 25 36 39分辨率，由于没有发生错判的情况，所以误判率为零。问题（2）的分类结果：（如图二）A类：1 2 3 4 5 6 8 9 10 11 13 14 15 16 17 18 19 20 21 22 23 25 26 27 29 31 32 33 35 36 37 38 39 40 41 42 44 45 46 47 48 49 51 52 53 54 55 56 57 58 59 60 61 62 64 65 66 67 68 69 70 71 72 73 74 76 77 78 79 81 82 83 84 87 88 89 90 91 93 94 95 96 97 98 99 100 101 104 105 106 108 109 111 112 113 115 117 118 120 124 125 126 129 130 132 133 134 135 136 137 138 139 140 141 142 145 146 147 148 149 150 152 153 154 155 157 158 160 164 165 167 168 169 171 172 174 175 176 177 178 179 180 181B类：7 12 43 75 85 86 92 102 107 110 114 116 119 121 122 123 128 151 159 161 162 163 170 182不能判断的点：24 28 30 34 50 63 80 103 127 131 143 144 156 166 173这个方案简单、灵活，使用方便，计算量小。但由于存在判断盲区，而使它分辨率降低。所以这种方案适合于对基因片段的粗略的估计。模型二的求解及结果分析：(1) 自然基因中的分类结果为:A类：1 2 3 5 6 8 9 11 13 14 15 16 17 18 19 20 21 25 27 31 32 33 35 37 38 39 40 41 42 44 45 46 47 49 51 52 55 56 57 58 59 60 61 62 63 64 66 67 68 69 70 71 72 73 77 78 79 81 82 87 89 90 91 93 94 95 96 100 101 104 105 106 108 109 111 112 113 115 117 118 120 124 132 133 134 135 136 137 139 141 145 147 148 150 152 153 154 155 157 158 164 165 167 171 172 176 177 178 B类：4 7 10 12 22 23 24 26 28 29 30 34 36 43 48 50 53 54 65 74 75 76 80 83 84 85 86 88 92 97 98 99 102 103 107 110 114 116 119 121 122 123 125 126 127 128 129 130 131 138 140 142 143 144 146 149 151 156 159 160 161 162 163 166 168 169 170 173 174 175 179 180 181 182仔细分析模型一和模型二的结果，发现在模型二被判别为B类的点的集合中都包括了模型一中被判别为B类的所有点。说明模型二是在模型一的基础上更进一步的精确。这种方法简单易懂，非常直观。模型三的求解及结果分析：在这个模型中，我们将64维空间分成了两个分别