基因组序列比对算法介绍(一)

基因组重测序中序列比对介绍

重测序基因组数据比对,是指将测序仪下机fastq数据(NGS read序列,通常100-150bp),与人类参考基因组(reference)进行匹配,允许错配(mismatch),插入缺失(indel),目的是在参考基因组找到序列最相似的位置,通常是基因组分析(包括 variation calling,ChIP-seq,RNA-seq,BS-seq)流程的第一步。

常用算法

image

图一

汉明距离(Hamming distance)表示两个(相同长度)字对应位置不同的数量,我们以d(x,y)表示两个字x,y之间的汉明距离。对两个字符串进行异或运算,并统计结果为1的个数,那么这个数就是汉明距离。图中read1最佳位置的方法,就是通过查找最小汉明距离的实现的。

编辑距离(Edit distance)是针对二个字符串(例如英文字)的差异程度的量化量测,量测方式是看至少需要多少次的处理才能将一个字符串变成另一个字符串。图中read3最佳位置,通过查找最小编辑距离的方法实现。

image

图二

全局比对(Global alignment):全局比对是指将参与比对的两条序列里面的所有字符进行比对。全局比对在全局范围内对两条序列进行比对打分,找出最佳比对,主要被用来寻找关系密切的序列。其可以用来鉴别或证明新序列与已知序列家族的同源性,是进行分子进化分析的重要前提。其代表是Needleman-Wunsch算法。图一中,read3使用全部比对。

局部比对(Local alignment):与全局比对不同,局部比对不必对两个完整的序列进行比对,而是在每个序列中使用某些局部区域片段进行比对。其产生的需求在于、人们发现有的蛋白序列虽然在序列整体上表现出较大的差异性,但是在某些局部区域能独立的发挥相同的功能,序列相当保守。这时候依靠全局比对明显不能得到这些局部相似序列的。其次,在真核生物的基因中,内含子片段表现出了极大变异性,外显子区域却较为保守,这时候全局比对表现出了其局限性,无法找出这些局部相似性序列。其代表是Smith-Waterman局部比对算法。图一中,read2使用局部比对。

image

图三

Smith-Waterman算法介绍

Smith-Waterman是由Temple F. Smith和Michael S. Waterman于1981年提出的一种进行局部序列比对(相对于全局比对)的算法,用于找出两个核苷酸序列或蛋白质序列之间的相似区域。该算法的目的不是进行全序列的比对,而是找出两个序列中具有高相似度的片段。S-W算法基于动态规划,它接受任意长度、任意位置、任意序列的对齐,并确定是否能找到最优的比对。

简单地说就是,动态规划找到问题中较小部分的解,然后把它们放在一起,形成整个问题的一个完整的最优最终解。

它优于BLAST和FASTA算法,因为它搜索了更大的可能性,具有更高的敏感性。

S-W算法不是一次查看整个序列,而是对多个长度的片段进行比较,寻找能够最大化得分的片段。算法本身本质上是递归的:

image

图四

算法步骤如下:

  1. 确定置换矩阵及空位罚分方法。

    如果两个位点的碱基匹配(MATCH),得3分,出现不匹配(MISMATCH)为-3分。

    如果出现空缺(GAP)惩罚分数做一个线性增长惩罚:W(k) = 2k (k指的是出现GAP的次数,初始值为2,也就是每次出现一个GAP惩罚的分数依次为-2,-4,-6,递增)。

  2. 初始化得分矩阵。

    矩阵的行列同样是是两个序列的碱基排列,第一行和第一列为置0,这是和Needleman-Wunsch算法不同的地方。

    得分矩阵的长度和宽度分别为两序列的长度+1。其首行和首列所有元素均设为0。额外的首行和首列得以让一序列从另一序列的任意位置开始进行比对,分值为零使其不受罚分。

  3. 打分。对得分矩阵的每一元素进行从左到右、从上到下的打分,考虑匹配或错配(对角线得分),引入空位(水平或垂直得分)分别带来的结果,取最高值作为该元素的分值。如果分值低于0,则该元素分值为0。打分的同时记录下每一个分数的来源用来回溯。

    image
                    ![image](https://upload-images.jianshu.io/upload_images/26826111-bc131680108a8dae?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)
    
                                                        图五
    
  4. 回溯。通过动态规划的方法,从得分矩阵的最大分值的元素开始回溯直至分数为0的元素。具有局部最高相似性的片段在此过程中产生。具有第二高相似性的片段可以通过从最高相似性回溯过程之外的最高分位置开始回溯,即完成首次回溯之后,从首次回溯区域之外的最高分元素开始回溯,以得到第二个局部相似片段

image
                                                                            图六

基因组分析***** 微信公众号推出 《50篇文章深入理解NGS》系列文章, 第三篇文章 《基因组序列比对算法介绍(一)》,争取每周更新一篇高质量生信干货帖子。

关注 "基因组分析" 微信公众号,了解最新最全生信分析知识。

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

推荐阅读更多精彩内容