rHAT,国内首个三代序列比对算法

本文来自“生信算法”公众号。

        对于以Illumina为代表的二代测序技术,研究者们开发出了许多优秀的序列比对算法,如BLAST、Bowtie2 、BLAT等软件。随着测序技术的不断发展,尤其是近几年以单分子测序技术为代表的三代测序技术的出现,测得的序列长度长达十万碱基数,远远高于二代序列长度,且同时具有较高的测序误差(错误率~15%)。因此,绝大多数针对二代测序数据的序列比对工具不适合处理三代测序数据。所以,开发基于三代测序数居的序列比对算法尤为迫切。

不同序列比对算法
不同序列比对算法

        在上一篇文章中(三代测序序列比对利器-BLASR,更小更快更方便),算法哥介绍了BLASR三代测序比对算法,BLASR是2012年开发的首个三代序列比对算法,但其运行速度慢。最近又有好多新发表的三代测序比对算法。其中rHAT(rHAT: fast alignment of noisy long reads with regional hashing)作为国内第一篇三代序列比对算法,由哈尔滨工业大学的王亚东教授团队开发,16年发表在Bioinformatics期刊上。博哥作为一个算法开发者,自然需要仔细拜读。因此,本次简单介绍rHAT算法的主要思想,希望对大家有所启发,可以用rHAT来处理自己的序列。

        rHAT(Regional Hashing-based Alignment Tool )算法主要是采用了区域哈希(regional hashing)策略。算法主要分为两部分,首先是定位(在整个基因组中找到与待比对序列相似的区域),通过区域哈希将待比对序列定位到基因组中与其相似的区域,如下图左图所示。然后根据区域内的种子构建有向无环图,找到比对路径进行详细比对,如下图右图所示。

rHAT算法框架(左图定位,右图详细比对)

哈希表定位

        首先看一下如何通过区域哈希表在基因组序列中对序列进行定位。

在基因组中通过区域哈希表对序列定位

  1. 如上图a所示,首先将基因组分割成长度为L的窗口,其中相邻两个窗口重叠L/2部分,这就是区域哈希表中区域的意思,一个窗口就是一个区域。
  2. 然后构建区域哈希表,如图b所示,将基因组每个k-mer(长度为k的子片段,又称作种子)进行哈希编码,每个哈希表中存储的是出现这个k-mer的窗口的位置。
  3. 然后对于待比较序列,如图c所示,提取待比较序列的每个k-mer,通过图b中的哈希表就可以找到包含每个k-mer的窗口,然后对窗口中的种子个数,即匹配个数,进行排序,找出种子个数最大的窗口,即为与待比较序列相似性较高的区域,用于接下来的序列比对。

序列路径选择与比对

        通过上面步骤在基因组中选取待比对区域后,利用种子扩展策略对非种子区域进行比对。

  1. 首先根据种子构建有向无环图,通过打分策略选取比对路径,如上图中a的红色圆圈所示。
  2. 找出比对路径中的种子,对种子间的间隙采用经典的动态规划算法进行详细比对
  3. 最后将种子间的比对结果与种子拼接,合并得到最终的长序列比对结果。

结果比较

        rHAT算法与目前BLASR和BWA-MEM两种方法进行了比较,测试数据为人类基因组数据,如下图所示,可以看出rHAT方法在比对碱基数、比对序列条数以及运行时间均有改善。希望可以对序列分析者有所帮助。详细的比对结果大家可以参考其原文。

代码下载地址

https://github.com/HIT-Bioinformatics/rHAT

原文链接

https://academic.oup.com/bioinformatics/article/32/11/1625/1742681

参考文献:

Bo Liu, Dengfeng Guan, Mingxiang Teng, Yadong Wang; rHAT: fast alignment of noisy long reads with regional hashing, Bioinformatics, Volume 32, Issue 11, 1 June 2016, Pages 1625–1631

        

### 高斯近似消息传递算法概述 高斯近似消息传递(Gaussian Approximate Message Passing, GAMP)是一种用于求解大规模线性逆问题的有效算法。它通过利用贝叶斯推断框架,在稀疏信号恢复等领域表现出优异性能。以下是关于该算法的实现原理、应用场景以及伪代码的具体说明。 --- #### 实现原理 高斯近似消息传递的核心思想在于将复杂的概率推理分解为一系列简单的局部计算,并假设这些计算中的变量服从高斯分布。这种方法能够显著降低复杂度并提升效率[^1]。 具体而言,GAMP 的主要目标是对以下形式的观测模型进行估计: \[ y = \Phi x + n, \] 其中 \( y \in \mathbb{R}^M \) 是观测向量,\( \Phi \in \mathbb{R}^{M \times N} \) 是测量矩阵,\( x \in \mathbb{R}^N \) 是待估信号,而 \( n \sim \mathcal{N}(0,\sigma_n^2I_M) \) 表示加性白噪声[^3]。 为了简化上述问题,GAMP 利用了因子图结构来进行消息传递操作。每一步的消息更新可以分为两部分:前向传播和反向传播。在每一层节点上,均值和方差被作为核心统计量进行交换,从而逐步逼近真实信号的概率密度函数。 --- #### 应用场景 由于其高效性和灵活性,GAMP 广泛应用于多个领域,主要包括但不限于以下几个方面: - **压缩感知**:当面对欠定线性系统时,可以通过先验知识重建原始稀疏信号。 - **通信信道均衡**:解决多径效应引起的干扰问题,改善数据传输质量。 - **图像处理与计算机视觉**:例如去噪、超分辨率重建等任务中也有出色表现[^2]。 --- #### 伪代码 下面是基于 Python 的简单版本 GAMP 算法伪代码描述: ```python def gamp_algorithm(y, Phi, A_prior_func, B_prior_func, max_iter=100, tol=1e-6): """ Gaussian Approximate Message Passing Algorithm Parameters: y (numpy array): Measurement vector of size M. Phi (numpy matrix): Sensing matrix of shape (M,N). A_prior_func (function): Function returning prior mean given scalar input rhat_jt. B_prior_func (function): Function returning variance scaling factor beta_jt from rhat_jt. max_iter (int): Maximum number of iterations allowed. tol (float): Convergence tolerance level. Returns: Estimated signal x_hat after convergence or reaching maximum iteration limit. """ import numpy as np # Initialization steps M, N = Phi.shape tau_v = 1 / np.diag(Phi.T @ Phi).mean() # Initial noise precision estimate v_tilde = np.zeros(N) u_bar = np.linalg.pinv(Phi) @ y x_hat_prev = np.random.randn(N) for t in range(max_iter): # Compute forward messages at variable nodes r_hat = u_bar - Phi.T @ ((v_tilde * (Phi @ x_hat_prev))[:, None]).squeeze() # Update posterior statistics using priors a_jt = A_prior_func(r_hat) b_jt = B_prior_func(r_hat) # Backward message computation at check node side zeta_t = 1 / (tau_v + sum(b_jt)) eta_t = zeta_t * (sum(a_jt) + tau_v * (y - Phi @ a_jt)) # Estimate new solution and compute residual error term x_hat_new = a_jt + b_jt * eta_t delta_x = np.abs(x_hat_new - x_hat_prev).max() if delta_x < tol: break # Prepare next round variables v_tilde = b_jt * zeta_t u_bar += Phi.T @ (zeta_t * (Phi @ (a_jt - x_hat_new))) x_hat_prev = x_hat_new.copy() return x_hat_new ``` 此代码片段展示了如何定义一个通用型 GAMP 函数,适用于不同类型的先验分布情况下的信号重构过程。 --- #### 示例 考虑这样一个例子——我们试图从一组随机投影中恢复出一幅自然图片的内容特征。假定已知原图为 K-sparse 类型,则可采用如下策略构建对应的 `A_prior_func` 和 `B_prior_func` 来反映这种特性约束条件: ```python from scipy.stats import norm K = ... # Sparsity parameter known beforehand lambda_ = float(K)/len(x_true) def soft_threshold(z, threshold): """Soft-threshold operator.""" abs_z = np.abs(z) signum = np.sign(z) rescaled_abs = np.maximum(abs_z - threshold, 0.) return signum * rescaled_abs def sparse_A_prior(rhats): global lambda_ thresholds = rhats.mean(axis=-1)*np.sqrt(lambda_) return soft_threshold(rhats, thresholds.reshape(-1,1)) def sparse_B_prior(_): global lambda_ return lambda_*np.ones_like(_.shape[:-1]) ``` 随后只需调用前述定义好的 `gamp_algorithm()` 方法即可完成整个流程执行工作流。 ---
评论
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符
 
红包 添加红包
表情包 插入表情
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值