60.《Bioinformatics Data Skills》之建立FASTA文件索引

编写代码提取参考基因组上特定区域序列并不是一件困难的事情,但是存在一个问题:速度。有时我们关注的是很小的一个片段,比如说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
©著作权归作者所有,转载或内容合作请联系作者
【社区内容提示】社区部分内容疑似由AI辅助生成,浏览时请结合常识与多方信息审慎甄别。
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

相关阅读更多精彩内容

友情链接更多精彩内容