ゲノムの任意の位置の塩基配列の取得
ポイント
- 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.
今回はここまで