バクテリアの遺伝子配列を比較する LS-BSR

2021 1/18 わかりにくい説明を修正

　細菌単離株からの全ゲノム配列データが安価に入手できるようになるにつれ、配列データと生物学的観察結果を相関させる計算手法が必要とされている。ここでは、数百から数千の細菌ゲノムの遺伝的内容を迅速に比較し、調査した全ゲノムの全コーディング配列（CDS）の関連性を記述したマトリックスを返す large-scale BLAST score ratio (LS-BSR) パイプラインを紹介する。このマトリックスは、細菌ゲノム間の遺伝的関係を特定するために簡単に解析することができる。配列類似度によってペプチドをグループ化するパイプラインは公開されているが、LS-BSRで実施されているような迅速で大規模なフルゲノム比較解析を行うソフトウェアは他にない。

　この方法の有用性を実証するために、LS-BSRパイプラインを96のEscherichia coliおよびShigellaゲノム上でテストしたところ、16プロセッサを用いて163分で実行された。各CDSのBSR値は、相対的な関連性のレベルを示すもので、独立したコアゲノム一塩基多型（SNP）ベースの系統樹に基づいて各ゲノムにマッピングされた。この比較により、クレード特異的なCDSマーカーを同定し、古典的な大腸菌病原性バリアント（病原性バリアント）指定を区別する分子マーカーに基づいてLS-BSRパイプラインを検証した。スケーラビリティ試験の結果、LS-BSRパイプラインは、アライメント方法にもよるが、16プロセッサを用いて1,000個の大腸菌ゲノムを27～57時間で処理できることが実証された。

　LS-BSRはオープンソースのBSRアルゴリズムの並列実装であり、大量のゲノムの遺伝的内容を迅速に比較することができる。パイプラインの結果は、ユーザーが定義した系統群間で特定のマーカーを識別したり、細菌分離株間での遺伝情報の喪失および/または獲得を識別するために使用することができる。分類群特異的な遺伝マーカーは、その後、臨床診断に変換することができ、また、広く保存されている治療候補を特定するために使用することができる。

　LS-BSR法は、定義された遺伝子セットを使用するか、またはクエリゲノムのセットからCDSをProdigal (Hyatt et al., 2010)で予測して実行することができる。Prodigalを使用する場合、すべてのCDSは連結され、その後0.9（同一性のしきい値は、ユーザーによって変更することができる）のペアワイズ同一性でUSEARCH（エドガー、2010）を使用してデレプリケートされる。それぞれのユニークなCDSは、次にBioPython (www.biopython.org)で翻訳され、参照ビットスコアを計算するためにTBLASTN (Altschul et al., 1997)でそのヌクレオチド配列に対してアラインメントされる; BLASTNまたはBLAT (Kent, 2002)が呼び出された場合、ヌクレオチド配列はアラインメントされる。次に、各クエリ配列は、BLAT、BLASTN、またはTBLASTNを用いて各ゲノムに対してアラインメントされ、クエリビットスコアが集計される。BSR値は、クエリビットスコアを参照ビットスコアで割って算出し、0.0～1.0の間のBSR値が得られる（TBLASTNで得られたビットスコア値が変動するため、1.0よりも若干高い値が観測される）。LS-BSRパイプラインの結果には、調査した各ゲノムに固有のCDS名とBSR値を含むマトリックスが含まれている。少なくとも1つのゲノムに1つ以上の有意なBSR値を持つCDSもまた、出力で識別される。少なくとも1つのゲノムにおいて、1つの重複が他の重複と有意に異なるCDSについては、別のファイルが生成される；これらの領域はパラログを表している可能性があり、さらに詳細な調査が必要な場合がある。LS-BSRマトリックスが生成されると、結果はMultiple Experiment Viewer (MeV) (Saeed et al., 2006)やR (R Core Team, 2013)を使ってヒートマップやクラスタとして簡単に可視化することができる。Interactive Tree Of Life project (Bork et al., 2008)を使用して、LS-BSRの出力からヒートマップを生成し、ヒートマップデータを提供された系統図と相関させることもできる。LS-BSRには、CDSの有無に応じて設定されたBSRしきい値の範囲を使って、ユーザー定義のサブグループ間でCDSを迅速に比較するためのスクリプト（compare_BSR.py）が含まれている。識別されたCDSのアノテーションは、RAST (Aziz et al., 2008) や prokka (http://www.vicbioinformatics.com/software.prokka.shtml) などのツールを使用して適用することができる。LS-BSRのソースコード、ユニットテスト、テストデータは、GNU GPL v3ライセンスのもと、https://github.com/jasonsahl/LS-BSR で自由に入手できる。

インストール

依存

Minimum requirements, see manual.md for version information

• Python >2.7 and <=3.5 (higher versions still work but tests may fail)
• BioPython
• Prodigal - Required for de novo gene prediction only
• VSEARCH - Optional
• mmseqs2- Optional
• CD-HIT - Optional
• Blast+ - Optional
• Blat - Optional
• Diamond - Optional

Github

#仮想環境の作成と依存ツールの導入
mamba create -n ls_bsr python=3.5 -y
conda activate ls_bsr
mamba install -c bioconda blast vsearch cd-hit prodigal ucsc-blat diamond biopython mmseqs2

git clone https://github.com/jasonsahl/LS-BSR.git
cd LS-BSR/
python setup.py install

>python ls_bsr.py

$ python ls_bsr.py 

Must provide directory.

Usage: ls_bsr.py [options]

Options:

  --version             show program's version number and exit

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

  -d DIRECTORY, --directory=DIRECTORY

                        /path/to/fasta_directory [REQUIRED]

  -i ID, --identity=ID  clustering id threshold (0.0-1.0), defaults to 0.9

  -f FILTER, --filter=FILTER

                        to use blast filtering or not, default is F or filter,

                        change to T to turn off filtering

  -p PROCESSORS, --parallel_workers=PROCESSORS

                        How much work to do in parallel, defaults to 2

  -g GENES, --genes=GENES

                        predicted genes (nucleotide) to screen against

                        genomes, will not use prodigal, must end in fasta (nt)

                        or pep (aa)

  -c CLUSTER_METHOD, --cluster_method=CLUSTER_METHOD

                        Clustering method to use: choose from mmseqs, mmseqs-

                        lin, vsearch, cd-hit

  -b BLAST, --blast=BLAST

                        use tblastn, blastn, blastp, blastn-short, diamond, or

                        blat (nucleotide search only), default is tblastn

  -l LENGTH, --length=LENGTH

                        minimum BSR value to be called a duplicate, defaults

                        to 0.7

  -m MAX_PLOG, --max_plog=MAX_PLOG

                        maximum value to be called a paralog, defaults to 0.85

  -n MIN_HLOG, --min_hlog=MIN_HLOG

                        minimum BLAST ID to be called a homolog, defaults to

                        75

  -t F_PLOG, --f_plog=F_PLOG

                        filter ORFs with a paralog from BSR matrix? Default is

                        F, values can be T or F

  -k KEEP, --keep=KEEP  keep or remove temp files, choose from T or F,

                        defaults to F

  -s FILTER_PEPS, --filter_short_peps=FILTER_PEPS

                        remove short peptides, smaller than 50AA?  Defaults to

                        F

  -e FILTER_SCAFFOLDS, --filter_scaffolds=FILTER_SCAFFOLDS

                        Filter any contig that contains an N? Defaults to F

  -x PREFIX, --prefix=PREFIX

                        prefix for naming output files, defaults to time/date

  -y INTERGENICS, --intergenics=INTERGENICS

                        Incoporate intergenic regions? T or F; Defaults to F

  --ml=MIN_LEN, --min_len=MIN_LEN

                        value for 's' flag in cd-hit, defaults to '-i' flag in

                        LS-BSR (float)

  -z DUP_TOGGLE, --dup_toggle=DUP_TOGGLE

                        Perform duplicate searching? T or F; Defaults to F

テストラン

python ls_bsr.py --version
python tests/test_all_functions.py

実行方法

fastaf配列が置かれているディレクトリを指定する。ここではテストデータを使う。LS-BSR/test_data/に移動。ゲノムファイルの拡張子は.fastaを認識する。

cd LS-BSR/test_data/

test_data/genomesを使う。

移動したらランする。上の画像のゲノムとクエリの遺伝子を指定する。クエリの遺伝子ecoli_markers.fastaには５つの遺伝子がmulti-fasta形式で記載されている。

python ../ls_bsr.py -d genomes/ -g genes/ecoli_markers.fasta

#proteinを指定するなら、拡張子は.pepである必要がある
python ../ls_bsr.py -d genomes/ -g queryproteins.pep

• -g    predicted genes (nucleotide) to screen against genomes, will not use prodigal, must end in fasta (nt) or pep (aa)

出力

bsr_matrix.txtを見てみる。タブをコンマに置換してtty-table（紹介）で整形表示する。

> sed s/$'\t'/','/g bsr_matrix.txt |tty-table 

$ sed s/$'\t'/','/g bsr_matrix.txt |tty-table 

  ┌───────┬──────────────┬───────────────────┬────────────┬──────────────┐

  │       │ E2348_69_all │ O157_H7_sakai_all │ H10407_all │ SSON_046_all │

  ├───────┼──────────────┼───────────────────┼────────────┼──────────────┤

  │ IpaH3 │    0.0000    │      0.1034       │   0.0529   │    0.9970    │

  ├───────┼──────────────┼───────────────────┼────────────┼──────────────┤

  │  LT   │    0.0000    │      0.0000       │   1.0000   │    0.0000    │

  ├───────┼──────────────┼───────────────────┼────────────┼──────────────┤

  │  ST2  │    0.0000    │      0.0000       │   0.9342   │    0.0000    │

  ├───────┼──────────────┼───────────────────┼────────────┼──────────────┤

  │ bfpB  │    1.0000    │      0.0329       │   0.0000   │    0.0000    │

  ├───────┼──────────────┼───────────────────┼────────────┼──────────────┤

  │ stx2a │    0.0000    │      0.0000       │   0.0000   │    0.0000    │

  └───────┴──────────────┴───────────────────┴────────────┴──────────────┘

他にもランパラメータlogや遺伝子名のファイルが出力される。

クエリを指定しないと、de novoで遺伝子予測が行われ、それから翻訳されたタンパク質配列を使って比較が実行される。入力ゲノム数に比例して相応の時間がかかる。

python ../ls_bsr.py -d genomes -c mmseqs

• d    /path/to/fasta_directory [REQUIRED]
• -c    Clustering method to use: choose from mmseqs, mmseqs-
  lin, vsearch,

出力されるmatrixファイルの先頭

GET HOMOLOGUESのcompare_clusters.plでも似た出力を得られる（GET HOMOLOGUESでは閾値を満たしたコア遺伝子のみ対象）。 

ヒートマップで可視化する例は伊藤さんが書かれたQiitaの記事を参考にして下さい。丁寧に説明されています。

追記

パラログは除く。20並列。クエリはプロテイン。タンパク質配列ファイルの拡張子は.pepを認識する。

python ls_bsr.py -d genomes/ -g query.pep -s T -t T -l 0.7 -p 20 -b tblastn

• -t    filter ORFs with a paralog from BSR matrix? Default is F, values can be T or F
• -s    remove short peptides, smaller than 50AA?  Defaults to F
• -l     minimum BSR value to be called a duplicate, defaults to 0.7
• -m   maximum value to be called a paralog, defaults to 0.85
• -n    minimum BLAST ID to be called a homolog, defaults to 75
• -p    How much work to do in parallel, defaults to 2
• -b    use tblastn, blastn, blastp, blastn-short, diamond, or blat (nucleotide search only), default is tblastn
   
   

引用
The large-scale blast score ratio (LS-BSR) pipeline: a method to rapidly compare genetic content between bacterial genomes

Jason W Sahl, J Gregory Caporaso, David A Rasko, Paul Keim

PeerJ. 2014 Apr 1;2:e332

関連