macでインフォマティクス

macでインフォマティクス

NGS関連のインフォマティクス情報についてまとめています。

耐性カセットや病原性遺伝子のタイピングを行う SRST2

 

病原性細菌の迅速な分子タイピングは、公衆衛生疫学、サーベイ ランスおよび感染制御にとって重要である[論文より ref.1,2]。そのような活動の2つの重要な目標は、(1)病原性遺伝子、抗菌性耐性遺伝子または血清型決定因子を含む臨床的に関連する表現型に関連する遺伝子の存在を検出すること、 (2)MLST [ref.3]またはクローン特異的または他の疫学マーカーの検出を介して、分離株をクローン群に分類すること。全ゲノムシーケンシング(WGS)またはゲノム疫学は、これらの作業にますます採用されており、PCR および/または制限酵素消化と、電気泳動 による配列決定またはサイズ分離とを組み合わせた現在の技術を置き換える可能性がある[ ref.4]。 WGSは、(1)生物または標的特異的試薬を必要としない多数の細菌分離株に同時に適用することができる。 (2)結果データは容易に共有可能であり、過去および将来のデータセットと容易に比較することができ、日常的なサーベイランス(遺伝子およびクローンのモニタリング)および詳細なアウトブレイク調査に利用可能である[ref.4]。

 WGSは病原体の研究に革命をもたらし、公衆衛生疫学、サーベイランス、感染管理の実践に革命をもたらす可能性がある[ref.4-10]。いくつかの実証研究[ref.11-16]にもかかわらず、WGSの日常的な使用は、公衆衛生および診断研究所にとって重要な課題を提起している。その主なものは、迅速で再現性のある情報抽出、解釈可能、rawシーケンシングデータの扱いである。

 現在利用可能な方法は、BLASTまたは他の検索アルゴリズムを用いて目的の遺伝子または対立遺伝子(例えば、ARG-Annot [ref.18]; ResFinder、PlasmidFinderおよびMLSTを同定することができるより長い連続した配列(コンティグ)へショートリードをアセンブリすることに依存している[ref.19-21]; BIGSdb [ref.22-23])。アセンブリへの依存は、ショートリードから高品質の細菌ゲノムアセンブリを生成するためのデータ、時間、および計算上の要求に起因する効率および感度の問題を引き起こす。数ギガバイトのメモリで数時間から数時間で細菌ゲノムアセンブリを生成できるアセンブラがいくつかある(例えば、Velvet [ref.24]、SPAdes [ref.25])。しかし、これらのツールを用いた高品質のアセンブリの作成には、質の高いフィルタリングやリードの前処理、kmerの長さや他のパラメータの最適化が必要であり、実際にはいくつかの代替アセンブリを生成して比較する必要がある。さらに、高度に最適化されたアセンブリの品質は、マルチプレックスで一緒にシーケンシングされたclosely relatedなゲノムであっても、非常に可変である。したがって、ショートリード技術でアセンブリベースで決定されたゲノムの解析では、公衆衛生および感染対策に使用するための堅牢で信頼性の高い再現性のあるアッセイを確実にするために重要な標準化と品質管理が非常に困難である。

 この論文では、ゲノム疫学の新しいツールであるSRST2について説明する。SRST2は、WGSのショートシークエンシングリードから直接的に遺伝子や対立遺伝子を迅速かつ正確に検出する。 SRST2は、任意の配列データベースを使用して読み込みを行い、MLST形式のデータベース[ref.3]で定義されたコンビナトリアルシーケンスタイプを計算することができる。(一部略)

 SRST2は、以前のツールSRST(Short Read Sequencing Typing)にちなんで名付けられたが、SRST2は完全に斬新なコードで、SRSTとは異なるリードマッピング、スコアリングとアルゴリズムを使用し、 MLSTと同様に、遺伝子の検出やアレルのタイピングを行う。

 

Labo HP

Holt Lab | microbial genomics

SRST2 HP

https://katholt.github.io/srst2/

 

インストール

ubuntu18.04のpython2.7.14でテストした。

依存

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)

本体 Github

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

#動作確認
srst2 --version
getmlst.py -h
slurm_srst2.py -h

> srst2 --version

$ srst2 --version

srst2 0.2.0

kazu@74a7df8f35c0:/data$ srst2 -h       

usage: srst2 [-h] [--version] [--input_se INPUT_SE [INPUT_SE ...]]

             [--input_pe INPUT_PE [INPUT_PE ...]] [--merge_paired]

             [--forward FORWARD] [--reverse REVERSE] [--read_type {q,qseq,f}]

             [--mlst_db MLST_DB] [--mlst_delimiter MLST_DELIMITER]

             [--mlst_definitions MLST_DEFINITIONS]

             [--mlst_max_mismatch MLST_MAX_MISMATCH]

             [--gene_db GENE_DB [GENE_DB ...]] [--no_gene_details]

             [--gene_max_mismatch GENE_MAX_MISMATCH]

             [--min_coverage MIN_COVERAGE] [--max_divergence MAX_DIVERGENCE]

             [--min_depth MIN_DEPTH] [--min_edge_depth MIN_EDGE_DEPTH]

             [--prob_err PROB_ERR]

             [--truncation_score_tolerance TRUNCATION_SCORE_TOLERANCE]

             [--stop_after STOP_AFTER] [--other OTHER]

             [--max_unaligned_overlap MAX_UNALIGNED_OVERLAP] [--mapq MAPQ]

             [--baseq BASEQ] [--samtools_args SAMTOOLS_ARGS] --output OUTPUT

             [--log] [--save_scores] [--report_new_consensus]

             [--report_all_consensus] [--use_existing_bowtie2_sam]

             [--use_existing_pileup] [--use_existing_scores]

             [--keep_interim_alignment] [--threads THREADS]

             [--prev_output PREV_OUTPUT [PREV_OUTPUT ...]]

 

SRST2 - Short Read Sequence Typer (v2)

 

optional arguments:

  -h, --help            show this help message and exit

  --version             show program's version number and exit

  --input_se INPUT_SE [INPUT_SE ...]

                        Single end read file(s) for analysing (may be gzipped)

  --input_pe INPUT_PE [INPUT_PE ...]

                        Paired end read files for analysing (may be gzipped)

  --merge_paired        Switch on if all the input read sets belong to a

                        single sample, and you want to merge their data to get

                        a single result

  --forward FORWARD     Designator for forward reads (only used if NOT in

                        MiSeq format sample_S1_L001_R1_001.fastq.gz; otherwise

                        default is _1, i.e. expect forward reads as

                        sample_1.fastq.gz)

  --reverse REVERSE     Designator for reverse reads (only used if NOT in

                        MiSeq format sample_S1_L001_R2_001.fastq.gz; otherwise

                        default is _2, i.e. expect forward reads as

                        sample_2.fastq.gz

  --read_type {q,qseq,f}

                        Read file type (for bowtie2; default is q=fastq; other

                        options: qseq=solexa, f=fasta).

  --mlst_db MLST_DB     Fasta file of MLST alleles (optional)

  --mlst_delimiter MLST_DELIMITER

                        Character(s) separating gene name from allele number

                        in MLST database (default "-", as in arcc-1)

  --mlst_definitions MLST_DEFINITIONS

                        ST definitions for MLST scheme (required if mlst_db

                        supplied and you want to calculate STs)

  --mlst_max_mismatch MLST_MAX_MISMATCH

                        Maximum number of mismatches per read for MLST allele

                        calling (default 10)

  --gene_db GENE_DB [GENE_DB ...]

                        Fasta file/s for gene databases (optional)

  --no_gene_details     Switch OFF verbose reporting of gene typing

  --gene_max_mismatch GENE_MAX_MISMATCH

                        Maximum number of mismatches per read for gene

                        detection and allele calling (default 10)

  --min_coverage MIN_COVERAGE

                        Minimum %coverage cutoff for gene reporting (default

                        90)

  --max_divergence MAX_DIVERGENCE

                        Maximum %divergence cutoff for gene reporting (default

                        10)

  --min_depth MIN_DEPTH

                        Minimum mean depth to flag as dubious allele call

                        (default 5)

  --min_edge_depth MIN_EDGE_DEPTH

                        Minimum edge depth to flag as dubious allele call

                        (default 2)

  --prob_err PROB_ERR   Probability of sequencing error (default 0.01)

  --truncation_score_tolerance TRUNCATION_SCORE_TOLERANCE

                        % increase in score allowed to choose non-truncated

                        allele

  --stop_after STOP_AFTER

                        Stop mapping after this number of reads have been

                        mapped (otherwise map all)

  --other OTHER         Other arguments to pass to bowtie2 (must be escaped,

                        e.g. "\--no-mixed".

  --max_unaligned_overlap MAX_UNALIGNED_OVERLAP

                        Read discarded from alignment if either of its ends

                        has unaligned overlap with the reference that is

                        longer than this value (default 10)

  --mapq MAPQ           Samtools -q parameter (default 1)

  --baseq BASEQ         Samtools -Q parameter (default 20)

  --samtools_args SAMTOOLS_ARGS

                        Other arguments to pass to samtools mpileup (must be

                        escaped, e.g. "\-A").

  --output OUTPUT       Prefix for srst2 output files

  --log                 Switch ON logging to file (otherwise log to stdout)

  --save_scores         Switch ON verbose reporting of all scores

  --report_new_consensus

                        If a matching alleles is not found, report the

                        consensus allele. Note, only SNP differences are

                        considered, not indels.

  --report_all_consensus

                        Report the consensus allele for the most likely

                        allele. Note, only SNP differences are considered, not

                        indels.

  --use_existing_bowtie2_sam

                        Use existing SAM file generated by Bowtie2 if

                        available, otherwise they will be generated

  --use_existing_pileup

                        Use existing pileups if available, otherwise they will

                        be generated

  --use_existing_scores

                        Use existing scores files if available, otherwise they

                        will be generated

  --keep_interim_alignment

                        Keep interim files (sam & unsorted bam), otherwise

                        they will be deleted after sorted bam is created

  --threads THREADS     Use multiple threads in Bowtie and Samtools

  --prev_output PREV_OUTPUT [PREV_OUTPUT ...]

                        SRST2 results files to compile (any new results from

                        this run will also be incorporated)

> getmlst.py -h

$ getmlst.py -h

usage: getmlst.py [-h] [--repository_url URL] --species NAME

                  [--force_scheme_name]

 

Download MLST datasets by speciesfrom pubmlst.org.

 

optional arguments:

  -h, --help            show this help message and exit

  --repository_url URL  URL for MLST repository XML index

  --species NAME        The name of the species that you want to download

                        (e.g. "Escherichia coli")

  --force_scheme_name   Flage to force downloading of specific scheme name

                        (e.g. "Clostridium difficile")

> slurm_srst2.py -h

$ slurm_srst2.py -h

usage: slurm_srst2.py [-h] [--walltime WALLTIME] [--memory MEMORY]

                      [--rundir RUNDIR] [--threads THREADS] --script SCRIPT

                      --output OUTPUT [--input_se INPUT_SE [INPUT_SE ...]]

                      [--input_pe INPUT_PE [INPUT_PE ...]] [--forward FORWARD]

                      [--reverse REVERSE] --other_args OTHER_ARGS

 

Submit SRST2 jobs through SLURM

 

optional arguments:

  -h, --help            show this help message and exit

  --walltime WALLTIME   wall time (default 0-1:0 = 1 h)

  --memory MEMORY       mem (default 4096 = 4gb)

  --rundir RUNDIR       directory to run in (default current dir)

  --threads THREADS     number of CPUs per job

  --script SCRIPT       path to srst2.py

  --output OUTPUT       identifier for outputs (will be combined with read set

                        identifiers)

  --input_se INPUT_SE [INPUT_SE ...]

                        Single end read file(s) for analysing (may be gzipped)

  --input_pe INPUT_PE [INPUT_PE ...]

                        Paired end read files for analysing (may be gzipped)

  --forward FORWARD     Designator for forward reads (only used if NOT in

                        MiSeq format sample_S1_L001_R1_001.fastq.gz; otherwise

                        default is _1, i.e. expect forward reads as

                        sample_1.fastq.gz)

  --reverse REVERSE     Designator for reverse reads (only used if NOT in

                        MiSeq format sample_S1_L001_R2_001.fastq.gz; otherwise

                        default is _2, i.e. expect forward reads as

                        sample_2.fastq.gz)

  --other_args OTHER_ARGS

                        string containing all other arguments to pass to srst2

 

またはdockerイメージを使う。対応するsamtoolsやbowtie2のバージョンが古いので、仮想環境で走らす方が無難です。

docker pull estrain/srst2

#ホストのカレントと/dataをシェアしてラン
docker run -itv $PWD:/data/ estrain/srst2
> srst2 -h

  

データベースの準備

fastaファイルのデータベースを準備する。準備されているスクリプトを使えば、PubMLST(リンク)から自動でダウンロードできる。検出対象の種を指定する。

getmlst.py --species "Escherichia coli#1" 
#e.coliは#1と#2の二つあると警告が出たのでここでは#1とした。詳細はPubMLSTで確認

複数のファイルがダウンロードされる。

f:id:kazumaxneo:20171218115905j:plain

Escherichia_coli.fasta 

user$ head -15 Escherichia_coli.fasta

>adk_1

GGGGAAAGGGACTCAGGCTCAGTTCATCATGGAGAAATATGGTATTCCGCAAATCTCCAC

TGGCGATATGCTGCGTGCTGCGGTCAAATCTGGCTCCGAGCTGGGTAAACAAGCAAAAGA

CATTATGGATGCTGGCAAACTGGTCACCGACGAACTGGTGATCGCGCTGGTTAAAGAGCG

CATTGCTCAGGAAGACTGCCGTAATGGTTTCCTGTTGGACGGCTTCCCGCGTACCATTCC

GCAGGCAGACGCGATGAAAGAAGCGGGCATCAATGTTGATTACGTTCTGGAATTCGACGT

ACCGGACGAACTGATTGTTGATCGTATCGTAGGCCGCCGCGTTCATGCGCCGTCTGGTCG

TGTTTATCACGTTAAATTCAATCCGCCGAAAGTAGAAGGCAAAGACGACGTTACCGGTGA

AGAACTGACTACCCGTAAAGACGATCAGGAAGAAACCGTACGTAAACGTCTGGTTGAATA

CCATCAGATGACTGCACCGCTGATCGGCTACTACTCCAAAGAAGCGGAAGCGGGTA

>adk_2

GGGGAAAGGGACTCAGGCTCAGTTCATCATGGAGAAATATGGTATTCCGCAAATCTCCAC

TGGCGATATGCTGCGTGCTGCGGTCAAATCTGGCTCCGAGCTGGGTAAACAAGCAAAAGA

CATTATGGATGCTGGCAAACTGGTTACCGACGAACTGGTGATCGCGCTGGTTAAAGGGCG

CATTGCTCAGGAAGACTGCCGTAATGGTTTCCTGTTGGACGGCTTCCCGCGTACCATTCC

GCAGGCAGACGCGATGAAAGAAGCGGGCATCAATGTTGATTACGTTCTGGAATTCGACGT

 

実行方法

ペアエンドfastqを指定する。

srst2 --input_pe ERR024070*.fastq.gz --output shigella1 --log --save_scores --mlst_db Escherichia_coli#1.fasta --mlst_definitions ecoli.txt --gene_db ARGannot.fasta

複数のファイルが出力される。 

f:id:kazumaxneo:20171218121404j:plain

結果

f:id:kazumaxneo:20190221134321j:plain

拡大

f:id:kazumaxneo:20190221134350j:plain

 

オーサーらの解説には、複数データを指定する方法、結果をヒートマップで出力するRのコード(論文のFig.8)も説明されています。

引用
SRST2: Rapid genomic surveillance for public health and hospital microbiology labs
Inouye M, Dashnow H, Raven LA, Schultz MB, Pope BJ, Tomita T, Zobel J, Holt KE

Genome Med. 2014 Nov 20;6(11):90