原创：hxj7
关键词：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两端往往会出出现“截取过度”或者“截取不足”的情况。
image
另外，有时候一条reads的非末端区域也会出现测序质量很差的碱基序列，那么这种从两头截取序列的策略就显得捉襟见肘了。综上，我们需要一种更为精细的截取方法。
Phred-based-trimming
该截取方法被phred软件所采用，又被称为modified Mott trimming algorithm，其描述如下：
image
简单来说，该方法可分为两步：
- 将原始测序质量phred值都减去一个阈值（默认0.05）得到一系列新的数值（该新的序列有正值也有负值）。
- 在该序列中找到和最大的子序列。
其中第2步是关键，稍有经验的读者应该发现第2步其实就是来自于“最大子序列和问题”，它表述如下：
（来自《数据结构与算法分析：C语言描述》）
该问题有多种解法，包括暴力解法（brute force）-O(n2)；递归解法-O(nlogn）以及线性时间解法-O(n) 等。其中效率最高的线性解法代码如下：
image
如果想知道目标子序列的下标，代码如下：
image
最后
小结一下，我们经常要对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后用一行命令实现：
<font color="green">seqtk trimfq your_fastq</font>
当然，你也可以自行实现该算法，不过seqtk已经非常棒了，毕竟是lh3的作品嘛！
(公众号：生信了)
image
附注：Phred-based-trimming的更多细节（来源）：
A somewhat simplified way of describing the algorithm is as follows: PHRED finds the longest fragment in a sequence where the estimated error rate is below the cutoff value. The trimmed sequence cannot be extended at either end without adding segments that have an error rate above the cutoff values. However, small segments within the trimmed sequence may have higher error rates &endash; sequencing artifacts like stops and compressions are an example.
The average estimated error rate of the trimmed sequence will generally be lower than the cutoff value. Statistically, the actual error rates of trimmed sequences will also be lower than the cutoff. However, any single sequence (or small number of sequences) may have an actual error rate that is higher, due to statistical sampling errors. (For example, if 10 bases have a 10% error probability, than 9 of the 10 will have an actual error rate of 0%, while one of the 10 will have a 100% actual error rate).
Finally, the trimming requires that the trimmed sequence is at least 20 bases long. Any sequence that is shorter after trimming will be completely trimmed, with no bases left over. The reason for this is that even traces that are clearly artifacts and do not contain any real sequence will often have short stretches (5-10 bases) where the base call quality is high. By enforcing the minimum length requirement, PHRED makes sure that garbage is always identified as such when trimming is used.
网友评论