ゲノムの任意の位置の塩基配列の取得

ポイント

bgzipでFastaファイルを圧縮しておく
- fastaファイルの取得時に圧縮がかかっている場合は、まず展開する
samtools faidx コマンドを利用してFasta ファイルのインデックスファイルを作成する
- 拡張子が.faiと拡張子が.gziのファイルが作成される
samtools faidx コマンドを利用して領域を取得する

Ensemblからヒトゲノムのリファレンスを取得

まずはEnsemblからヒトゲノムのリファレンスを取得しよう

GRCh37の場合

$ curl -L -O ftp://ftp.ensembl.org/pub/grch37/current/fasta/homo_sapiens/dna/Homo_sapiens.GRCh37.dna.primary_assembly.fa.gz

gzip -d を利用して圧縮されたファイルを展開
bgzip を利用して展開されたファイルを再度圧縮（理由は後述）

$ gzip -d Homo_sapiens.GRCh37.dna.primary_assembly.fa.gz
$ bgzip Homo_sapiens.GRCh37.dna.primary_assembly.fa

samtools faidxコマンドを利用して圧縮ファイルのインデックスファイルを作成しよう

$ samtools faidx Homo_sapiens.GRCh37.dna.primary_assembly.fa.gz

拡張子が.faiと拡張子が.gziのファイルが新たに作成された事を確認

$ ls     
Homo_sapiens.GRCh37.dna.primary_assembly.fa.gz  Homo_sapiens.GRCh37.dna.primary_assembly.fa.gz.fai  Homo_sapiens.GRCh37.dna.primary_assembly.fa.gz.gzi

任意の位置の塩基配列の抽出

準備ができたので GRCh37のゲノム上の塩基配列を抽出してみよう

利用するコマンドはsamtools faidxでインデックスを作成するコマンドと同じ

samtools faidx 対象とするfasta配列クロモゾムのID:開始点(1-base)-終了点(1-base)

ヒト変異のHGVS表記の変換(vep api)で取得したゲノム上の変異を表すhgvs表記"NC_000007.13:g.150696111T>G"を利用する

上記のhgvs表記はGRCh37の7番 Chromosomeの150696111の塩基TがGに変化している事を表している
そこで、リファレンス上の該当塩基がTかどうかを確認してみよう

$ samtools faidx Homo_sapiens.GRCh37.dna.primary_assembly.fa.gz 7:150696111-150696111
>7:150696111-150696111
T

変異位置の前20ベース及び変異の位置の後20ベースを取得

$ samtools faidx Homo_sapiens.GRCh37.dna.primary_assembly.fa.gz 7:150696091-150696110 
>7:150696091-150696110
CTGCTGCAGGCCCCAGATGA
$ samtools faidx Homo_sapiens.GRCh37.dna.primary_assembly.fa.gz 7:150696112-150696131
>7:150696112-150696131
CCCCCAGAACTCTTCCTTCT

GRCh38の場合

GRCh38のリファレンスゲノムをEnsemblから取得してインデックスファイルを作成する
取得するためのアドレスが異なる以外はGRCh37のそれと同様

$ curl -L -O ftp://ftp.ensembl.org/pub/release-95/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
$ gzip -d Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz 
$ bgzip Homo_sapiens.GRCh38.dna.primary_assembly.fa 
$ samtools faidx Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz

ヒト変異のHGVS表記の変換(vep api)で取得したGRCh38のゲノム上の変異を表すhgvs表記"NC_000007.14:g.150999023T>G"を利用する

$ samtools faidx Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz 7:150999023-150999023
>7:150999023-150999023
T
$ samtools faidx Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz 7:150999003-150999022
>7:150999003-150999022
CTGCTGCAGGCCCCAGATGA
$ samtools faidx Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz 7:150999024-150999043
>7:150999024-150999043
CCCCCAGAACTCTTCCTTCT

bgzipについて

bgzipで圧縮したファイルにはこんな特徴がある

ブロック毎に圧縮されているので、必要な部分を含むブロックを展開できる
- そのために拡張子が.gziのファイルを利用する
bzgip -d, gzip -d または gunzipコマンドを利用して展開できる

bgzipは、ディスクアクセスが十分高速な場合には複数のスレッドを利用する事で高速化できる

--threads スレッド数のオプションを利用する事で圧縮時間を削減できる
- ただしディスクの読み書きが十分に高速で無ければスレッド数を増やしても実質の処理時間は変わらない
- ホストのCPUの個数(スレッド数)以上に設定しても効率は良くならない

良くある問題

圧縮されているfastaファイルにインデックスを作成する場合には、必ずbgzipで圧縮したものでなければだめ
例えば、Ensemblからダウンロードした圧縮されているfastaファイルをそのまま利用して、インデックスを作成しようとすると、次のようなエラーとなる

$ samtools faidx Homo_sapiens.GRCh37.dna.primary_assembly.fa.gz
[E::fai_build3] Cannot index files compressed with gzip, please use bgzip
Could not build fai index Homo_sapiens.GRCh37.dna.primary_assembly.fa.gz.fai

各ツールをubuntu(Debian系列のLinux)に導入する方法

実行環境は以下の通り

$ grep 'VERSION' /etc/*release 
/etc/os-release:VERSION="18.04.1 LTS (Bionic Beaver)"
/etc/os-release:VERSION_ID="18.04"
/etc/os-release:VERSION_CODENAME=bionic

基本的にはapt-getでインストールすれば良い

bgzip と samtoolsのインストール

# apt-get update
# apt-get install -y tabix
# apt-get install -y samtools

どうしても各ツールの最新版を導入したい場合

samtoolsの最新版をソースコードからインストールする方法(少なくともfaidxコマンドを利用するだけなら、最新版を無理にインストールする必要は無い)

githubのreleasesを確認して、最新版のAssetsを確認する
https://github.com/samtools/samtools
https://github.com/samtools/htslib

# apt-get install -y build-essential
# apt-get install -y zliblg-dev libbz2-dev liblzma-dev libncurses5-dev
# curl -L -O https://github.com/samtools/samtools/releases/download/1.9/samtools-1.9.tar.bz2
# tar -jxf samtools-1.9.tar.bz2
# pushd samtools-1.9
# ./configure
# make
# make test
# make install
# popd
# curl -L -O https://github.com/samtools/htslib/releases/download/1.9/htslib-1.9.tar.bz2
# tar -jxf htslib-1.9.tar.bz2 
# pushd htslib-1.9.tar.bz2
# ./configure
# make
# make test
# make install
# popd

make testを実行してツールが正常に作成されたことを確認してからインストールを実行している

上記で samtools, bgzip(tabixも含む) が/usr/local/binにインストールされる

$ /usr/local/bin/bgzip --version
bgzip (htslib) 1.9
Copyright (C) 2018 Genome Research Ltd.
$ /usr/local/bin/samtools --version
samtools 1.9
Using htslib 1.9
Copyright (C) 2018 Genome Research Ltd.

今回はここまで