前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >算法(一)截取reads的算法

算法(一)截取reads的算法

作者头像
一只羊
发布2019-07-27 19:01:45
1.1K0
发布2019-07-27 19:01:45
举报
文章被收录于专栏:生信了

关键词:phred; trim; mott;

NGS(二代测序)分析的起点往往是fastq文件。fastq文件其实就是一条条的记录,每个记录包含4行。其中比较重要的是第二行和第四行:第二行是测序得到的碱基序列,第四行是每个碱基相应的测序质量,测序质量越高代表该碱基被测错的概率越低,反之越高。

正因为二代测序是有一定的错误率的,所以我们在进行下游分析之前,常常要对fastq文件中的reads进行修剪(trim),将一条reads中测序质量不高的部分截掉。那么截取reads常用的策略有两种,Fixed-length-trimming以及Phred-based-trimming。我们一一介绍:

Fixed-length-trimming

顾名思义,该方法就是截取固定长度的序列。一般来说,一条reads的头几个碱基和末尾几个碱基的测序质量比较差,所以你可以不加区分地将所有reads的前m个碱基以及后n个碱基去除。这种方法简单直接,但是不够精细。为什么这么说呢?因为每条reads测序质量差的区域长度并不固定,用一个固定的参数去截取reads两端往往会出出现“截取过度”或者“截取不足”的情况。

另外,有时候一条reads的非末端区域也会出现测序质量很差的碱基序列,那么这种从两头截取序列的策略就显得捉襟见肘了。综上,我们需要一种更为精细的截取方法。

Phred-based-trimming

该截取方法被phred软件所采用,又被称为modified Mott trimming algorithm,其描述如下:

简单来说,该方法可分为两步:

  1. 将原始测序质量phred值都减去一个阈值(默认0.05)得到一系列新的数值(该新的序列有正值也有负值)。
  2. 在该序列中找到和最大的子序列。

其中第2步是关键,稍有经验的读者应该发现第2步其实就是来自于“最大子序列和问题”,它表述如下:

给定整数A1, A2, ..., AN(可能有负数),求

的最大值(为方便起见,如果所有整数均为负数,则最大子序列和为0)。

例:

输入-2, 11, -4,13, -5, -2时,答案为20(从A2到A4)。

(来自《数据结构与算法分析:C语言描述》)

该问题有多种解法,包括暴力解法(brute force)-O(n2);递归解法-O(nlogn)以及线性时间解法-O(n) 等。其中效率最高的线性解法代码如下:

如果想知道目标子序列的下标,代码如下:

最后

小结一下,我们经常要对fastq文件中的reads进行修剪(trim),而常用的修剪策略有两种:Fixed-length-trimming以及Phred-based-trimming。其中,Fixed-length-trimming是截取固定长度的序列,该方法缺点是不够精细;而Phred-based-trimming通过寻找最大子序列和来进行截取,效果往往更精细。

值得注意的是,bwa以及seqtk中也采用了Phred-based-trimming算法。如果你有一个fastq文件想利用Phred-based-trimming算法进行修剪,可以安装seqtk后用一行命令实现:

seqtk trimfq your_fastq

当然,你也可以自行实现该算法,不过seqtk已经非常棒了,毕竟是lh3的作品嘛!

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2018-10-13,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 生信了 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档