ゲノムの任意の位置の塩基配列を取得したい

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

##ポイント

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

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

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

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

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

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

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

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

```teminal
$ 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)](https://qiita.com/percipere/items/7dcc714a57161c5dc0b4)で取得したゲノム上の変異を表すhgvs表記"NC_000007.13:g.150696111T>G"を利用する

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

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

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

```terminal
$ 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のそれと同様

```terminal
$ 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)](https://qiita.com/percipere/items/7dcc714a57161c5dc0b4)で取得したGRCh38のゲノム上の変異を表すhgvs表記"NC_000007.14:g.150999023T>G"を利用する

```terminal
$ 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)に導入する方法

実行環境は以下の通り

```terminal
$ 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のインストール

```terminal
# 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

```terminal
# 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`にインストールされる

```terminal
$ /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.
```

今回はここまで:smile: