DNA序列分类模型

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个。具体情况见文中答案。 我们还对三种分类方法进行了类比，认为模型二、三方法新颖独特，结果稳定，它们是一种较好的分类方法。并且对各种计算结果进行误差分析和检验等工作。 1 一、 问题的重述 本问题为一个DNA序列分类问题。假定已知两组人工已分类的DNA序列(20个已知类别的人工制造的序列)，其中序列标号1—10 为A类，11-20为B类。要求我们从已经分类了的DNA序列片段中提取共同特征构造分类方法，并评价所用分类方法的好坏，从而构造或选择一种较好的分类方法。测试对象是20个未标明类别的人工序列(标号21—40)和182个自然DNA序列。 二、模型的假设及符号说明 1、名词解释: 碱基: 在生物学中，用A，T，C，G四个字符代表组成DNA序列的四种碱基; 密码子:在遗传学中每三个碱基的组合被称为一个密码子，可以编码一个氨基酸，共有64个，还可以由密码子组成20个氨基酸。 特征密码子:若类中某个密码子出现的概率p与类中同一个密码子出现的概率ABA p之差的绝对值大于0.05，并且，p> p，定义该密码子为类特征密码子。同理，ABA B 可定义B类特征密码子。 2、模型假设: ? 所给的DNA序列片段中没有断句和标点符号。对密码子的计数没有固定的起点 和终点。 ? 假定各密码子近似呈现正态统计规律。 ? DNA序列的某些片段之间具有一定的相关性。 ? DNA序列不构成环状，划分密码子时，碱基剩余个数小于3时舍去不考虑。 3、符号说明: p(k)------------------- 表示第k种密码子在第i个序列片段中出现的频率;i H---------------------- 类特征密码子; Ai R---------------------- 类特征密码子; Bi P(H)--------------------最适A类密码子的概率; i P(R)---------------------最适B类密码子的概率; i ,-------------------------误判率; W(i)---------------------空间向量夹角余弦之差; 2 三、模型的分析和建立 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(p) B(p) |p-p| A(p) B(p) |p-p| ABABABAB1 AAA 0.045 0.045 0.0004 33 CAA 0.0138 0.0065 0.0073 2 AAT 0.0156 0.0426 0.0270 34 CAT 0.0083 0.0083 0.0001 3 AAC 0.0193 0.0102 0.0091 35 CAC 0.0101 0.0037 0.0064 4 AAG 0.0257 0.0093 0.0164 36 CAG 0.0092 0.0037 0.0055 5 ATA 0.0156 0.0185 0.0029 37 CTA 0.0037 0.0111 0.0074 6 ATT 0.0165 0.0741 0.0576* 38 CTT 0.0083 0.0241 0.0158 7 ATC 0.0073 0.0157 0.0084 39 CTC 0.0037 0.0083 0.0047 8 ATG 0.0156 0.0111 0.0045 40 CTG 0.0092 0.0046 0.0046 9 ACA 0.0147 0.0083 0.0064 41 CCA 0.0064 0.0019 0.0046 10 ACT 0.0064 0.0139 0.0075 42 CCT 0.0009 0.0083 0.0074 11 ACC 0.0055 0.0056 0.0001 43 CCC 0.0009 0.0037 0.0028 12 ACG 0.0285 0.0074 0.0211 44 CCG 0.0092 0.0019 0.0073 13 AGA 0.0055 0.0130 0.0075 45 CGA 0.0046 0.0028 0.0018 14 AGT 0.0028 0.0083 0.0056 46 CGT 0.0046 0.0056 0.0010 15 AGC 0.0037 0.000 0.0037 47 CGC 0.0092 0.0009 0.0083 16 AGG 0.0597 0.0093 0.0504* 48 CGG 0.0734 0.0056 0.0679* 17 TAA 0.0101 0.0463 0.0362 49 GAA 0.0367 0.0120 0.0247 3 18 TAT 0.0128 0.0574 0.0446 50 GAT 0.0138 0.0130 0.0008 19 TAC 0.0083 0.0194 0.0112 51 GAC 0.0174 0.0019 0.0156 20 TAG 0.0073 0.0130 0.0056 52 GAG 0.0294 0.0046 0.0248 21 TTA 0.0138 0.0972 0.0835* 53 GTA 0.0055 0.0102 0.0047 22 TTT 0.0156 0.1787 0.1631* 54 GTT 0.0110 0.0241 0.0131 23 TTC 0.0129 0.0139 0.0010 55 GTC 0.0028 0.0056 0.0028 24 TTG 0.0092 0.0093 0.0001 56 GTG 0.0055 0.0009 0.0046 25 TCA 0.0018 0.0083 0.0065 57 GCA 0.0193 0.0037 0.0156 26 TCT 0.0037 0.0241 0.0204 58 GCT 0.0138 0.0019 0.0119 27 TCC 0.0046 0.0065 0.0019 59 GCC 0.0064 0.0009 0.0055 28 TCG 0.0165 0.0046 0.0119 60 GCG 0.0386 0.0009 0.0377 29 TGA 0.0046 0.0074 0.0028 61 GGA 0.0854 0.0083 0.0770* 30 TGT 0.0064 0.0093 0.0028 62 GGT 0.0110 0.0083 0.0027 31 TGC 0.0018 0.0028 0.0009 63 GGC 0.0634 0.0046 0.0587* 32 TGG 0.0275 0.0065 0.0211 64 GGG 0.0202 0.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， TATT，TAT，按A类中每一行DNA序列(1-10)统计出它们的特征密码子频率，记为:P,P,？P(i,1,2,？10)()其中。 i1i1ij 4 P(H),P计算: 定义为最适A类密码子的概率; ,iij1j, 8 P(R),P 定义为最适B类密码子的概率; ,iij5j, P(H)P(R)P(H),P(R)建立平面直角坐标系，以为横坐标，为纵坐标，以()，其ii(i,1,2,？10)中,描绘出A类的样本点(见图一，用“*”表示);同理，描绘出B类的样本 4 点坐标，用“?”号表示。建立平面直角坐标系，以为横坐标，为纵坐标，以P(H)P(R)()，其中,描绘出A类的样本点(见图一，用“*”表示);同(i,1,2,？10)P(H),P(R)ii 理，描绘出B类的样本点坐标，用“?”号表示。再计算出未知DNA序列(21-40)的八个密码子(GGA，CGG，GGC，AGG，TTA，TTT，ATT，TAT)的统计概率，以平面坐标()，P(H),P(R)ii(如图一)。 P(H)P(R)从图一中可清楚的看出，在由、构成的二维空间内，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序列进行分类了。但它最主要的缺点就是存在判断盲区，而且， 5 当我们把它应用在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的一点作为起点，肯定能够 找到除主干本身外的没有回路的链，边树最多的就称为“子主干”，而它的边数就 是该点的“深度”。 6 3、 找“颈”，在的主干上找出“深度”大于2的点，如果这些点之间存在“深度”为 零的点，则这些“深度”为零的点之间的边就称为“颈”。 4、 根据最小生成树，去掉“颈”，就得到基因序列的分类。(如图三所示) 通过以上步骤，找出了三个“颈”，分别是(72,,105)、(64,,99)、(28,,159)，自然分成四类，再将相邻两类合并，得到A、B两类的分类结果，该方法也非常直观，易懂。该方法的缺点是误判概率较难估计，因为不能判别的点数以在“颈”附近相邻点数记，相对来讲，应该很少。 模型三 多维向量空间的判别分析模型: 前面两个模型都有自身的局限性，判别分析就是要将与已知类别特征相近的物质归到一类，其关键就是要寻求待分样品和已知类别的相似特征。在这个模型中我们利用多维向量之间的夹角余弦差构造判别函数W(i)，对待判DNA序列片段进行归类。下面我们以二维向量空间为例介绍建模的原理: 如图四所示在二维空间有两个类别已知的Y 向量和向量，以及待判OA(a,a)OB(b,b)1212 (b,b)B12 (OCOAOB类别的向量，与、之OC(c,c)12) C1,C2C() ,,,,,间的夹角分别为和。若，则认为 ,A(a1,a2) 7 , X 0 图4 二维空间向量夹角示意图 OCOAOCOAOA更加靠近，也就是在性质上与相近。若其中，向量、OB分别代表 OC了A、B两类向量的重心位置，则可以断定属于A类;反之，则属于B类。 当 时， ,,,cos,,cos,,0 当 时， ,,,cos,,cos,,0 根据平面几何的知识: ac，acbc，bc11221122cos,,cos,, 22222222a，a,c，cb，b,c，c12121212 其中(i=1，2)分别为A、B、C各点的平面坐标值。a,b,ciii 令W,cos,,cos, OC则的判别式可写为: ac，acbc，bc11221122W,, 22222222a，a,c，cb，b,c，c12121212 OCW,0当 时， ,,,，属于A类; OCW,0,,,当 时，，属于B类。 若有i个待归类的向量则W可以表示为W(i)，它是向量坐标的函数，称之判别函数。 对于DNA序列片段的分类问题，虽然没有二维向量分类这么简单，但如果我们能用向量表示每个片段的结构，并且找出A、B两类片段的重心向量，就可以用上面同样的方法对DNA片段进行归类了。 我们知道，DNA的排序是千变万化的，但是反映遗传信息的是由4种碱基所组成的64种密码子的有效排列，据此可以用一个64维向量表示一个DNA片段,向量中的每个元素对应一个密码子在这个片段中出现的频率，第i个片段的向量表示为: P,[p(k)] k=1,2,3„64 ii p(k)表示第k种密码子在第i序列片段中出现的频率;i 那么所有属于A类和属于B类的DNA片段向量的集合就将整个64维空间分为两个互不 8 重叠的域，我们姑且称之为A域和B域。若已知这两个域的重心向量，则可以分别求出待判别向量和A域重心向量间的夹角以及它和B域重心向量间的夹角，夹角越小表明它离,, 所对应的重心向量越近。如果则表明这个向量在空间上更加接近A域，也就是这个,,, 向量所代表的DNA片段在特征上更类似于A类DNA片段，则可以认为这个片段是属于A类的。反之，则是属于B类的。值得注意的是，在这里“夹角”已经失去了它在二维和三维空间里的具体含义，而仅仅是表示两个向量间关系的“关系系数”，因此用夹角余弦所表达的意义和用夹角表达的意义是完全一致的，夹角越小，夹角余弦值越大，二者的关系也就越密切。所以我们用关系函数来判断待判向量与哪一类序列的关系更加密W(j),cos,,cos, 切。 现在我们就以给出的A、B两类片段来估计这两个域的中心向量。对已给出的两类序列中的各个密码子的频率求平均,由这些平均频率所组成的向量就是这两个域的中心向量估计。两个域的重心向量和的估计为: PPAB 101010，，1ˆ ,Pp(1),p(2),...,p(64),,,Ammm,,10m,1m,1m,1，， 202020，，1ˆ ,Pp(1),p(2),...,p(64),,,Bnnn,,10n,11n,11n,11，， 第j(j=21,22,23,„40)个待判序列的向量和A域重心向量的夹角余弦: 64ˆp(k)P(k),jAk,1cos, ,164642，，22ˆp(k)P(k),,jA,,kk,,11，， 第j(j=21,22,23,„40)个待判序列的向量和B域重心向量的夹角余弦: 64ˆp(k)P(k),jBk,1cos, , 164642，，22ˆp(k)P(k),,jB,,kk,,11，， 第j个待判序列的判别函数为: 6464ˆˆp(k)P(k)p(k)P(k),,jAjBk,k,11W(j),, (j=21,22,„40)116464646422，，，，2222ˆˆp(k)P(k)p(k)P(k),,,,jAjB,,,,kk,,kk,,1111，，，， 9 ˆ当 时,, 向量与的夹角较小 ,第j个DNA片段属P,[p(k)]W(j),0,,,PjA 于A类; ˆ当 时, , 向量与的夹角较小 ,第j个DNA片段属P,[p(k)]W(j),0,,,PjB 于B类。 四、结果分析 模型一的求解及结果分析: 问题(1)的分类结果: A类:22 23 27 29 34 35 37 B类:21 24 26 28 30 31 32 33 38 40 不能判断的点: 25 36 39 17分辨率，由于没有发生错判的情况，所以误判率为零。,,，100%,85%20 问题(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 181 B类: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 168 ,,，100%,92.3%182 这个方案简单、灵活，使用方便，计算量小。但由于存在判断盲区，而使它分辨率降低。所 以这种方案适合于对基因片段的粗略的估计。 模型二的求解及结果分析: (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 10 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维空间分成了两个分别以为重心的、互斥的域，P和PBB 从而对基因片段进行分类。利用Matlab软件，我们编制计算机程序20个人工合成基因片段以及182个自然DNA序列片段分别进行了判类计算(有关程序请参见附录程序清单～)。同时为了评估这种分类方法的优劣，我们还对已知类别的20个基因片段进行了重新判类分析，结果如下: 判别函数的值: W(i)=[ 0.6526 0.6317 0.7481 -0.0923 0.7198 0.6145 0.5194 0.5331 0.6313 0.6246 -0.6162 -0.5975 -0.6681 -0.6274 -0.7024 -0.5414 -0.0072 -0.6082 -0.6681 -0.6309] (i)>0时，第i个序列片段属于A类; 根据判类准则:当W 当W(i)<0时，第i个序列片段属于B类 已知人工序列片段的重新分类结果如下: A类:1， 2， 3， 5， 6， 7， 8， 9，10; B类:4，11，12，13，14，15，16，17，18，19，20; 假设原来给定的分类是准确的话，显然新方案在第四段处发生了误判。由判别函数值我们也看到W(4)和W(17)的绝对值相对其他值而言要小的多，因此，虽然第十七片段并没有发生误判，但我们仍然认为是不可靠的。同时我们也必须看到，由于样本非常有限，每类给出的样本分别只有10个，从这么少的样本中抽取的共同特征是有限的。因此，这种分类 2方法的误判率。和模型二相比较，这个误判率是可以接受的。而且,,，100%,10%20 我们相信，随着样本数的增加，重心向量的估计值会更靠近真实值，从而使误判率进 一步降低，因此，多维空间的判别分析不失是一种对基因片段分类好的方法。 采用与上述相同的计算程序我们对另外二十个不知类别的人工基因片段进行了分类计算。为了更直观的反映判类函数在分类计算中的作用，我们在坐标纸上描出了每一片段的判别函数值，如图五所示，凡是处于零线以上的均属于A类，以下的则属于B类，零线两侧非常近的点代表的则是那些可能发生误判的片段。由图可知第21到第40人工片段 11 图5 人工DNA判别函数的点阵图 分类结果: A类: 22, 23，25，27，29，34，35，36，37, 39; B类: 21 ,24 ,26 ,28 ,30 ,31 ,32 ,33 ,38 ,40。 对182种自然DNA序列片段的分类结果如下: A类共120个片段: 1 2 3 4 5 6 8 9 11 13 14 15 16 17 18 19 20 21 25 27 29 31 32 33 35 36 37 38 39 41 42 44 45 46 47 49 51 52 53 55 56 57 58 59 60 61 62 63 64 66 67 68 69 70 71 72 73 77 78 79 81 82 83 87 88 89 90 91 93 94 95 96 97 99 100 101 104 105 106 108 109 111 112 113 115 117 118 120 124 125 126 132 133 134 135 136 137 139 141 145 147 148 149 150 152 153 154 155 157 158 160 164 165 167 169 171 172 176 177 178 B类共62个片段: 7 10 12 22 23 24 26 28 30 34 40 43 12 48 50 54 65 74 75 76 80 84 85 86 92 98 102 103 107 110 114 116 119 121 122 123 127 128 129 130 131 138 140 142 143 144 146 151 156 159 161 162 163 166 168 170 173 174 175 179 180 181 182 鉴于篇幅有限，其判别函数值列于附录上，这里仅将数据点阵绘出，如图六所示。同 样位于零线以下的点代表属于B类的片段，而位于零线以上的点则代表了属于A类的点。由 图上看出这些点中决大多离零线较远，只有少数几个点位于零线附近，是可能的误判点。我 们取0.02作为判别可靠性的阀值，即当|W(i)|>0.02 (i=1，2，…182)时，我们就认为 对它的判类是可靠的，否则就是不可靠的。在图上做出这两条阀值线，可以看到落在这两条 线之间的点很少。从实际计算得到的判别函数值上看，仅有五条基因片段所对应的判类函数 值落在这两个阀值之间，假使它们全都被误判了类别，那么总误判率 5,,10%，这比我们开始预计的要小的多。这又一次证明了,,，100%,2.75%182 采用多维向量来辩析DNA序列片段是一种非常有效的方法 图6 自然DNA判别函数的点阵图 13 六、模型的优缺点 模型的优点: (1)模型一 如前所述，用模型一对DNA序列片段进行分类方便，快洁，计算量小，虽存在辨析盲区，但不会发生误判的情况，故用第一中方案进行初步的估计分类是能够满足要求的。 (2)模型二 采用图论的方法，对图中的点画出最小支撑树，分类直观，判别率高达90%以上，误判率控制在0.05-0.1的范围内。 (3)模型三采用多维向量空间的夹角余弦值来判别序列的类别。从序列的全局密码子概率分析数据，保留了数据的完整性，而且判别准确率高，更值得一提的是它的判断结果和模型二的结果惊人的相似，也从另一个角度证明了这两个方案的准确性。 参考文献: 《基因》 [美] B.卢因 《大学生数学建模竞赛辅导教材》 《数学的实践与认识》 《概率统计计算方法》 肖如云编 南开大学出版社 14