ショートリード配列から高速に薬剤耐性や病原性のタイピングをするSRST2を使う

## SRST2とは
whole genome sequencing (WGS)のショートリード配列を入力データとして、遺伝子/対立遺伝子、MLSTを高速かつ正確に検出するためのマッピングベースのツール。

[論文](https://genomemedicine.biomedcentral.com/articles/10.1186/s13073-014-0090-6)内では、一般的な病原体から得られた900以上のゲノムを用いて検証したところ、SRST2は高精度であり、遺伝子の検出と対立遺伝子のアサインの両方の点でアセンブリベースの手法よりも優れていると述べている。

またSRST2では薬剤耐性や病原性の迅速な検出も可能なため、ビクトリア州の微生物診断ユニット公衆衛生研究所（MDU PHL）内でのSRST2の検証を行い、病院での微生物ゲノムサーベイランスにSRST2を利用できることを実証した。

## インストール
SRST2は、Python2.7を指定しているので、condaで仮想環境を作成する。

```bash
# condaでpython2.7.5の仮想環境を作成
conda create -n srst2 python=2.7

git clone https://github.com/katholt/srst2

pip install srst2/
```

依存関係をインストールする。

- python2.7
- scipy
- bowtie2 v2.1.0
- SAMtools v0.1.18

SAMtoolsもbowtie2も、とても古いバージョンを指定している。
エラー出力により、bowtie2は、2.1.0から2.2.9まで対応していることがわかった。
> 12/28/2020 23:55:04 Incorrect version of bowtie installed.
12/28/2020 23:55:04 bowtie versions compatible with SRST2 are 2.1.0, 2.2.3, 2.2.4, 2.2.5, 2.2.6, 2.2.7, 2.2.8, 2.2.9

```bash
pip install scipy

# bowtie2 2.2系でcondaでインストールできる最新バージョンが2.2.5だった。
conda install bowtie2=2.2.5

# SAMtoolsはバージョンを指定する。
conda install samtools==0.1.18
```

SRST2の実行テスト。

```bash
srst2 --version
getmlst.py -h
slurm_srst2.py -h
```

## ケース1. MLST

### 1.入力ファイル
- fastqファイル

```bash
# Vietnam Shigella（赤痢菌）のPaired-endをENAからダウンロードする
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR024/ERR024070/ERR024070_1.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR024/ERR024070/ERR024070_2.fastq.gz
```

- fasta配列のデータベース

PubMLSTにて、[MLST allelic profiles and sequences](https://pubmlst.org/data/)から、E.coliのfastaファイルをダウンロードするが、#1と#2が存在するのでどちらかを指定する。

```bash
# 対立遺伝子配列のfastaファイルをダウンロードする
getmlst.py --species 'Escherichia coli#1'
```
コマンドを実行した際、以下のようなエラーが出た場合は、こちらを参考にしてエンコーディングを変更。
[Python スクリプト実行時に UnicodeDecodeError が出る場合の対処方法](https://shu223.hatenablog.com/entry/20111201/1328334689)
> UnicodeDecodeError: 'ascii' codec can't decode byte 0xef in position 16: ordinal not in range(128)

### 2. MLSTの実行

```bash
# MLSTの実行
srst2 \
--output test \
--input_pe ERR024070_* \
--mlst_db Escherichia_coli#1.fasta \
--mlst_definitions profiles_csv \
--mlst_delimiter '_'
```

### 3. 出力結果
以下のファイル`test__mlst__Escherichia_coli\#1__results.txt`が解析結果。

|Sample|ST|adk|fumC|gyrB|icd|mdh|purA|recA|mismatches|uncertainty|depth|maxMAF|
|---|---|---|---|---|---|---|---|---|---|---|---|---|
|ERR024070|NF*?|-|-|-|-|14*?|489*?|524*?|mdh_14/42holes;purA_489/36holes;recA_524/32holes|mdh_14/edge1.0;purA_489/edge1.0;recA_524/edge0.0|2.49966666667|0.5|

> so typing calls should be made based on the presence of genes rather than allele assignments (i.e. it is generally safe to ignore *? characters)

## ケース2. 薬剤耐性遺伝子の検出

### 1.入力ファイル
- fastqファイル（先ほどと同様のShigella）

- fasta配列の薬剤耐性遺伝子データベース
`data/ARGannot_r3.fasta`を用いる。薬剤耐性遺伝子のデータベースである[ResFinder database](https://cge.cbs.dtu.dk/services/ResFinder/)と[CARD](https://card.mcmaster.ca/)に基づいて作成されている。

### 2. 遺伝子検出の実行

```bash
# 遺伝子の検出
srst2 --input_pe ERR024070_* \
--gene_db ../data/ARGannot_r3.fasta \
--output test \
--log
```

### 3. 出力結果
以下のファイル`test__genes__ARGannot_r3__results.txt`が解析結果。

|Sample|AmpC1_Ecoli_Bla|StrB_AGly|
|---|---|---|
|ERR024070|AmpC1_1670*?|StrB_1614*?|

## 出力結果の説明

**「*」**は、最高スコアの対立遺伝子に1個以上のミスマッチがあることを示している。「mismatches」のカラムに詳細が出力される。この場合はたいてい、新規の対立遺伝子があること示している。
**「?」**は、最高スコアの対立遺伝子にdepthがあるため、出力結果に不確実性があることを示している。
**「-」**は、対立遺伝子がアサインされなかったということ。リードで90%以上がカバーされている対立遺伝子がなかったことを示している。

### MLSTの出力結果
**「ST」**カラムは、解析で、対立遺伝子の組み合わせが`--mlst_definitions`によって定義しているデータにヒットした場合、このSTはST列に記載される。**「NF」**は対立遺伝子の組み合わせが見つからなかったことを示し、**「ND」**はSTの計算が行われなかったことを示す（STの定義がされなかったため）。
**「mismatches」**カラムは、最高スコアの対立遺伝子に対するミスマッチが記載される。
**「uncertainty」**カラムは、coverageのdepthが低すぎた最高スコアの対立遺伝子の詳細を記載している。
**「depth」**カラムは、最高スコアの対立遺伝子の番号がアサインされた、すべての対立遺伝子の平均read depthを示している（つまり、'-'として記録されたものを除く）。
**「maxMAF」**カラムは、MLSTのアライメントで見つかったバリアントの中で最も高いMinor allele frequency (MAF)を記載する。この値は最小値0から最大値0.5。

## 複数ジョブを同時に実行
ジョブ管理システムである「Slurm」に`slurm_srst2.py`を用いることで、複数のPaired-endリードセットで解析を実行することもできる。
`--other_args`で、まとめて引数をSRST2に渡す。

```bash
slurm_srst2.py --script srst2 \
--output test \
--input_pe *.fastq.gz \
--other_args '--gene_db resistance.fasta --mlst_db Escherichia_coli.fasta --mlst_definitions ecoli.txt --save_scores'
--walltime 0-1:0 \
> job_sub_list.txt
```

## 参考リンク
[SRST2: Rapid genomic surveillance for public health and hospital microbiology labs](https://genomemedicine.biomedcentral.com/articles/10.1186/s13073-014-0090-6)
November 20, 2014

[SRST2 github](https://github.com/katholt/srst2/blob/master/README.md)

[全ゲノムシーケンスの決定から検出アッセイの確立まで](https://www.thermofisher.com/content/dam/LifeTech/Documents/PDFs/jp/materials/appnote-ion012.pdf)

> MLST 法は、数種類の細菌のハウスキーピング遺伝子由来の配列をシーケンシングし菌株の同定を行います。このフラグメントの長さは約450～500塩基対で、ここでシーケンシングした配列はそれぞれの遺伝子を確実にプロファイルしています。これらの遺伝子座で同定される対立遺伝子から、標的分離株の対立遺伝子プロファイルすなわち遺伝子型（sequence type: ST）を定義し、以前に同定された株のデータベースと比較して、推定同定を行うことができます [4]。

[耐性カセットや病原性遺伝子のタイピングを行う SRST2](http://kazumaxneo.hatenablog.com/entry/2019/02/21/112350)
Feburary 22, 2019

> bowite2の対応バージョンは 2.1.0, 2.2.3, 2.2.4, 2.2.5, 2.2.6, 2.2.7, 2.2.8, 2.2.9のみ。
> samtoolsの対応バージョンは 0.1.18, 1.0, 1.1, 1.2, 1.3, (0.1.18 is recommended)