生信了

本文介绍了BWT算法。bwa是目前最流行的二代测序比对工具，其中就用到了BWT算法。BWT（Burrows-Wheeler Transform）算法是一种数据转换算法，它将一个字符串中的相似字符放在相邻的位置，以便于后续的压缩。简要回顾BWT算法可以分为编码和解码两部分。编码后，原始字符串中的相似字符会处在比较相邻的位置；解码就是将编码后的字符串重新恢复成原始字符串的过程。BWT的一个特...

前文已经介绍过KMP算法和Trie树，本文将在此基础上介绍AC自动机。之前的序列比对文章大都在利用动态规划算法解决字符串的非精准匹配（允许错配、插入和缺失），比如全局比对和局部比对问题。当然，后来我们还介绍了模序发现和中间字符串问题，并初次学习了如何用分支定界法解决这一类问题。既然有非精准匹配问题，就有精准匹配问题。所谓精准匹配，就是两个字符串在比对时，不允许错配、插入和缺失。其实，笔者之...

本文介绍两个字符串的编辑距离并给出代码。编辑距离所谓编辑距离，就是给定两个字符串后，将一个字符串变为另一个字符串所需要花费的最少步骤。这个改变包括“插入一个字符”、“删除一个字符”，“替换一个字符”。比如：v=TGCATAT\bm{v}=TGCATATv=TGCATAT与w=ATCCGAT\bm{w}=ATCCGATw=ATCCGAT这两个字符串的编辑距离为4。编辑距离的求解过程和全局比...

本文介绍如何求解两个字符串的最长公共子序列。最长公共子序列问题前文《序列比对（23）最长公共子字符串》介绍了如何求解两个字符串的最长公共子字符串，本文将介绍如何求解两个字符串的最长公共子序列。二者听起来很像，所以我们首先得说明一下子字符串和子序列的区别。在本文语境中，子字符串是由原字符串中的连续字符组成；而子序列是从原字符串中选择字符，只需要保持其次序不变（可以说，子字符串都是子序列，但...

本文介绍如何求解两个字符串的最长公共子字符串。其实这个问题可以放在序列比对专题的最开始，只是笔者是个新手，所以当初只是照《生物序列分析》教材的进度写的，教材是直接从全局比对开始讲的。Anyway，我们在本文介绍也不迟。给定两个字符串v\bm{v}v和w\bm{w}w，长度分别为mmm和nnn，如何找到这两个字符串的最长公共子字符串呢？所谓最长公共子字符串，顾名思义，很好理解。举例来说：如果...

前文介绍了中间字符串的算法和代码，但是使用分支定界策略时所使用的界限是很宽松的。本文给出了一个更紧的界限。对分支定界法的简单回顾前文《序列比对（21）中间字符串问题的算法及实现代码》介绍了中间字符串的算法和代码，但是使用分支定界策略时所使用的界限是很宽松的。分支定界法的伪代码如下：引自《生物信息学算法导论》其中关键的是这几句：optimisticDistance←TotalDist...

前文介绍了基序发现问题和中间字符串问题。本文给出了中间字符串的算法和实现代码。中间字符串问题的简单算法及伪代码前文《序列比对（19）基序发现和中间字符串问题》介绍了基序发现问题和中间字符串问题；《序列比对（20）基序发现问题的算法及实现代码》给出了基序问题的算法和实现代码。本文将介绍中间字符串问题的算法，并给出实现代码。简单回顾一下，中间字符串问题其实就是要找到使得总距离最小的一个lll...

前文介绍了基序发现问题和中间字符串问题，本文给出了基序发现问题的具体算法和实现代码。基序发现问题的简单算法及伪代码前文《序列比对（19）基序发现和中间字符串问题》介绍了基序发现问题和中间字符串问题，本文将介绍基序发现问题的算法，并给出实现代码。简单回顾一下，基序发现问题其实就是要找到使得共有序列得分最大的一组起始位点。argmaxs⃗ Score(s⃗,DN...

本文介绍了基序发现问题和中间字符串问题。引言：DNA调控元件我们知道，DNA调控元件往往是一段相似的DNA序列。理想情况下这些序列完全一致，比如下面这样：图片引自《生物信息学算法导论》但实际上，这些序列不会完全一样，总会有若干位点发生“变异”，从而不同，比如下面这样：图片引自《生物信息学算法导论》如果给定一组DNA序列（暂且假定它们长度相等），那么如何找出这些相似的序列呢？由此...

前文介绍了重复匹配问题的动态规划算法，但是遗留了重复结果输出的问题。本文对该问题进行了补充说明。前文《序列匹配(五)——重复匹配问题的动态规划算法》介绍了重复匹配问题的动态规划算法。重复匹配问题就是从序列xxx中找到序列yyy的多个拷贝，并且使这多个拷贝的“标准比对分值”之和最大。《生物序列分析》一书中给出的迭代公式是：F(i,0)=max⁡{F(i−1,0),F(i−1,j)−T,&nb...

序列比对的系列文章第二部分主要介绍了HMM（隐马尔科夫模型），包含了八篇文章：《序列比对（九）从掷骰子说起HMM》《序列比对（十）viterbi算法求解最可能路径》《序列比对（11）计算符号序列的全概率》《序列比对（12）：计算后验概率》《序列比对（13）后验解码》《序列比对（14）viterbi算法和后验解码的比较》《序列比对（15）EM算法以及Baum-Welch算法的推导》《...

原创：hxj7本文介绍了如何用Baum-Welch算法来估算HMM模型中的概率参数。Baum-Welch算法应用于HMM的效果前文《序列比对（15）EM算法以及Baum-Welch算法的推导》介绍了EM算法和Baum-Welch算法的推导过程。Baum-Welch算法是EM算法的一个特例，用来估算HMM模型中的概率参数。其具体步骤如下： 图片引自《生物序列分析》本文给出了Baum-...

原创：hxj7本文介绍了 EM算法 和 Baum-Welch算法的推导过程。一般地，估算概率模型参数可以用最大似然法，即找到θ^=argmaxθ P(x∣θ)\hat{\theta} = \mathrm{argmax}_\theta \, P(x|\theta)θ^=argmaxθ​P(x∣θ)有时为了简便运算，也可以用最大对数似然代替最大似然，即以下公式...

本文比较了viterbi算法求解最可能路径以及后验解码这两种不同的解码方法。前文《序列比对（十）viterbi算法求解最可能路径》介绍了用viterbi算法求解最可能路径：在符号序列已知而状态序列未知时，最可能路径是：本文将这两种方法比较了以下，看它们各自求解的路径差异是否显著。分两种情况：一、如前面几篇文章一样，从公平骰子转为作弊骰子的概率是0.05。效果如下：（其中Rolls一行...

原创: hxj7本文介绍了如何利用后验概率进行解码，可称为后验解码。前文《序列比对（12）计算后验概率》介绍了如何计算某一位置可能状态的后验概率。那么可以据此找到某一位置最有可能的状态。即从上面的公式可以看出，有两种方法求解某一位置的最可能状态。如果是依据公式(1)，先计算出后验概率，然后找到其中最大后验概率对应的状态；如果是依据公式(3)，无需计算后验概率，比较简单。本文所用代码...

原创: hxj7本文介绍如何计算状态的后验概率。前文《序列比对(十一)——计算符号序列的全概率》介绍了如何使用前向算法和后向算法计算符号序列的全概率。但是很多情况下我们也想了解在整条符号序列已知的情况下，某一位置符号所对应的状态的概率。也就是说要计算的概率。很明显，此概率为一后验概率。要计算上述后验概率，可以经过以下推导：其中：根据公式(1),(4),(5),(6)，可以重新...

原创：hxj7前文介绍了在知道符号序列后用viterbi算法求解最可能路径。本文介绍了如何使用前向算法和后向算法计算符号序列的全概率。如果一个符号序列中每个符号所对应的状态是已知的，那么这个符号序列出现的概率是容易计算的：但是，如果一个符号序列中每个符号所对应的状态未知时，该怎么求取这条序列的概率呢？我们知道：如果我们用穷举法求出所有的P(x,π)是不现实的，因为随着序列长度的增长...

原创: hxj7本文介绍了如何使用viterbi算法得到概率最大的路径。前文《序列比对（九）从掷骰子说起HMM》已经通过投骰子的例子介绍了HMM的基本概念，引用如下：那么如何通过已知的符号序列来预测未知的状态序列呢？首先想到的自然是最大似然法，即看不同的状态序列下的概率，选取使整体概率最大的那个状态序列。如果将一串状态序列称为一条路径的话，我们感兴趣的就是找到一条最大似然的路径。如《...

从本文开始，到了序列比对的第二部分，主要是介绍HMM以及其在序列比对中的基础应用。前面八篇文章介绍了动态规划在序列比对中的基础应用。从本文开始，开始介绍HMM（隐马尔可夫模型）。###Markov链首先我们先介绍马尔可夫链（Markov链），这个大家可能都熟悉，高中数学中就学过。简单来说，Markov链就是一系列状态构成，每一个状态出现的概率只与它前一个状态有关。具体地说，假设有一系列...

原创: hxj7不知不觉中，已写了七篇关于序列比对的文章：《序列比对（一）全局比对Needleman-Wunsch算法》《序列比对（二）Needleman-Wunsch算法之仿射罚分》《序列比对（三）局部联配Smith-Waterman算法》《序列比对（四）Smith-Waterman算法之仿射罚分》《序列匹配（五）重复匹配问题的动态规划算法》《序列比对（六）交叉匹配问题》以及《序...

原创: hxj7一般而言，运用动态规划算法进行序列比对对内存空间的要求是 O(mn) 阶的，本文介绍了一种线性空间要求的序列比对方法。前文如《序列比对（一）全局比对Needleman-Wunsch算法》所介绍的运用动态规划算法进行序列比对时，对内存空间的要求是 O(mn) 阶的。如果只是要求最大得分而不要求最大得分所对应的序列（也就是不进行回溯）的话，其实只需要保留上一行的得分就够了，这一...

原创: hxj7之前几篇文章介绍了全局匹配以及局部匹配，本文介绍交叉匹配问题并给出代码。交叉匹配所谓交叉匹配（overlap alignment 或者叫 glocal alignment），就是两条序列中至少有一条的头部序列要参加比对并且至少有一条的尾部序列要参加比对。一般而言，就是下面两种情形：一种是两条序列有重叠的部分，但互不包含。比如x序列的头部与y序列的尾部匹配。第二种是...

原创： hxj7前言：蛋白质序列中常有重复的功能域（domain）或模体（motif）拷贝，由此衍生出一个抽象的序列多重匹配的问题，即如何从一个序列中找出另一个序列的某部分（如功能域或模体）的多个无交叠（non-overlapping)拷贝。本文给出了该问题的示例、关键计算公式以及C语言实现代码。问题及算法描述更具体地描述上面的问题：有序列x和y，其中y是包含结构域的序列，x是要从中找...

原创: hxj7前言： 本文介绍的是采用仿射罚分模型的Smith-Waterman算法。关于全局联配，局部联配以及仿射罚分模型的介绍可参见前文：序列比对（一）全局比对Needleman-Wunsch算法序列比对（二）Needleman-Wunsch算法之仿射罚分序列比对（三）局部联配Smith-Waterman算法其中，《序列比对（二）Needleman-Wunsch算法之仿射罚分...

原创：hxj7关于全局联配的介绍可参见前文：序列比对（一）全局比对Needleman-Wunsch算法序列比对（二）Needleman-Wunsch算法之仿射罚分所谓局部联配，就是两条序列的子序列的联配。局部联配算法就是找到联配得分最高的子序列。其中最常见的就是Smith-Waterman算法。Smith-Waterman算法与Needleman-Wunsch算法类似，只是在计算得分...

原创： hxj7关于全局比对Needleman-Wunsch算法的介绍及线性罚分的实现代码请参考序列比对（一）全局比对Needleman-Wunsch算法。概念引用如下：线性罚分，即penalty＝g*d，其中g为空位长度，d为一个空位的罚分。仿射型罚分，即penalty=d+(g-1)*e， 其中g为连续空位的长度，d为连续空位中第一个空位的罚分，e为连续空位中第2个及以后空位的罚分。...

原创：hxj7前言序列比对是生信领域的一个古老课题，在这一波NGS的浪潮中重新引起大家的广泛关注。由于生物序列的特殊性，在比对的时候允许插入缺失，所以往往是一种不精确匹配。从本文开始，我们陆续学习前人开发出的流行算法。全局比对算法所谓全局比对算法，就是根据一个打分矩阵（替换矩阵）计算出两个序列比对最高得分的算法。关于它的介绍网上已经非常多了，我们只需看看其中的关键点及实现代码。关键点1...

关键词：trie; prefix; search; match;字典树，又称单词查找树，是一个典型的一对多的字符串匹配算法。“一”指的是一个模式串，“多”指的是多个模板串。字典树经常被用来统计、排序和保存大量的字符串。它利用字符串的公共前缀来减少查询时间，最大限度地减少无谓的字符串比较。那它一般应用在什么地方呢？我们举一个例子说明：假设有一个单词表，里面有10w个单词。如果别人给你2000个...

关键词：KMP; string; match;字符串匹配是一个既古老又现代的问题，历久弥新。生信领域中字符串处理更是daily work。诸如bwa这般神一样的软件，本质上也是在解决字符串非精准匹配的问题。所以，从本文开始，我们陆续会分享一些对我们有用的字符串匹配算法。目前计划是先讲三个算法：KMP算法字典树算法AC自动机算法今天介绍的是KMP算法，一种常用的字符串精准匹配算法。这...

原创：hxj7关键词：k-mer; recursive; trick;什么是k-mer？比如，“ATGC”的所有1-mer是：’A’, ‘T’, ‘G’, ‘C’。共4^1=4种组合。而“ATGC”的所有2-mer是“AA”, “AT”, “AG”, “AC”“TA”, “TT”, “TG”,“TC”“GA”, “GT”, “GG”,“GC”“CA”, “CT”, “CG”,“CC...

原创：hxj7关键词：蓄水池算法;fastq文件往往都很大，出于测试目的，我们经常要从fastq文件中随机抽取reads，生成一个小一点的fastq文件，以加快测试效率。假设我们要从一个包含大约100M reads的fastq文件中随机抽取1M reads，该怎么办呢？我们将问题简单化：假设我们要从一个txt文件中（不知道总共多少行）随机抽取M行（fastq文件的处理与之类似，只不过fast...

原创：hxj7关键词：phred; trim; mott;NGS（二代测序）分析的起点往往是fastq文件。fastq文件其实就是一条条的记录，每个记录包含4行。其中比较重要的是第二行和第四行：第二行是测序得到的碱基序列，第四行是每个碱基相应的测序质量，测序质量越高代表该碱基被测错的概率越低，反之越高。正因为二代测序是有一定的错误率的，所以我们在进行下游分析之前，常常要对fastq文件中的r...