Clumpify:能使 Fastq 压缩文件再缩小 30% 并加速后续分析流程
由于微信不允许外部链接,你需要点击文章尾部左下角的 "阅读原文",才能访问文中链接。
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 的操作原理与使排序的 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 只需要一次读取和写入文件所需的时间。如果数据不能适合内存,则需要大约两倍的时间。 但是,我们的集群上有一个高性能的并行文件系统。 在某些文件系统或单个旋转磁盘上,当数据不适合内存时,可能需要几倍的时间,因为一次只能读取和写入多个文件。
有很多进程是具有 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)通常不会在速度方面受益太多(尽管仍然有点因为改进了缓存局部性)。
如果您想要对压缩数据进行进一步压缩,请尽早进行(例如,在 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).
在实际测试用,使用 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 确实能显著减少 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 下提升生物信息分析工作效率的神器(持续更新)》, 简书
JGI,《BBTools User Guide: Clumpify Guide》
·end·
—如果喜欢,快分享给你的朋友们吧—
我们一起愉快的玩耍吧
本文分享自微信公众号 - 生信科技爱好者(bioitee)。
如有侵权,请联系 support@oschina.cn 删除。
本文参与“OSC源创计划”,欢迎正在阅读的你也加入,一起分享。