60.《Bioinformatics Data Skills》之

编写代码提取参考基因组上特定区域序列并不是一件困难的事情，但是存在一个问题：速度。有时我们关注的是很小的一个片段，比如说8:123,407,082-123,407,182区域，如果将整个基因组都导入内存再查找的话将非常低效（磁盘的读写速度很慢），有可能我们导入了上百兆的数据只是为了那几KB的序列。一个更为巧妙的解决方案是为FASTA文件建立索引，它会帮助我们定位指定区域，提取特定区域时会变得更加快速。

为了方便展示，这里以小鼠的8号染色体参考基因组文件（Mus_musculus.GRCm38.75.dna.chr omosome.8.fa.gz）为例子：

本节数据下载

首先我们对数据进行解压(samtools faidx可以直接处理压缩文件，但是只对BGZF压缩文件有效)：

gunzip Mus_musculus.GRCm38.75.dna.chromosome.8.fa.gz

接着使用samtools工具建立FASTA文件的索引：

samtools faidx Mus_musculus.GRCm38.75.dna.chromosome.8.fa

当前目录下将会出现Mus_musculus.GRCm38.75.dna.chr omosome.8.fa.fai文件

之后我们可以使用如下命令提取指定的区域：

$ samtools faidx Mus_musculus.GRCm38.75.dna.chromosome.8.fa 8:123,407,082-123,407,182
>8:123,407,082-123,407,182
GAGAAAAGCTCCCTTCTTCTCCAGAGTCCCGTCTACCCTGGCTTGGCGAGGGAAAGGAAC
CAGACATATATCAGAGGCAAGTAACCAAGAAGTCTGGAGGT

注：指定区域时注意你的fasta文件的染色体编号方式是纯数字还是"chr"开头，不正确的命名不会返回结果。

也可以同时查看多个区域，以空格隔开

$ samtools faidx Mus_musculus.GRCm38.75.dna.chromosome.8.fa 8:123,407,082-123,407,182 8:123,407,282-123,407,382 8:123,407,482-123,407,582
>8:123,407,082-123,407,182
GAGAAAAGCTCCCTTCTTCTCCAGAGTCCCGTCTACCCTGGCTTGGCGAGGGAAAGGAAC
CAGACATATATCAGAGGCAAGTAACCAAGAAGTCTGGAGGT
>8:123,407,282-123,407,382
AGAGCTCAGGGTCCTCAGCAGGAAGTGTCTATGCCATGCCGAGGCTGGCCTGTCCAGCCA
GAAAGAACACAAGTGTAAAGGAAAATCGGAGCGTGCCTGTA
>8:123,407,482-123,407,582
AGACCGCTTCCTACTTCCTGACAAGACTATGTCCACTCAGGAGCCCCAGAAGAGTCTTCT
GGGTTCTCTCAACTCCAATGCCACCTCTCACCTTGGACTGG