79.《Bioinformatics Data Skills》之

某些基因组文件（包括GFF，BED与VCF类型）特别大，无法读入到内存。然而需要随机地读取某些区域的信息的话，不读入内存进行处理又特别费时。借助bgzip与tabix命令可以帮助我们解决这个问题。

这个过程包括以下三步：

1. 基于染色体与起始位置信息对文件内容进行排序
2. 使用bgzip工具对文件进行压缩
3. 使用tabix工具对文件建立索引

该方法的原理是利用bgzip可以按模块对文件进行压缩与解压的特点，可以在不解压整个文件的情况下读取特定模块的内容。

使用bgzip对文件进行压缩

以Mus_musculus.GRCm38.75.gtf.gz文件（下载）为例，排序之前需注意文件开头几行为注释：

$ zcat Mus_musculus.GRCm38.75.gtf.gz | head -n5
#!genome-build GRCm38.p2
#!genome-version GRCm38
#!genome-date 2012-01
#!genome-build-accession NCBI:GCA_000001635.4
#!genebuild-last-updated 2013-09

因此先将注释与内容拆开，对内容排序完再进行合并（详细原理见subshell），代码如下：

$ (zgrep "^#" Mus_musculus.GRCm38.75.gtf.gz; \
zgrep -v "^#" Mus_musculus.GRCm38.75.gtf.gz | sort -k1,1 -k4,4n) \
| bgzip > Mus_musculus.GRCm38.75.gtf.bgz

使用tabix建立索引

通过如下命令建立文件的索引。

$ tabix -p gff Mus_musculus.GRCm38.75.gtf.bgz

这里-p参数指定文件类型，Tabix也可以用于自定义的tsv文件，只要有染色体、起始与终止信息位置列即可。

之后当前目录下会生成了一个tbi文件：

$ ls *.tbi
Mus_musculus.GRCm38.75.gtf.bgz.tbi

查看文件内容

通过tabix <文件> <位置>来查看文件某个区域的内容，如下：

$ tabix Mus_musculus.GRCm38.75.gtf.bgz 16:23146536-23158028
16      protein_coding  exon    23146536 ...
16      protein_coding  gene    23146536 ...
16      protein_coding  transcript      23146536 ...
16      protein_coding  UTR     23146536 ...
...

结果可以重定向或者管道流处理，例如使用awk保留外显子区域：

$ tabix Mus_musculus.GRCm38.75.gtf.bgz 16:23146536-23158028 | awk '$3 ~ /exon/ {print}' 
16      protein_coding  exon    23146536 ...
16      protein_coding  exon    23146646
16      protein_coding  exon    23155217
16      protein_coding  exon    23155217
16      protein_coding  exon    23157074
16      protein_coding  exon    23157074

最后补充一点，tabix可以用于HTTP或者FTP服务器，将以上步骤处理的文件托管在服务器供别人远程使用。这里不再展开，感兴趣的话可以自行了解一下。