《序列分析与联配》由会员分享，可在线阅读，更多相关《序列分析与联配（79页珍藏版）》请在金锄头文库上搜索。
1、第三章 序列分析与联配第一节第一节 序列组成和单一序列分析序列组成和单一序列分析第二节第二节 序列联配序列联配第三节第三节 数据库搜索引擎数据库搜索引擎BLAST和和FASTA应用应用第四节第四节 寡核苷酸设计寡核苷酸设计序列分析是生物信息学最主要的研究内容之一，它可以分为两个主要部分：一是序列组成（特别是涉及到基因组层次上）分析，二是序列之间的比较分析。两条序列或多条序列间的比对或联配(alignment)的目的，是对它们的序列相似性进行评估，找出这些序列中结构或功能相似性区域等。通过联配未知序列与已知序列(其功能或结构等已知)的相似程度，我们可以判断或推测未知序列的结构与功能。 第一节第一
2、节 序列组成及单一序列分析序列组成及单一序列分析一、碱基组成一、碱基组成二、碱基相邻频率二、碱基相邻频率三、同向重复序列分析三、同向重复序列分析四、四、DNA序列的几何学分析序列的几何学分析Z曲线曲线一、碱基组成一、碱基组成DNA序列一个显而易见的特征是四种碱基类型的分布。尽管四种碱基的频率相等时对数学模型的建立可能是方便的，但几乎所有的研究都证明碱基是以不同频率分布的。表3.1包含了9条完整DNA分子序列的资料，表3.2的数据来自两个胎儿球蛋白基因(Gr和Ar)，每个基因具有三个外显子和两个内含子(shen等1981)。这两个例子说明序列内和序列间碱基具有不同的频率。在基因每一侧的500 个
3、任意碱基区域被称为“侧翼”，基因间区域是指两个基因间的其余序列。 二、碱基相邻频率二、碱基相邻频率分析DNA序列的主要困难之一是碱基相邻的频率不是独立的。碱基相邻的频率一般不等于单个碱基频率的乘积：如果Pu是序列中碱基u的频率，且Puv为两个相邻碱基u和v的频率，则 Puv PuPvNussinov(1984)研究了两碱基相邻的频率(表3.3)。数据来自166个脊椎动物的DNA序列，总长136731个碱基。表中的比值为16种二个碱基相邻的频率除以相应的单个碱基频率的乘积。 作为一个特别的例子，图3.1给出了鸡血红蛋白链的mRNA编码区的438个碱基。表3.4列出了4种碱基和16种两碱基的数目。
4、将该表看作44的表，计算行列独立性的卡方统计量，得到x2=59.3（x20.05,9=16.92）表明行(第一碱基)列(第二碱基)之间存在明显的关联。 在编码区，存在某种约束来限制DNA序列编码氨基酸。在密码子水平上，这一约束与碱基相邻频率有关。表3.5列出了遗传密码和图3.1序列中各密码子数量。尽管数目很小，难以作出有力的统计结论，但编码同一氨基酸的不同密码子(同义密码子)好像不是等同存在的。这种密码子偏倚必定与两碱基相邻频率水平有关。表3.5还清楚地表明，由于密码子第3位置上碱基的改变常常不会改变氨基酸的类型，因而对第3位置上碱基的约束要比第 2位碱基小得多。 相邻碱基之间的关联将导致更远
5、碱基之间的关联，这些关联延伸距离的估计可以从马尔科夫链(Markov chain)理论得到(Javare和Giddings，1989)。在不援引任何生物学机制的情况下，第k阶马尔科夫链假定在序列中某一位置上碱基的存在只取决于前面k个位置上的碱基。一阶链假定一个特定碱基存在于位置i的概率只取取决于在位置i-1的4种碱基概率。相互独立的碱基所组成的序列将与0阶马尔科夫链相对应。阶可以通过似然法估计。同时，马尔科夫链分析更适应于基因组水平，而非单一序列(基因)。三、同向重复序列分析三、同向重复序列分析 除了分析整个序列碱基关联程度的特征外，我们常对寻找同向重复序列(direct repeats)之类
6、的问题感兴趣。Karlin等(1983)给出了完成这一分析的有效算法。该法采用由特定的几组碱基字母组成的不同亚序列或称为字码(word)。只需要对整个序列搜索一次。给一碱基赋以值,例如A、C、G、T的值为0、1、2、3。由X1、X2、.、Xk共k个字母组成的每一种不同的字码按 计算字码值。这些值的取值范围为1到4k。例如，5字码TGACC的值为1+344+243+ 042+141+140=459。可先从低k值的字码开始搜索。记录序列中每一个位置k字码的字码值。只有在发现k字码长度重复的那些位置考虑进行长度大于k的字码搜索。 表3.6列出了序列TGGAAATAAAACGTAAGTAG中所有碱基2
7、字码(k=2)的初始位置和字码值。对于完全重复、长度大于2的同向重复或亚序列的搜索可只限于2字码重复的初始位置。在本例中只有4个重复的2碱基重复序列。例如，在位置4、5、8、9、10和15均发现了字码值为1的碱基重复序列。从有重复的第2个碱基为起点的3字码值及位置列于表3.7，其中发现字码值为1、45和49的序列有重复。以每一重复的3碱基为起点的4字码搜索未能发现更长的重复序列。因此最长的同向重复为4、8、9位置上的AAA，13、17位置上的GTA以及7、14位置上的TAA。同样对图3.1鸡球蛋白DNA序列进行同向重复序列搜索，一些最长同向重复序列列于表3.8。 Karlin等(1983)提出
8、了序列内存在的最长同向重复序列的统计显著性评价方法。在核苷酸的位置为独立的假定下(相当于阶次为0的马尔科夫链)，长度为n的序列中，最长同向重复Ln的期望长度和方差为：其中，P为序列中碱基频率的平方和：用尽可能接近最大长度的期望均值的字码(即RL) 来开始同向重复序列的搜索计算可能节省计算量。 方程3.1可以用一个近似方法来验证以上统计假说。假定同向重复序列的长度呈正态分布。对于图3.1鸡蛋白序列，A、C、G、T四个碱基的次数分别为87、144、118和89，因而P=0.2614，最长重复序列的期望长度为8.13且具有期望方差0.9138。根据95%的正态分布概率，理论上可以预期最长同向重复序列
9、不超过10。 四、四、DNA序列的几何学分析序列的几何学分析Z曲线曲线DNA序列实际上是一种用4种字母表达的“语言”，只是其“词法”和“语法”规则目前还没有搞清楚。人类的语言有文字、声音两种基本表现形式，此外还有手语、旗语甚至图画语等特殊表达形式。同样，DNA序列作为一种语言，其表达形式也不是唯一的。传统上，DNA序列是用4种字母符号表达的一维序列。这是一种抽象形式，适合于存储、印刷和代数算法的处理，包括比较、排列和查找特殊序列等。我国学者张春霆等开展了DNA序列三维空间曲线表示形式，即DNA序列几何表示形式的研究。几何形式虽然与符号形式完全等价，但显示了DNA序列的新特征。两种形式各有其特点
10、，相互补充。这一新方法，为解读DNA序列信息提供了崭新的手段。 他们的研究始于对4种碱基对称性的观察，提出了用正面体表示碱基对称性。1994年，他们利用这种形式来表示任意长度的DNA序列。现将这种序列表示方法简述如下。 考察一个长为L的单股DNA序列，方向(53或35)不限。从第一个碱基开始，依次考察此序列，每次只考察一个碱基。当考察到第n个碱基时(n=1,2，.,L)，数一下从1到n这个子序列中四种碱基各自出现的次数。设4种碱基A、C、G、T出现的次数分别以An、Cn、Gn、Tn表示之，这里下标“n”是表明这些整数是从1到n这个子序列中数出来的，如图3.2所示。显然，它们都是正整数。根据正四
11、面体的对称性可以证明，在正面体内存在唯一的一个点Pn与这四个正整数对应。点Pn构成了四个正整数的一一对应映射。点Pn坐标可用四正整数表达： 方程3.2其中xn，yn和zn为点Pn的三个坐标分量。当n从1跑到L时，我们依次得到P1，P2，P3，.，PL共L个点。将相邻两点用适当的曲线连接所得到的整条曲线，就称为表示DNA序列的Z曲线。可以证明，Z曲线与所表示的DNA序列是一一对应的，即给定一DNA序列，存在唯一的一条Z曲线与之对应；反之，给定一条Z曲线，可找到唯一的一个DNA序列与之对应。换言之，Z曲线包含了DNA序列的全部信息。Z曲线是与符号DNA序列等价的另一种表示形式，一种几何形式。可以通
12、过Z曲线对DNA序列进行研究。 Z曲线的三个分量(方程3.2)具有明确的生物学意义：xn表示嘌呤/嘧啶碱基沿序列的分布。当从1到n的这个子序列中(图3.2)嘌呤碱基多于嘧啶碱基时，xn0，否则，xn0，否则，yn0，否则zn0，当两者相等时，zn=0。这三种分布是相互独立的，表现在以下事实上：任何一种分布不能由其它两种分布的线性叠加表示出来。给定的DNA序列唯一地决定了这三种分布；三种分布唯一地描述了DNA序列。对DNA序列的研究就是通过对这三种分布的研究来进行。从方法学的角度来看，这是DNA序列的一种几何学研究途径。 图3.3给出了大肠杆菌ayoP基因族序列Z曲线的三个分量，即三种分布图。该
13、基因族包含了大肠杆菌5 个基因aroP，A，aceFE，aceF和lpd，总长度为9501bp，分别编码芳香族氨基酸运输蛋白aroP，蛋白质A(功能不详)和三种酶，即丙酮酸脱氢酶，二氢硫辛酰基转移酶和二氢硫辛酰脱氢酶。它们位于此序列的0039-1406，1947-2654，2870-5527，5545-7434，7759-9183区间。在图中X轴的下方的基因排列图上已分别用阴影标出相应基因。在这些基因之间有三个启动子区（pm1、pm2和pm3），其中aceE和aceF基因属于ace操纵子，共用一个启动子。三个启动子区亦在图中标出。非常令人感兴趣的是，在5个编码区，Z曲线的z分量基本上都是单调下
14、降的，而在三个启动子区基本上都是单调上升的。x，y分量亦有变化，但不如z分量明显。在上升、下降的交界处，Z曲线均发生了重大的转折，据此有可能用Z曲线识别这些位置。由此图可见，用Z曲线这种几何方法显示DNA序列不仅直观，而且作为一种识别序列中的不同基因和功能区的新方法，展现了广阔的应用前景。第二节第二节 序列联配序列联配一、一、Needleman-Wunsch算法算法二、二、Smith-Waterman算法算法三、序列相似性统计特征三、序列相似性统计特征四、替换矩阵四、替换矩阵五、多序列联配五、多序列联配一、一、Needleman-Wunsch算法算法有2种经典方法可以计算两条序列间的最适联配。
15、Needleman-Wunsch算法是一种整体联配(global alignment)算法，最佳联配中包括了全部的最短匹配序列。Smith-Wateman算法是在Needleman-Wunsch算法基础上发展而来的，它是一种局部联配(Local alignment)算法。这二种算法均可以用于核酸和蛋白质序列。在给定空位罚值和替换矩阵情况下，它们总是能给出具有最高（优）联配值的联配。但是，这个联配并不需要达到生物学意义上的显著水平。GCG软件包中，BESFIT和GAP程序，EMBOSS的needle等可用于该联配。一些网站可以通过递交序列进行两条序列的联配分析。 从整体上分析两个序列的关系，即考
16、虑序列总长的整体比较，用类似于使整体相似 (global similarity)最大化的方式，对序列进行联配。两个不等长度序列的联配分析必需考虑在一个序列中圈掉一些碱基或在另一序列作空位(gap)处理。Needleman 和Wunsch(1970)的法则为这些步骤提供了实例。这一算法是为氨基酸序列发展的，但也可以用于核苷酸序列。算法最初寻求的是使两条序列间的距离最小。尽管这类距离的元素是以一种特定的方式定义的，但该算法的良好特性在于它确定了最短距离。这是一个动态规划(dynamic programming)的方法。将两条联配的序列沿双向表的轴放置，两条序列的所有可能的联配方式都将在它们所形成的
17、方形图中（见下图）。从任一碱基对，即表中的任一单元开始，联配可延三种可能的方式延伸：如果碱基不匹配，则每一序列加上一个碱基，并给其增加一个规定的距离权重；或在一个序列中增加一个碱基而在另一序列中增加一个空位或反之亦然。引入一个空位时也将增加一个规定的距离权重。因此，表中的一个单元可以从(至多)三个相邻的单元达到。我们把达左上角单元距离最小的方向看作相似序列延伸的方向。等距离时意味着存在两种可能的方向。将这些方向记录下来，并在研究了所有的单元之后，沿着记录的方向就有一条路径可从右下角(两个序列的末端)追踪到左上角 (两个序列的起点)。由此所产生的路径将给出具有最短距离的序列联配。以两个短序列CT
18、GTATC和CTATAATCCC为例，将上述过程说明于图3.4。设碱基错配时距离权重为1，引入一个空位时距离权重为3。该图边缘的行和列作为起始条件增加到表中。在单元5行3列，即相应较短序列(第二序列)的第2个T碱基和较长序列(第一序列)的第1个T碱基位置，有三种可能的距离增量。设在各序列中增加碱基T时 (从4行2列移动)对距离的贡献为0。从5行2列的位置作水平移动(等价于增加第二序列的碱基T而在第一序列引入一个空位)，在本例中增加一个罚值3。从3列4行向该单元作垂直移动，使第一序列增加碱基T而第二序列引入一个空位，结果也得到一个罚值3。因此从该单元(5行3列)所得到的最小距离的延伸方向是沿对角
19、线和水平方向。在表中这两个方向用箭头表示。这两种最短方向都使从左上角到该单元的距离为6。在上述6种联配中，距离均为10，即在较短序列中有6个匹配碱基、1个错配碱基和3个空位。 沿箭头所指方向在表中从右下角向左上角追踪，得到6种可能的联配：该算法可以用代数形式来描述。设具有碱基ai和bj的两个序列a和b，这两个序列间距离为d(a,b)。通过评价序列a中前i个位置和序列b前j位置的距离，递归地得到距离d(ai,bj)。如果a和b的长度为m和n，则其期望距离为d(am,bn) 。上表中引入的第1行1列单元的距离为0(相当于空序列)，在单元(i,j)内，使到达该单元距离增加的三种可能事件为： 1.从单
20、元(i-1,j)向(i,j)的垂直移动，相当于在b序列中插入一个空位使相似序列延伸。换言之，b序列由a序列中ai的缺失所产生，这一事件的权重记作w_(ai)。2.从单元(i-1,j-1)向(i,j)的对角线移动，相当于增加碱基ai和bj使相似序列延伸。换言之，b序列由a序列中的ai被bj取代所产生，这一事件的权重记为w_(ai,bj)。3.从单元(i,j-1)向(i,j)的水平移动，相当于在序列b中插入一个空位使相似序列延伸。换言之，b序列由bj插入a序列所产生，这一事件的权重记为w +(bj)。因此，单元(i,j)的距离可看成三个相邻单元的距离加上相应权重后的最小者，即方程3.3且初始条件为
21、在图3.4的实例中 当两个序列被联配时，通过计算其重排序列(shuffed version)的联配距离，可以得到这两个序列间的最小距离估计。如果实际得到的联配距离小于重排序列距离的95%，则表明实际的联配距离达到了5%的显著水平，是不可能由机误造成的。 二、二、Smith-Waterman算法算法由于亲缘关系较远的蛋白质序列可能只有一些相互独立的相同片段，所以进行局部相似性分析有时可能比整体相似性分析更合理。Smith和Waterman描述了一种查找具有最高相似性片段的算法。对于序列A=(a1,a2,am)和 B=(b1,b2,bn)，Hij被定义为以ai和bj碱基对结束的片段(亚序列)的相似
22、性值。与Needle-Wunsch算法一样，Smith-Waterman算法也要利用递推关系来确定H值，H的初始值为： 相似性计算中包括2个统计量：碱基对(序列因子) ai,bj的相似性值S(ai,bj)和空位权重wk=v+uk(k 为空位长度)。Smith-Waterman算法可以给出2条序列的最大相似性值。以ai,bj碱基对结束的片段可以由以ai-1和bj-1结束片段增加碱基(因子)来获得，或者ai可以删除k长度的碱基片段，bj可删除l长度碱基片段。具体算法如下： 方程3.4则方程3.5其中该算法可以确保具有最大Hij值的序列片段是相似性最好的。从(ai,bj)为起点，向后 追踪矩阵，直到
23、到达某一负值。对于具有最大相似性片段以外部分的差异性不会影响到该片段的H值。 举例说明了这一算法。我们同样以上节Needleman-Wunsch算法中的两条短序列为例。两条序列(CTGTATC和CTATAATCCC)排于表3.9的两侧，相应的 和 值分别列入表中。本例的权重等根据Smith和Waterman(1981)以前的例子设定为： 方程3.6对于4个碱基具有相同频率的随机长序列，S(ai,bj)值的平均值为零。wk值应至少不小于匹配与不匹配权重的差值。 表3.9的最大Hij为4.33(8行与7列相交处)，星号(*)表示出具有最大相似性的片段匹配方式： 三、序列相似性统计特征三、序列相似性
24、统计特征到目前为止，对局部联配的统计学问题已基本搞清楚，特别是那些不含有空位(gap)的局部联配更是如此。我们不妨首先考虑不含有空位的局部联配问题，BLAST最初的搜索程序便是以此为基础的。 无空位局部联配涉及的是等长度的一对序列片段，两个片段的各部分彼此比较。一种Smith-Waterman或Sellers算法的改进算法可以找到所有高比值片段对(high-scoring segment pairs,HSPs)，即这些片段对的比较分值不会因片段的延伸而进一步升高。 为了分析上述分值随机性产生的几率大小，需要建立一个随机序列模型。对于蛋白质而言，最简单的序列模型可通过从一条序列中随机地选取氨基酸
25、残基，当然这一条序列中各种残基的频率必需一定。另外，一对随机氨基酸的联配期望值必需为负值，否则不论联配片段是否相关的，都会得到高比值，统计理论也将派不上用场。 就象独立随机变量之和总是倾向于正态分布(normal distribution)一样，独立随机变量的最大值倾向于极值分布(extreme value distribution)。在研究最佳局部联配时，主要涉及的是后一种情况。在一定的序列长度m和n限定下，HSP的统计值可由2个参数(k和)确定。最简单的形式，即不小于比较值为S的HSP个数，可由下列公式算得其期望值： 我们称该期望值为比值S的E值(E-Value)。方程3.7上述公式非常灵
26、敏。在给定比值的情况下，将比较序列长度加倍，则HSP数(即E值)也将加倍，同样，S值为2X的某个HSP长度必是S值为X的两倍，所以E值将随着s值的增大急剧减少。参数K和可分别被简单地视为搜索步长(search spacesize)和计分系统(scoring system)的特征数。 1、二进制值或标准比值、二进制值或标准比值(Bit score) 最初获得的比值(S)在没有计分系统或统计量K和的辅助下，没有什么意义。单独的比值就如同没有单位(米或者光年)的距离。可使比值按下式标准化： 方程3.8获得S 值就如同得到了具有标准单位的数值。 E值因此可简化为： 方程3.9二进制值使所使用的计分系统赋予了统计学意义，使除了可以确定搜索步长外，同样可以计算相应的显著水平。
