Clumpify 是 BBMap 工具包中的一个组件,它与其他工具略有不同的是 Clumpify 不会实际更改你的数据,只是重新排序以最大化 gzip 压缩。因此,输出文件仍然是完全兼容的 gzipped fastq 文件,而 Clumpify 除了使其更快之外,对下游分析也不会造成影响。 它使用起来非常简单:
clumpify.sh in=reads.fq.gz out=clumped.fq.gz reorder
This command assumes paired, interleaved reads or single-ended reads. But Clumpify also now works with paired reads in twin files like this:
clumpify.sh in1=r1.fq.gz in2=r2.fq.gz out1=clumped1.fq.gz out2=clumped2.fq.gz reorder
Clumpify 工作原理
Clumpify 的操作原理与使排序的 bam 文件小于未排序 bam 文件的原理相似 —— reads 被重新排序,以便具有相似序列的 reads 都集中在附近,从而使得 gzip 压缩更有效。但是与 bam 排序不同,在此过程中,配对的 pairs 会保持在一起,这样交错的文件(interleaved file)将保持交错,配对完好无损。另外一个与 bam 排序不同的是,Clumpify 不需要比对(mapping)和参考序列(reference),除非在非常特殊的情况下,Clumpify 可以使用任意少量的内存。 因此,与 mapping 比对相比,它的速度非常快且内存效率高,并且可以在不知道读取来自哪些生物体的情况下完成。
在内部,Clumpify 形成了大量共享特殊 ‘pivot’ kmers 的 reads 簇(团块),也就是说这些 reads 是重叠的(overlap)。接着,通过 reads 中的 kmer 位置进一步对这些团块进行分类,使得在一个团块内的 reads 都是位置分类的(position-sorted)。 最终结果是得到一个排序的 reads 块列表,在排序的 bam 的百分比左右产生压缩(yielding compression within a percent or so of sorted bam)。
Clumpify 耗时
Clumpify 运行会非常快。如果所有数据都适合内存时,Clumpify 只需要一次读取和写入文件所需的时间。如果数据不能适合内存,则需要大约两倍的时间。 但是,我们的集群上有一个高性能的并行文件系统。 在某些文件系统或单个旋转磁盘上,当数据不适合内存时,可能需要几倍的时间,因为一次只能读取和写入多个文件。
为什么能提升速度
有很多进程是具有 I/O 限制的。 例如,在多核处理器上,使用 BBDuk,BBMerge,Reformat 等去解压缩一个 gzipped fastq 文件通常会受到 gzip 解压缩的速率限制(rate-limited)(即使你使用 pigz,它在解压缩方面比 gzip 快得多)。 Gzip 解压缩似乎会受到每秒输入字节数的速率限制而不是输出,这意味着一个给定大小的原始文件,如果压缩 Y%,将更快地解压缩 X%(meaning that a raw file of a given size will decompress X% faster if it is compressed Y% smaller);这里 X 和 Y 是成比例的,但不是 1 比 1。 在我的测试中,使用 Spades 和 Megahit 进行组装可以减少使用 Clumpified 输入所需的时间,而不仅仅是用于运行 Clumpify 所需的时间,这主要是因为两者都是 multi-kmer 组装软件(assemblers),它们多次读取输入文件(根据 Megahit 的作者, 由于缓存局部性)。 单纯的受 CPU 限制(CPU-limited)的程序和软件,如比对(mapping) 通常不会在速度方面受益太多(尽管仍然有点因为改进了缓存局部性)。
什么时候怎么使用 Clumpify
如果您想要对压缩数据进行进一步压缩,请尽早进行(例如,在 raw reads 时)。 然后运行所有下游处理步骤,确保维持 read 顺序(例如,如果使用 BBDuk 进行 adapter-trimming 时,则使用 “ordered” 标志),以便保持 clump 顺序; 因此,所有中间文件都将受益于增加的压缩和增加的速度。 我建议在所有将进入长期存储的数据上运行 Clumpify,或者在有多个步骤和中间 gzip 压缩文件的长管道时运行 Clumpify。 此外,即使数据不会进入长期存储,如果正在使用共享文件系统或需要通过互联网发送文件,尽早运行 Clumpify 将节省带宽。 下面列举了我不会聚集数据的唯一时间。
- 对于具有非常低 kmer 深度的 reads,由于覆盖率非常低(如 1x WGS)或超高错误率(如原始 PacBio 数据)。 它不会损坏任何东西,但也不会完成任何事情。
- 对于大量扩增子数据。 这可能有效,也可能无效; 但如果你的所有的 reads 都被期望共享相同的 kmers,它们可能都会形成一个巨大的丛,再也没有任何东西可以实现。 同样,它不会损坏任何东西,如果从可变区域中随机选择 pivots ,它可能会增加压缩。
- 当你的过程依赖 reads 的顺序时。如果你总是从文件中获取前一百万条 reads,假设它们在文件其余部分具有良好代表性,那么 Clumpify 将导致你的假设无效 —— 就像从分类的 bam 文件中获取前一百万条 reads 不具有代表性。 幸运的是,这绝不是一个好习惯,所以如果你现在正在这样做,现在无论如何都是改变你的流程的好机会。 随机子采样(Randomly subsampling)是一种更好的方法。
- If you are only going to read data fewer than ~3 times, it will never go into long-term storage, and it’s being used on local disk so bandwidth is not an issue, there’s no point in using Clumpify (or gzip, for that matter).
Clumpify使用与数据一致性检验
在实际测试用,使用 Hiseq2500 的 fastq 压缩数据 A,经 clumpify.sh 处理后得到压缩数据 B,但 A、B 解压出来的 fastq 数据 md5 是不一致的:
# 本测试所所用 bbmap 版本:38.20
$ cp B1_CNE2T1_R1.fastq.gz test_R1.fastq.gz
$ md5sum B1_CNE2T1_R1.fastq.gz test_R1.fastq.gz
0dc5da3a25d4f27e1ff02d2924497f5e B1_CNE2T1_R1.fastq.gz
0dc5da3a25d4f27e1ff02d2924497f5e test_R1.fastq.gz
$ sh clumpify.sh in=test_R1.fastq.gz out=clumpify_test_R1.fastq.gz
$ du -sh *
701M B1_CNE2T1_R1.fastq.gz
515M clumpify_test_R1.fastq.gz
701M test_R1.fastq.gz
$ gzip -d test_R1.fastq.gz
$ gzip -d clumpify_test_R1.fastq.gz
$ ll
total 8011828
-rw-r--r-- 1 shenweiyan bioinfo 734384226 Aug 9 16:39 B1_CNE2T1_R1.fastq.gz
-rw-r--r-- 1 shenweiyan bioinfo 3734854983 Aug 9 16:45 clumpify_test_R1.fastq
-rw-r--r-- 1 shenweiyan bioinfo 3734854983 Aug 9 16:40 test_R1.fastq
$ md5sum test_R1.fastq clumpify_test_R1.fastq
a3c411044914a671f88dd0976f29b237 test_R1.fastq
cf3b381d95725f1c2ddd352e50f60832 clumpify_test_R1.fastq
其中的一个原因是,clumpify 其实只是将 Fastq 文件根据序列相似性进行了位置重排(有点类似 CD-hit 工作原理),以便使得文件压缩率达到最大,它并没有对文件内容做任何改动。它除了会使后续分析流程变快之外,没有任何副作用。为了验证这一点,我们把原始的 Fastq 和经 clumpify 处理后解压出来的 Fastq 进行重新排序,并对排序后的结果进行一致性检测:
$ sort test_R1.fastq >sorted_test_R1.fastq
$ sort clumpify_test_R1.fastq >sorted_clumpify_test_R1.fastq
$ ll
total 15306476
-rw-r--r-- 1 shenweiyan bioinfo 734384226 Aug 9 16:39 B1_CNE2T1_R1.fastq.gz
-rw-r--r-- 1 shenweiyan bioinfo 3734854983 Aug 9 16:45 clumpify_test_R1.fastq
-rw-r--r-- 1 shenweiyan bioinfo 3734854983 Aug 9 17:43 sorted_clumpify_test_R1.fastq
-rw-r--r-- 1 shenweiyan bioinfo 3734854983 Aug 9 17:26 sorted_test_R1.fastq
-rw-r--r-- 1 shenweiyan bioinfo 3734854983 Aug 9 16:40 test_R1.fastq
$ md5sum *
0dc5da3a25d4f27e1ff02d2924497f5e B1_CNE2T1_R1.fastq.gz
cf3b381d95725f1c2ddd352e50f60832 clumpify_test_R1.fastq
240f0aa5ec8be86b38908f8bb09fb21d sorted_clumpify_test_R1.fastq
240f0aa5ec8be86b38908f8bb09fb21d sorted_test_R1.fastq
a3c411044914a671f88dd0976f29b237 test_R1.fastq
可以看出来原始的 Fastq 和经 clumpify 处理后解压出来的 Fastq 重新排序后得到的结果是一致的,说明 clumpify 的确没有对文件内容做任何改动。
因此,我们可以基于此建议对于数据需要长期存储,或者后续分析流程耗时较长的场景,数据下机之后首先用 clumpify 对 fastq 进行排序压缩,节省数据存储空间。
Clumpify使用注意事项
内存设置
实测 Clumpify 确实能显著减少 fastq 压缩文件的体积,这里注意一个参数 -Xmx:限制 JVM 使用的最大内存。一定要根据自己机器配置情况设置一下,不然 Clumpify 默认自动检测可用内存,如果处理的 Fastq 文件太大,会把机器内存都占满,就别再想跑别的程序了。
重复数据删除
Clumpify 的 37.24 版本起,增加了一些很好的光学重复数据删除的改进。 它现在更快(many times faster in certain degenerate cases),并且提高了 NextSeq tile-edge duplicates 的精度。 具体来说,现在建议使用下面的方式进行删除处理:
clumpify.sh in=nextseq.fq.gz out=clumped.fq.gz dedupe optical spany adjacent
这将删除所有正常的光学 duplicates,以及所有的 tile-edge duplicates,但是如果它们位于相邻的 tile 中并且共享它们的 Y 坐标(within dupedist),而不是之前,它会仅仅将 reads 视为 tile-edge duplicates。 它们可以在任何 tiles 中并且可以共享它们的 X 坐标。 这意味着假阳性较少(PCR 或巧合重复被分类为光学/tile-edge 重复)。 这可能是因为在 NextSeq 上,tile-edge duplicates 仅存在于 tile X-edges 上,并且 duplicates 仅在相邻 tiles 之间。
正常光学复制删除(HiSeq,MiSeq 等)应使用此命令:
clumpify.sh in=nextseq.fq.gz out=clumped.fq.gz dedupe optical
参考资料
- Brian Bushnell,《Tool: Introducing Clumpify: Create 30% Smaller, Faster Gzipped Fastq Files. And remove duplicates.》,Biostars
- wangpeng905,《Linux 下提升生物信息分析工作效率的神器(持续更新)》,简书