生信课程笔记1-序列算法

这篇笔记是我蹭课《生物信息学》的课堂笔记,本人为非生信专业,主要做实验,所以此文仅供个人学习。


序列搜索算法的进化

虽然现代高通量生物分子技术产生了各种类型的数据,但生物序列数据(biosequence data)仍然是生物信息学分析的核心。

NCBI在2007年底推出了SRA(Sequence Read Archive)数据库,用于存储二代测序数据,目前包含大约15PB碱基,其中几乎一半在开放访问中。SRA数据库地址https://www.ncbi.nlm.nih.gov/sra/。SRA是国际核苷酸测序数据库合作(INSDC)项目的一部分,INSDC包括NCBI的SRA,欧洲生物信息研究所(EBI)和日本DNA数据库(DDBJ),SRA数据库是美国国立卫生研究院(NIH)存储高通量测序数据的主要数据库。数据上传到以上任何一个组织,其他两个组织都是可以共享的。

序列比对的意义:相似的序列可能会有相似的结构,从而有相似的功能。序列之间的相似性可以帮助我们推断一个未知新序列可能的功能。不同物种中相似的序列往往意味着其具有共同的祖先,即同源。事实上,序列间的相似性是在演化分析中用来构建演化树的重要依据之一。(序列同源是指两条序列有共同的祖先,然后经过突变变得不同。同源无法考证,只能通过相似性来推断同源。同源不一定相似,相似也不一定同源。)

相似性矩阵(similarity matrix)或打分矩阵:两条序列做全局比对,然后计算全局比对中一致字符的个数和相似字符的个数,再除以全局比对的长度,就可以得到两序列的一致度(identity)和相似度(similarity)。替换记分矩阵是反映残基之间相互替换率的矩阵,也就是说它描述了残基两两相似的量化关系。DNA的替换计分矩阵有等价矩阵(相同1分,不同0分)、转换-颠换矩阵(转换-1分,颠换-5分)、BLAST矩阵(相同5分,不同-4分,是实践经验所得)。蛋白质的常用替换计分矩阵PAM 矩阵基于进化原理,如果两种氨基酸替换频繁,说明自然界容易接受这种替换,那么这一对氨基酸替换的得分就高。基础的PAM-1矩阵反应的是进化产生的每一百个氨基酸平均发生一个突变的量值,由统计方法得到。PAM-1自乘n次,可以得到PAM-n,表示两序列的亲缘关系更远。目前最常用的蛋白序列比对打分矩阵是BLOSUM 62矩阵,这个矩阵是由一致度≥62%的序列计算而来的。≥0 的得分代表对应的一对氨基酸为相似氨基酸,<0 的是不相似的氨基酸。

早期算法:动态规划(dynamic programing

检测两条短序列(如基因)的相似性(similarities),通常只考虑点突变(point mutations),即字符替换(character substitutions)、插入或删除(insertions or deletions,indel)。最简单的两序列比对方法有打点法,用对角线表示相同区域。序列比对就是运用特定的算法找出两个或者多个序列之间产生最大相似度。根据比对序列的个数可以把序列比对分为双序列比对和多序列比对。根据序列比对的算法不同,双序列比对又分为全局比对和局部比对。

编辑距离(edit distance)是指两个字符串之间,由一个转成另一个所需的最少编辑操作次数。可进行替换、插入或删除操作。一般来说,编辑距离越小,两个字符串的相似度越大。简单的编辑距离不足以反映生物的距离,被加权编辑距离(weighted edit distance)所取代,通过相似性评分(similarity score)来定义,其中不同的字符替换可能有不同的罚分(penalties)。然后将两个序列的整体比较形式化为最优全局比对(optimal global alignment)问题,并用生物信息学中的经典动态规划算法Needleman-Wunsch算法求解。

动态规划是求解最优化问题的一种途径,其基本思想是将过程分成若干个互相联系的阶段,在它的每一阶段都作出最优决策,从而使整个过程达到最好的活动效果。动态规划算法对于序列长度m、n具有时间复杂度(time complexity)O(nm),因此对于大规模序列搜索的计算要求太高。虽然全局比对的空间复杂度(space complexity)可以通过一个简洁的分而治之的技巧(divide-and-conquer trick)使之线性化。

最优全局比对算法——Needleman-Wunsch算法:使用先前对于较小子序列的最优比对结果,来进行下一步的最优比对。按照线性空位罚分(linear gap penalty)(1个gap -2分,k个gap -2k分)计算序列比对的得分。

举个栗子:右下图比对,绘制多一行多一列的表格,从左上角0开始,向下移动(在横向序列引入一个gap)得-2分,向右移动(在纵向序列引入一个gap)得-2分,对角线移动得到match得1分,mismatch得-1分。每个位置得分由它的上一个位置得分加上移动得分计算得到。最终右下角最高得分为S(i,j)=-1,有3种最优全局比对方式(从右下角到左上角回溯的路径)。

最优局部序列比对(Local sequence alignment)算法——Smith-Waterman算法:寻找输入的两个序列V和M的子序列(subsequences)的最优全局比对。适用于寻找相似的小片段,如蛋白质的共有结构域、不同序列的保守区域。计算方法与全局比对相似,只是整个表格中没有负分,最小为0分。在整个计分表中找到最大的分数,并向前追溯到0停止。

Gotoh’s算法进一步将Needleman-Wunsch算法和Smith-Waterman算法扩展到仿射空位罚分(affine gap penalties)的情况,对k个连续(consecutive)插入或删除的罚分是空位开始罚分g(gap opening penalty)与空位延伸罚分e(k gap extension penalties)之和,penalty = g + (k − 1)e,|e|≤|g|,往往第一个空位罚分较多。gap开头罚分大,gap延长罚分小的时候,gap集中连成长串出现。

长度分别为m和n的两条序列的所有可能的比对个数为

Filtration-based heuristics and database search 基于过滤的试探法和数据库搜索

基本局部比对搜索工具BLAST(Basic Local Alignment Search Tool)是一种启发式的(heuristic)序列比对方法,它遵循两种一般的思想:过滤和索引(filtration and indexing)。与给出最优解的动态规划算法不同,BLAST是一种近似算法。

BLASTp是用蛋白质序列搜索蛋白质序列数据库,BLASTn是用核酸序列搜索核酸序列数据库,BLASTx是将核酸序列按6条链翻译成蛋白质序列后搜索蛋白质序列数据库。tBLASTn 是用蛋白质序列搜核酸序列数据库,核酸数据库中的核酸序列要按6条链翻译成蛋白质序列后再被搜索。tBLASTx是将核酸序列按6条链翻译成蛋白质序列后搜索核酸序列数据库,核酸数据库中的所有核酸序列也要按6条链翻译成的蛋白质序列后再被搜索。除了按照搜索内容分类,BLAST还可以根据搜索算法不同分为标准BLAST,PSI-BLAST(位点特异性迭代,Position-Specific Iterated)和PHI-BLAST(模式识别,Pattern-Hit Initiated)等。

BLAST的思想是真正匹配的比对很有可能包含得分非常高的配对片段,即seeds。片段对是指两个给定序列中的一对子序列,它们的长度相等,且可以形成无空位的完全匹配。首先在搜索的序列中选取固定长度的小片段(如3个氨基酸或11个核苷酸),称为seeds/anchors,也就是在序列中从头到尾一个残基一个残基地移动,每3个作为一个seeds。然后从打分矩阵计算各种3氨基酸组合的得分,设置一个阈值,大于阈值的3氨基酸作为比对上的,数据库中不能匹配的部分都从考虑中丢弃(过滤filtered out)。然后从seeds开始向两端扩展,得到高分值片段对(high scoring segment pairs,HSPs)。注意,低复杂度的序列(重复序列)应该首先过滤掉。这样计算复杂度是序列长度n的一次方,如果做双序列比对话需要构建一个n乘以n的表格,计算复杂度是n的二次方。

搜索使用数据库的预构建索引(pre-built index),通常是哈希表(hash table),加快定位seed patterns的位置。这种方法的性能被定义为灵敏性(Sensitivity)与特异性(specificity)之间的权衡(trade-off),并由用于过滤的种子类型决定。由于分数小于阈值的seeds被直接过滤掉了,所以其灵敏性不高。灵敏性是指真实为阳性的样本中,用这种方法能够检测出来的比率。

由于连续11个核苷酸完全相同(k consecutive matching nucleotides)要求太高,又提出了一种改进的算法,选择根据给定模式的k个间隔的匹配核苷酸作为种子,即间隔的种子(spaced seed),可以改善灵敏性与特异性之间的权衡。

索引是为了加速对表中数据行的检索而创建的一种分散的存储结构。索引是针对表而建立的,它是由数据页面以外的索引页面组成的,每个索引页面中的行都会含有逻辑指针,以便加速检索物理数据。哈希表(Hash table,也叫散列表),是根据关键码值(Key value)而直接进行访问的数据结构,通过把关键码值映射到表中一个位置来访问记录,以加快查找的速度,这个映射函数叫做哈希函数或散列函数。Hash(散列、哈希),就是把任意长度的输入通过算法变成固定长度的输出。直观解释起来,就是对一串数据m进行杂糅,输出另一段固定长度的数据h,作为这段数据的特征(指纹)。MD5和SHA都是历史悠久的Hash算法。


(参考文献:Kucherov,G. (2019) Evolution of biosequence search algorithms: a brief survey. Bioinformatics, 35, 3547–3552.)

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 156,907评论 4 360
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 66,546评论 1 289
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 106,705评论 0 238
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 43,624评论 0 203
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 51,940评论 3 285
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 40,371评论 1 210
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,672评论 2 310
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 30,396评论 0 195
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 34,069评论 1 238
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 30,350评论 2 242
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 31,876评论 1 256
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 28,243评论 2 251
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 32,847评论 3 231
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 26,004评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,755评论 0 192
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 35,378评论 2 269
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 35,266评论 2 259