M1を搭載した Macbook Proでbwaを実行してみる

# M1チップを搭載したMacbook proでbwa

ペアエンドのショートリードをbwaを用いてマッピングを行ってみた

手順は次の通り

- bwaのコンパイル
- ヒトゲノムのリファレンス配列のダウンロード
- bwaによるindexの作成
- Illumina のシーケンサーで取得されたリードデータをダウンロード
- マッピングの計算時間測定

コンパイルにはXcode Command Line Toolsで提供されるコンパイラーを利用 

```terminal
$ clang --version
Apple clang version 13.0.0 (clang-1300.0.29.3)
Target: arm64-apple-darwin21.1.0
Thread model: posix
InstalledDir: /Library/Developer/CommandLineTools/usr/bin
```

## bwaのコンパイル

`work`ディレクトリを作成し、bwa-0.7.17のソースコード(intel x86用)をダウンロード

```terminal
$ mkdir work
$ cd work
$ curl -LO https://github.com/lh3/bwa/archive/v0.7.17.tar.gz
```

ダウンロードしたソースコードを展開

```terminal
$ tar -zxf v0.7.17.tar.gz
```

bwaではIntel x86のSSE2を利用したコードが使われている
M1(Arm)のNEONを利用するため[SIMD Everywhere](https://github.com/simd-everywhere/simde)を使うが、[SIMD Everywhere](https://github.com/simd-everywhere/simde)のテストコードは不要なので、テストコード無しの[simde-no-tests](https://github.com/simd-everywhere/simde-no-tests)を使用する

```terminal
$ cd bwa-0.7.17
$ git clone https://github.com/simd-everywhere/simde-no-tests.git
```

sse2のコードをNEONのコードに対応させるため、`ksw.c`でインクルードしているヘッダーファイル`emmintrin.h`を`simde-no-tests/x86/sse2.h`に置き換える

```terminal
$ sed -i -e 's/<emmintrin.h>/"simde-no-tests\/x86\/sse2.h"/' ksw.c
```

適切なオプションを指定してコンパイル

```terminal
$ make clean
$ make -j 8 all CFLAGS="-Wall -Wno-unused-function -Ofast -mtune=native -fopenmp-simd　-DSIMDE_ENABLE_OPENMP -DSIMDE_ENABLE_NATIVE_ALIASES"
```

fileコマンドを利用してバイナリを確認

```terminal
$ file bwa
bwa: Mach-O 64-bit executable arm64
```

M1(arm)用のバイナリコードが作成できていることがわかる

##ヒトゲノムのリファレンス配列のダウンロード

一つ上の階層に`sample`ディレクトリを作成
`sample`ディレクトリ内に
ヒトゲノムのリファレンスファイル(GRCh37)をダウンロードして展開

```
$ mkdir ../sample
$ pushd ../sample
$ curl -LO ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/human_g1k_v37.fasta.gz
$ gunzip human_g1k_v37.fasta.gz
```

## bwaによるindexの作成

先ほどコンパイルした`bwa`を利用してリファレンスファイルのインデックスを作成する
複数のCPUを活用することはできないので、25分程必要

```
$ ../bwa-0.7.17/bwa index human_g1k_v37.fasta
```

次のファイルができていればindexが無事作成されたことになる

```
$ ls human_g1k_v37.fasta*
human_g1k_v37.fasta	human_g1k_v37.fasta.amb	human_g1k_v37.fasta.ann	human_g1k_v37.fasta.bwt	human_g1k_v37.fasta.pac	human_g1k_v37.fasta.sa
```

## Illumina のシーケンサーで取得されたリードデータをダウンロード

マッピングに必要な時間を計測するため、サンプルのリードをダウンロードする

- NIST7035_TAAGGCGA_L001_R1_001.fastq.gz
- NIST7035_TAAGGCGA_L001_R2_001.fastq.gz

ファイルサイズが大きいので注意
何かの問題でダウンロードの途中で失敗しても、下記のコマンドを繰り返せば、途中から再開してくれる
`-C -`のオプションは継続してダウンロードするためのもの

R1のリード

```
$ curl -LO -C - ftp://ftp-trace.ncbi.nih.gov/ReferenceSamples/giab/data/NA12878/Garvan_NA12878_HG001_HiSeq_Exome/NIST7035_TAAGGCGA_L001_R1_001.fastq.gz
```

R2のリード

```
$ curl -LO ftp://ftp-trace.ncbi.nih.gov/ReferenceSamples/giab/data/NA12878/Garvan_NA12878_HG001_HiSeq_Exome/NIST7035_TAAGGCGA_L001_R2_001.fastq.gz
```

ダウンロードしたリードのfastqファイルを展開

```
$ gunzip NIST7035_TAAGGCGA_L001_R1_001.fastq.gz
$ gunzip NIST7035_TAAGGCGA_L001_R2_001.fastq.gz
```

bwaをコンパイルしたディレクトリに戻る

```
$ popd
```

## マッピングの計算時間測定

準備ができたので、bwaを利用してマッピングを行う
測定条件は下記の通り

- Performanceコア(8個)のみ利用して 8スレッド
- Performanceコア(8個)とEfficient コア(2個)を利用して 10スレッド

10スレッドの際は、`-t 8`を`-t 10`に変更して計測

```
$ ./bwa mem -t 8 -o NIST7035_TAAGGCGA_L001.sam \
../sample/human_g1k_v37.fasta \
../sample/NIST7035_TAAGGCGA_L001_R1_001.fastq \ ../sample/NIST7035_TAAGGCGA_L001_R2_001.fastq \
```

時間の測定はそれぞれ2回行っている

### 8 スレッド: 平均 532秒

```
[main] Version: 0.7.17-r1188
[main] CMD: ./bwa mem -t 8 -o NIST7035_TAAGGCGA_L001.sam ../sample/human_g1k_v37.fasta ../sample/NIST7035_TAAGGCGA_L001_R1_001.fastq ../sample/NIST7035_TAAGGCGA_L001_R2_001.fastq
[main] Real time: 533.574 sec; CPU: 4281.706 sec
./bwa mem -t 8 -o NIST7035_TAAGGCGA_L001.sam ../sample/human_g1k_v37.fasta    4256.32s user 25.40s system 802% cpu 8:53.63 total
```

```
[main] Version: 0.7.17-r1188
[main] CMD: ./bwa mem -t 8 -o NIST7035_TAAGGCGA_L001.sam ../sample/human_g1k_v37.fasta ../sample/NIST7035_TAAGGCGA_L001_R1_001.fastq ../sample/NIST7035_TAAGGCGA_L001_R2_001.fastq
[main] Real time: 531.318 sec; CPU: 4263.938 sec
./bwa mem -t 8 -o NIST7035_TAAGGCGA_L001.sam ../sample/human_g1k_v37.fasta    4239.77s user 24.19s system 802% cpu 8:51.37 total
```

### 10 スレッド: 平均 477秒

```
[main] Version: 0.7.17-r1188
[main] CMD: ./bwa mem -t 10 -o NIST7035_TAAGGCGA_L001.sam ../sample/human_g1k_v37.fasta ../sample/NIST7035_TAAGGCGA_L001_R1_001.fastq ../sample/NIST7035_TAAGGCGA_L001_R2_001.fastq
[main] Real time: 478.127 sec; CPU: 4708.004 sec
./bwa mem -t 10 -o NIST7035_TAAGGCGA_L001.sam ../sample/human_g1k_v37.fasta    4680.20s user 27.84s system 984% cpu 7:58.18 total
```

```
[main] Version: 0.7.17-r1188
[main] CMD: ./bwa mem -t 10 -o NIST7035_TAAGGCGA_L001.sam ../sample/human_g1k_v37.fasta ../sample/NIST7035_TAAGGCGA_L001_R1_001.fastq ../sample/NIST7035_TAAGGCGA_L001_R2_001.fastq
[main] Real time: 476.498 sec; CPU: 4691.805 sec
./bwa mem -t 10 -o NIST7035_TAAGGCGA_L001.sam ../sample/human_g1k_v37.fasta    4663.68s user 28.17s system 984% cpu 7:56.56 total
```

### 他の環境との比較

参考までに同じマッピングを次のマシンで実行すると26分から27分の実行時間

- iMac Core(TM) i5-6600 CPU @ 3.30GHz
- MacBook Air  Core(TM) i7-1060NG7 CPU @ 1.20GHz

今回はこれまで:smile: