北邮生物信息基础 大作业报告

一、 疾病自动诊断问题

1. 问题分析

问题中已知100个腹泻患者的20项检查指标样本，并将其作为观测向量

x(x1,x2,...,x19,x20)T

其中xj为诊断结果值，要求据此诊断出新的腹泻病例的类型是细菌性还是病毒性。这属于模式识别中的非参数判别分类问题，一般而言需要根据观测样本提取特征，之后设计分类器并应用于新的数据的决策。非参数判别分类方法有线性分类器、非线性分类器、近邻法、人工神经网络、支持向量机等。在此我仅提出两个模型供参考。

模型一：利用线性分类器。假设判别函数是线性函数，用已知的100个训练样本去估计线性判别函数的参数。在选择和提取特征时，由于这20项指标对判别的贡献程度不同，所以应该根据重要性相应加权，这一步需要利用主成分分析的方法。随后我们需要求解判别函数，这一步我们采用较为成熟的Fisher线性判别函数分析方法，通过降低维度的方式，把20维空间映射到1维空间，找到一个最合适的投影轴，使两类样本（细菌性和病毒性）在该轴上投影的交迭部分最少，从而使分类效果为最佳。当维数和样本数都比较大时（比如题中维数=20，样本数=100），可以采用贝叶斯决策规则，从而获得一种在一维空间的“最优”分类器。

模型二：利用支持向量机（SVM）的方法。支持向量机是基于机器学习理论的一种新型的模式识别方法。在模式识别等领域获得了广泛的应用。其主要思想是：找到一个超平面，使得它能够尽可能多的将两类数据点正确地分开，同时使分开的两类数据点距离分类面最远（如图1.1.1）。我们采用该模型对问题进行建模。

图1.1.1 最佳超平面示意图

1

2. 模型建立

记n（这里n =100）个已知观测样本为(t1,g1),(t2,g2),...,(tn,gn)，其中tiR20，即为问题分析中的xi，为了不混淆改用ti。gi1为病毒性，gi1为细菌性。当样本线性或者非线性可分时，我们要找一个最优分类面wTxb0，其中

w,xR20，bR，w、b待定，满足如下条件

wTtib1,gi1 wtb1,g1ii即有gi(wtib)1，其中，满足方程wTtib1的样本成为支持向量。

要使两类总体到分类面的距离最大，则有

max212minw w2于是建立SVM的数学模型如下： 模型I

min12w 2s.t.gi(wtib)1，i1,2,...,n

求得最优值对应的w、b，可得分类函数

**g(x)sgn(w*Txb*)

模型I是一个二次规划模型。下面把模型I化为其对偶问题。 定义广义拉格朗日函数

n12T L(w,)wj1gwtibi2i1其中1,2,...,n，i0。

由Karush-Kuhn-Tucker互补条件，通过对w和b求偏导可得

nLwigiti0 wi1Lnigi0 bi12

得

wigiti

i1ngii1ni0

代入原始拉格朗日函数得

1nnLiijgigj(titj)

2i1j1i1其中(xixj)表示向量的内积。

于是模型I可以化为 模型II

n1nnmaxiijgigj(titj)

2i1i1j1ng0s.t.i1ii 0,i1,2,...,ni解此二次规划得到最优解，从而得权重向量wi*giti。

*n*ni1由KKT互补条件知

i*1giw*tib0

这意味着仅仅是支持向量（距离分类超平面为1）的输入点ti，使得i*为正，所有其它样本对应的i*均为0。选择*的一个正分量j*，并以此计算

bgjgii*titj

*i1n最终的分类函数表达式如下：

g(x)sgn(w*Txb*)

实际上，模型II中的(titj)是核函数的线性形式。核函数可以将原样本空间线性不可分的向量转化到高维特征空间中线性可分的向量。将模型II换成一般的核函数Kx,y，可得一般的模型：

3

模型III

1nnmaxiijgigjK(ti,tj)

2i1j1i1ng0s.t.i1ii 0,i1,2,...,ni分类函数表达式：

nn*g(x)sgnigiKti,xb*

i13. 模型测试

将新的测试样本t代入三个模型的分类函数g(x)中，按如下规则分类：

g(t)1，该样本点为病毒性； g(t)1，该样本点为细菌性。

~~~模型I其实已经可以较好地进行分类了，对已知样本点，支持向量机的分类方法错判率为0。模型II和模型III的计算需要较长时间，需要进行一定的优化。

二、 判断病毒是否变异

1. 问题分析

题中说CDC保留了15份诺如病毒的RNA序列，现又得到了一份来自患者P的诺如病毒的RNA序列。为判断该RNA序列是否变异，需要与之前的15份RNA序列进行多重序列比对。通过序列的多重比对，可以得到一个序列家族的序列特征。当给定一个新序列时，根据序列特征，可以判断这个序列是否属于该家族。

多重序列比对的最终目标是通过处理得到一个得分最高的序列对比排列，从而分析各序列之间的相似性和差异性。一般情况下是利用动态规划算法，但是对于序列条数比较多的时候往往需要大量时间。因此，可以考虑采用星形比对或树形比对等启发式方法。这类方法在绝大多数情况下计算结果接近于最优结果，但却可以大大减少计算时间，因此实际应用广泛。

由于比对整条RNA序列的复杂度较大，所需时间较长，因此我们可以比对病毒关键部分的RNA序列，或者说是与致病能力相关的RNA序列，从而确定是否变异。如果是小规模变异，由于比对的序列并非整条RNA序列，也比较容易找到变异的位点。

4

2. 模型建立

模型I：星形比对

星形比对的基本思想是：在给定的若干条序列中，选择一个核心序列，通过该序列与其他序列的两两比对，形成所有序列的多重比对α，从而使得α在核心序列和任何一个其他序列方向的投影是最优的两两比对。

星形比对的基本步骤如下： 1. 选择核心序列；

2. 计算与核心序列的两两比对；

3. 逐对聚合两两比对的结果，获得多重比对。

选择核心序列有两种方法：一种是尝试将每一个序列分别作为核心序列，按上述过程进行，去结果最好的一个；另一种方法是计算所有的两两比对，去下式值最大的一个

sim(s,s)

icic算法详细叙述如下：设s1,s2,...,sk是k条待比对的序列。假设一直核心序列是sc，c介于1到k之间，则可以利用标准的动态规划算法求出所有si与sc的最优两两比对，这个工作需要时间为O(kn2)（假设所有序列的长度为n）。将这些序列的两两比对聚集起来，并根据“只要是空位，则永远是空位”的原则。聚集过程从某一个两两比对开始，如s1与sc，然后逐步加上其他的两两比对。在这个过程中，逐步增加sc中的空位字符，以适应其他的比对，但决不删除sc中已存在的字符。假设在上述过程中的某一时刻，有一个由sc指导的、已经建立好的部分序列的多重比对，接下来就是加入一个新的、与sc两两比对的序列，如果需要，则插入新的空位字符。这个过程一直进行到所有的两两比对都加入后结束，所需的计算量为O(k2n)。因此，算法总的时间复杂度为O(kn2k2n)。 模型II：树形比对

在这里只给出树形比对的基本思路：把k个待比对的序列转化成具有k个叶节点的树，每个叶节点对应一个序列，如果给树的每个内部节点各指派一条序列，则可以计算树中每条边的权值，权值代表对应分支连接的两个序列之间的相似性，所有权值的和就是这棵树的得分。我们需要寻找一种树的内部节点序列赋予方式，使得树的得分最大。

5

3. 模型测试

对给定的15+1条序列（或关键部分序列）进行多重序列比对后，各序列已经对齐。对齐之后即可进行人工比对，可以很方便地看出碱基的变化情况。如果连续变化的碱基过多，则应判定为变异。至于究竟怎样的程度算是变异，则需要利用生物学的相关知识来确定。

在这里我们仅对星形比对法进行一个简单的举例验证。假设有5个序列分别为

s1= ATTGCCATT

s2= ATGGCCATT s3= ATCCAATTTT s4= ATCTTCTT s5= ACTGACC

经过计算，得到核心序列scs1，通过星形比对得到的结果为：

A T T G C C A T T - - A T G G C C A T T - - A T C – C A A T T T T A T C T T C – T T - - A C T G A C C - - - -

对齐后的序列很容易比较出哪里不一致，对于小规模变异，则更容易看出变异位点在哪里。

三、 基因编码区域识别

1. 问题分析

基因识别问题仍然是一个模式识别中的分类问题，理论上可以采用第一题中所述的几种方法。但是基因识别因其生物性，内在存在着一定的关联性，相邻密码子之间可能存在一定的依存关系，因此考虑使用马尔科夫模型完成编码区域和非编码区域的判别。而且题目中已经说明了RNA序列中基因编码区域的碱基排列顺序具有特定规律，那么就将每个碱基作为一个单位（而不是每个密码子），构建马尔科夫链。对于阶数问题，由于每个片段只有13个碱基，所以考虑采用一阶马尔科夫连进行建模（由于样本较少，高阶模型的转移矩阵会出现大量的0，且矩阵过大不便于计算）。

6

2. 模型建立

首先确定初始转移概率。根据基因编码区域和非编码区域的样本，统计出A、U、G、C在编码区和非编码区中出现的概率pC0和pN0（此处用频率估计概率）。之后计算相邻两个碱基转移的频率，以此作为一阶马氏链的概率转移矩阵。最后，若序列SA1A2...An利用公式

pC(S)pC0(A1)p(Ai|Ai1,SC)

i2nnpN(S)pN0(A1)p(Ai|Ai1,SN)

i2计算序列s1、s2在编码区和非编码区出现的概率pC(S)和pN(S)。若pC(S)>

pN(S)，则判定该序列为编码区域，否则为非编码区域。

3. 模型测试

基因编码区域片段为： AUGGGCAAAUAGC AUGGAUAGCGCAA GCGCAAGAAUGUA AGCGCAAGAAUGU UAGCGCAAGAAUG 待判断的序列片段为：

基因非编码区域片段为： CUCUCUCACACGU CACACGUCUCUCU CACACGUCUCUCU CACACUCUGUCUC GUCUCACACCUCU

S1：GAAUGUAGCGCAA S2：CUCUCUCACACGU

编码区A、U、G、C出现的频数为 碱基 频数

A 25

U 10

G 20

C 10

总数为65，得编码区初始转移概率pC0为 碱基 概率 碱基 频数

A 0.3846 A 10

U 0.1538 U 20

7

G 0.3077 G 5

C 0.1538 C 30

非编码区A、U、G、C出现的频数为

得非编码区初始转移概率pN0为 碱基 概率

A 0.1538

AU 0.3077

UGG 0.0769

CC 0.4615

编码区一阶马氏链的概率转移矩阵为

A0.39130.30430.30430 U0.444400.55560PCG0.21050.10530.15790.5263C0.555600.44440非编码区一阶马氏链的概率转移矩阵为

AUGCA0001.0000 U000.06250.9375PNG01.000000C0.34880.51720.10340.0345本题中，n=13。将S1、S2分别代入两个公式，得

pC(S1)pC0(A1)p(Ai|Ai1,S1C)

i2n=0.3077*0.2105*0.3913*0.3043*0.5556*0.1053*0.4444*0.3043*0.5263 *0.4444*0.5263*0.5556*0.3913 =1.6329e-06 =1.6329106

pN(S1)pN0(A1)p(Ai|Ai1,S1N)=0（因为连乘中有0，故结果为0）

i2nnpC(S2)pC0(A1)p(Ai|Ai1,S2C)=0（因为连乘中有0，故结果为0）

i2npN(S2)pN0(A1)p(Ai|Ai1,S2N)

i2=0.1538*0.5172*0.9375*0.5172*0.9375*0.5172*0.9375*0.3488*1.0000 *0.3488*1.0000*0.1034*1.0000 =2.2056e-04 =2.2056104

8

因为pC(S1)pN(S1)，pC(S2)pN(S2)，所以我们判定S1序列属于编码区，

S2序列属于非编码区。

四、 课堂内容回顾和建议

本学期我选修了《生物信息基础》这门课。这门课以分子生物学为背景，以信息和计算机技术为手段，研究生物学信息的组织、传递和表达规律等问题，并从中发现生物遗传、变异、进化等规律。

首先，老师先对生物方面的基本知识进行了讲解，主要是遗传方面的知识。比如DNA转录成为mRNA，mRNA翻译成为蛋白质，以及有关DNA、RNA、蛋白质序列的一些问题。这些基础的问题我们在高中阶段已经有所涉及，相当于借此复习一遍，为后续的课程打下基础。

之后的课上，老师介绍了一些常用的生物数据库。这些数据库包括基因组数据库、核酸序列数据库、蛋白质序列数据库、蛋白质结构数据库等。我们可以从这些数据库中找到各种生物的原始数据。接下来便自然想到要对这些数据进行分析。由于生物遗传信息一般存在于DNA、RNA、蛋白质序列中，所以对序列进行分析便是这门课的重点之一。

在序列分析中，经常需要判断序列的相似性，所以，对序列进行比对便是最基本的操作。一般来说，序列是由碱基构成的。序列比对时，可以通过字符匹配、替换、插入、删除字符等操作，使两序列长度相等，便于对比出碱基的异同之处，从而分析生物学信息。为了判断比对的效果，还引入了打分函数和打分矩阵机制。除了两两比对，我们还可以进行多序列比对，以便于判断某个序列是否属于某一族，这在实际的生物学中有很高的应用价值。

序列比对的方法也有很多种，这里就涉及到算法的问题。老师向我们介绍了很多序列比对的算法，如BLAST搜索算法、FASTA算法等。这使我了解到不同学科之间可以存在这种交叉现象，计算机学科的算法可以应用到生物学中帮助分析遗传信息序列。

序列分析之后的另一个重点是基因组与基因识别问题。首先，老师介绍了原核生物和真核生物的基因组相关知识。之后提出了一个实际问题：我们经常需要分析某一个基因是属于功能区还是非功能区、编码区还是非编码区。这其实是一种模式识别问题。也就是说，我们需要根据已知样本做出一个分类器，来对新的样本进行分类判别。老师在课堂上提出了基于贝叶斯判别的朴素贝叶斯方法和马尔科夫链两种模型，并介绍了它们的优缺点。这使我了解到，在面对具体的生物学问题时，应选择合适的方法进行分析，这样才能达到比较好的效果。

除了刚才的两种模型外，在基因识别问题中还有一个常用的模型——隐马尔科夫模型，它用来描述一个含有隐含未知参数的马尔可夫过程。这种模型可以方

9

便地对DNA的编码区和非编码区建模，分析进化和发育问题，或者对蛋白质结构进行预测。

最后，老师简要介绍了系统生物学的一些内容和相关算法，包括如何构建进化树，以及对进化树进行置信检验的方法。还介绍了基因表达数据的相关概念，以及一些重要的算法——k均值聚类算法、最近邻法、主成分分析法等，这些都算模式识别中重要的算法。

我本人对各种模式识别的算法比较感兴趣，而这门课讲的许多算法，如动态规划、马尔科夫链、隐马尔科夫模型、聚类算法、主成分分析法等，都属于模式识别的范畴，只不过是应用在了生物信息的分析中。这些算法在建模中有重要的意义，是数学方法应用在具体问题中的典范。这门课使得我了解了许多算法的基本原理和操作步骤，以及它们的应用范畴。

由于各种算法本质上来说仍属于数学的范畴，所以本身具有一定的难度，有时会觉得没有头绪、难以理解，然后就容易听不下去了，这对于我们的学习是不利的。所以我希望老师在解释这些算法的时候，可以加入一些演示和提问的环节，充分与学生互动，让学生亲自到黑板上操作，参与到过程中来，这样更利于我们掌握课堂知识。

总的来说，我在这门课上收获颇多，不仅学习到了生物学的相关知识，还了解了许多模式识别方面的算法以及它们的实际应用。这些对我对算法的理解和实际建模都很有帮助。

五、 参考文献

[1]边肇祺，张学工等．模式识别（第二版）．北京：清华大学出版社．2000． [2]孙啸，陆祖宏，谢建明．生物信息学基础．北京：清华大学出版社．2005． [3]蔡禄．生物信息学教程．北京：化学工业出版社．2007．

[4]司守奎，孙玺菁．数学建模算法与应用．北京：国防工业出版社．2011．

10