排序后的BAM具有更高的压缩率
Preface
前一段时间的数据分析产生了大量的文件，考虑到使用超算集群时有义务节约集群的储存空间，为其他用户的存储提供便利。因此这两天我开始检查自己的中间文件，打算删除无用的或者重复的中间文件以腾出更多的存储空间。
BAM文件排序前后对存储空间的占用不同
在我检查中间文件期间，我发现了自己曾经没有注意到的但是很有意思的现象。我分析的很多NGS测序数据都会使用samtools处理，我发现同一个样本的测序数在BAM文件排序前和排序后大小显著不同。例如我有一个BAM文件在排序前占用了208GB的储存空间，但是在排序后仅仅只占用92GB的储存空间。后者连前者的1/2都不到。这让我很是吃惊。我一度怀疑自己是不是在排序过程中同时进行了序列过滤的操作。因此我先检查了两个文件序列数量
samtools view sample.bam | wc -l; samtools view sample.srt.bam | wc -l
1310365258
1310365258
通过初步检查，发先序列数量没有发生变化，说明在samtools的排序处理时并没有发生序列被删除的情况。我又考虑到有没有可能在排序处理的过程中遗漏原始文件中序列的tag信息或者其他信息。如果这些信息被遗漏，确实有可能发生在序列总数不变的情况下BAM文件被缩小的结果。于是我查询了排序前和排序后BAM文件内部实际的一些序列的全部信息
samtools view sample.bam | head -n 1
A01340:58:HJHVNDSX7:1:1101:2284:1000 0 chr2 100174320 255 91M * 0 0 NTCAGAAATGCTTTCTCAAGGACATTCTTGTGTCCAGGCTTCAACAGAGCCAAGATACTTCATCCTCTTAGGGCGCCTGGTGAAACCTCGA #:FFFFFFFFFFFFFFFFFFF::FFFFFFFF:FFFFFFF:FFFFFF:FFFFF:FFFF::FFFFFFFFFFFFFFFFF:FFFFF:FFFFFF:F AS:i:-5 XN:i:0 XM:i:2 XO:i:0 XG:i:0 NM:i:2 MD:Z:0T0A89 YT:Z:UU XE:i:1
samtools view sample.srt.bam | grep A01340\:58\:HJHVNDSX7\:1\:1101\:2284\:1000
A01340:58:HJHVNDSX7:1:1101:2284:1000 0 chr2 100174320 255 91M * 0 0 NTCAGAAATGCTTTCTCAAGGACATTCTTGTGTCCAGGCTTCAACAGAGCCAAGATACTTCATCCTCTTAGGGCGCCTGGTGAAACCTCGA #:FFFFFFFFFFFFFFFFFFF::FFFFFFFF:FFFFFFF:FFFFFF:FFFFF:FFFF::FFFFFFFFFFFFFFFFF:FFFFF:FFFFFF:F AS:i:-5 XN:i:0 XM:i:2 XO:i:0 XG:i:0 NM:i:2 MD:Z:0T0A89 YT:Z:UU XE:i:1
我检查了多条序列发现都是一样的，因此我确定了BAM文件在排序前后序列文件内部的信息不会丢失。但是为什么排序后文件存储空间占用变小了，这是一个很有趣的问题。
BGZF压缩算法
我查阅了一下samtools的BAM文件对应的压缩算法，了解到该软件使用一种叫BGZF（Blocked GNU Zip Format）的压缩算法，这种算法是一种专门为BAM文件设计的压缩算法，结合了gzip和block的特性。BGZF将每个压缩块的大小限制在64KB，将整个文件分成多个这样的块。这种设计使得在BAM文件中可以轻松进行随机访问，而不需要解压整个文件。BGZF压缩算法的主要特点包括：1.支持快速的随机访问：由于将文件分成多个块，每个块的数据可以独立解压，从而实现对文件的快速随机访问。2.较高的压缩效率：与gzip相比，BGZF在保持相当的压缩比的情况下，具有更好的性能。3.足够高的解压速度：BGZF算法在进行解压操作时速度较快，可以有效减少读取时间。
但是受限于资料有限，我暂时没有找到关于BGZF压缩算法相关的数学模型。不过在搜索BGZF算法的过程中，我发现国外有学者也发现了samtools对bam文件排序后会减少文件的大小。并且samtools软件的开发者李恒也给这位国外学者做出了简单回复,大概意思是说由于在reads排序前，BAM文件中不同reads在染色体上的分布和坐标不同，排序前的bam文件中的reads是散乱分布。而排序后的reads做到了在基因组上的位置按照染色体和坐标的顺序排序，这样sort后的bam文件中相邻reads的一定在染色体的临近位置。这样使得BAM文件中每个相邻的压缩块更相似。我理解的是文本空间经过排序后的熵减少了，从而使得储存空间占用更小。
后记
虽然暂时没有搞清楚BGZF算法的原理，不过通过这次清理空间，我还是发现了如何更高效更节省地储存BAM文件。
samtools sort用法
Usage: samtools sort [options...] [in.bam]
Options:
-l INT Set compression level, from 0 (uncompressed) to 9 (best)
-u Output uncompressed data (equivalent to -l 0)
-m INT Set maximum memory per thread; suffix K/M/G recognized [768M]
-M Use minimiser for clustering unaligned/unplaced reads
-R Do not use reverse strand (only compatible with -M)
-K INT Kmer size to use for minimiser [20]
-I FILE Order minimisers by their position in FILE FASTA
-w INT Window size for minimiser indexing via -I ref.fa [100]
-H Squash homopolymers when computing minimiser
-n Sort by read name (natural): cannot be used with samtools index
-N Sort by read name (ASCII): cannot be used with samtools index
-t TAG Sort by value of TAG. Uses position as secondary index (or read name if -n is set)
-o FILE Write final output to FILE rather than standard output
-T PREFIX Write temporary files to PREFIX.nnnn.bam
--no-PG
Do not add a PG line
--template-coordinate
Sort by template-coordinate
--input-fmt-option OPT[=VAL]
Specify a single input file format option in the form
of OPTION or OPTION=VALUE
-O, --output-fmt FORMAT[,OPT[=VAL]]...
Specify output format (SAM, BAM, CRAM)
--output-fmt-option OPT[=VAL]
Specify a single output file format option in the form
of OPTION or OPTION=VALUE
--reference FILE
Reference sequence FASTA FILE [null]
-@, --threads INT
Number of additional threads to use [0]
--write-index
Automatically index the output files [off]
--verbosity INT
Set level of verbosity