算法（二）蓄水池抽样算法快速随机抽取reads

原创：hxj7

关键词：蓄水池算法;

fastq文件往往都很大，出于测试目的，我们经常要从fastq文件中随机抽取reads，生成一个小一点的fastq文件，以加快测试效率。假设我们要从一个包含大约100M reads的fastq文件中随机抽取1M reads，该怎么办呢？

我们将问题简单化：假设我们要从一个txt文件中（不知道总共多少行）随机抽取M行（fastq文件的处理与之类似，只不过fastq文件是压缩过的，且其一条记录由4行组成），比较容易想到的是如下办法（伪代码）：

image

该方法有一个很大的缺点，就是需要读取文件两遍。一般来说，I/O操作是很耗时的，所以该方法效率不高，当文件很大时，慢得惊人。

如果只要读取一次文件就好了，比如下面这样（伪代码）：

image

该方法把文件整个读入内存，的确减少了程序读取文件的总次数。但是，当文件很大时，该方法消耗的内存就太大了（想像一下把一个8G的txt文件整个加载到内存时的糟糕情况）。所以，不光要减少读取文件的次数，还要消耗较少的内存才好！

蓄水池抽样算法（Reservoir Sampling）就可以较好地解决上述问题，伪代码如下：

image

蓄水池抽样方法只需读取文件一次，且消耗的内存只有M行大小，而不是整个文件。所以，程序运行的效率会大大提高。

蓄水池抽样算法适用于大数据随机抽样，其关键在于证明其抽样的步骤是等概率的。其实证明方法也不难，只需运用归纳法即可，具体证明过程可参照wiki。值得注意的是，lh3大神的seqtk工具集就实现了该算法，只需运用seqtk sample命令。

今天的分享就到这里，谢谢大家！

(公众号：生信了)

image