macでインフォマティクス

macでインフォマティクス

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

従来型のMLSTタイピングを行う mlst

 

mlstは、ゲノム配列を入力として伝統的な7つの遺伝子座に対するタイピングを行う。PubMLSTのタイピングスキームに従う。

 

インストール

依存

  • Perl >= 5.26
  • NCBI BLAST+ blastn >= 2.9.0
  • You probably have blastn already installed already.
  • If you use Brew or Conda, this will install the blast package for you.
  • Perl modules: Moo,List::MoreUtils,JSON
  • Debian: sudo apt-get install libmoo-perl liblist-moreutils-perl libjson-perl
  • Redhat: sudo apt-get install perl-Moo perl-List-MoreUtils perl-JSON
  • Most Unix: sudo cpan Moo List::MoreUtils JSON
  • any2fasta
  • Converts sequence files to FASTA, even compressed ones

Github

#bioconda(link)
conda create -n mlst -y
conda activate mlst
conda install -c conda-forge -c bioconda -c defaults mlst -y

#homebrew
brew install brewsci/bio/mlst

 > mlst -h

$ mlst -h

SYNOPSIS

  Automatic MLST calling from assembled contigs

USAGE

  % mlst --list                                            # list known schemes

  % mlst [options] <contigs.{fasta,gbk,embl}[.gz]          # auto-detect scheme

  % mlst --scheme <scheme> <contigs.{fasta,gbk,embl}[.gz]> # force a scheme

GENERAL

  --help            This help

  --version         Print version and exit(default ON)

  --check           Just check dependencies and exit (default OFF)

  --quiet           Quiet - no stderr output (default OFF)

  --threads [N]     Number of BLAST threads (suggest GNU Parallel instead) (default '1')

  --debug           Verbose debug output to stderr (default OFF)

SCHEME

  --scheme [X]      Don't autodetect, force this scheme on all inputs (default '')

  --list            List available MLST scheme names (default OFF)

  --longlist        List allelles for all MLST schemes (default OFF)

  --exclude [X]     Ignore these schemes (comma sep. list) (default 'ecoli_2,abaumannii')

OUTPUT

  --csv             Output CSV instead of TSV (default OFF)

  --json [X]        Also write results to this file in JSON format (default '')

  --label [X]       Replace FILE with this name instead (default '')

  --nopath          Strip filename paths from FILE column (default OFF)

  --novel [X]       Save novel alleles to this FASTA file (default '')

  --legacy          Use old legacy output with allele header row (requires --scheme) (default OFF)

SCORING

  --minid [n.n]     DNA %identity of full allelle to consider 'similar' [~] (default '95')

  --mincov [n.n]    DNA %cov to report partial allele at all [?] (default '10')

  --minscore [n.n]  Minumum score out of 100 to match a scheme (when auto --scheme) (default '50')

PATHS

  --blastdb [X]     BLAST database (default '/Users/kazu/anaconda3/envs/mlst/db/blast/mlst.fa')

  --datadir [X]     PubMLST data (default '/Users/kazu/anaconda3/envs/mlst/db/pubmlst')

HOMEPAGE

  https://github.com/tseemann/mlst - Torsten Seemann

 

実行方法

ゲノム配列、またはcontig配列のfastagenbankファイルを指定する。

mlst input.fa

#複数
mlst *fasta

gzip, zip、 bzip2圧縮状態にも対応している。

 

Streptomyces atratusのstrain SCSIO ZH16 完全長ゲノムに対して適用してみた。

mlst Streptomyces_atratus_strain_SCSIO_ZH16.fasta

$ mlst Salmonella_enterica_subsp._enterica_serovar_Enteritidis_str._RM2968.fasta 

[13:34:55] This is mlst 2.19.0 running on darwin with Perl 5.026002

[13:34:55] Checking mlst dependencies:

[13:34:55] Found 'blastn' => /Users/kazu/anaconda3/envs/mlst/bin/blastn

[13:34:55] Found 'any2fasta' => /Users/kazu/anaconda3/envs/mlst/bin/any2fasta

[13:34:55] Found blastn: 2.9.0+ (002009)

[13:34:55] Excluding 2 schemes: abaumannii ecoli_2

[13:34:57] Found exact allele match ecoli.icd-40

[13:34:57] Found exact allele match senterica.thrA-11

[13:34:57] Found exact allele match senterica.sucA-6

[13:34:57] Found exact allele match senterica.hisD-7

[13:34:57] Found exact allele match senterica.dnaN-2

[13:34:57] Found exact allele match senterica.aroC-5

[13:34:57] Found exact allele match senterica.hemD-3

[13:34:57] Found exact allele match senterica.purE-6

Salmonella_enterica_subsp._enterica_serovar_Enteritidis_str._RM2968.fasta senterica 11 aroC(5) dnaN(2) hemD(3) hisD(7) purE(6) sucA(6) thrA(11)

[13:34:57] If you like MLST, you're absolutely going to love wgMLST!

[13:34:57] Done.

見えにくいが、以下の部分が結果。

Salmonella_enterica_subsp._enterica_serovar_Enteritidis_str._RM2968.fasta senterica 11 aroC(5) dnaN(2) hemD(3) hisD(7) purE(6) sucA(6) thrA(11)

オートモードで実行したが、sentericaと検出され、7つの遺伝子座とタイピング結果が示されている。

 

スキームを自動指定せず、ユーザーが指定することも可能。 利用可能なスキームを表示。

 mlst --list

$ mlst --list 

cronobacter otsutsugamushi wolbachia sinorhizobium sthermophilus_2 ssuis csputorum cconcisus tpallidum neisseria ecloacae bcereus chyointestinalis senterica borrelia spyogenes mmassiliense vvulnificus saureus shaemolyticus pacnes pmultocida_multihost pfluorescens ganatis bpseudomallei sbsec bsubtilis kaerogenes vibrio aeromonas kkingae slugdunensis cmaltaromaticum bhampsonii bpilosicoli ypseudotuberculosis sdysgalactiae kpneumoniae vtapetis cbotulinum spneumoniae sthermophilus sepidermidis vparahaemolyticus clari yersinia mpneumoniae smaltophilia psalmonis leptospira mhaemolytica pgingivalis ecoli_2 hinfluenzae mhyorhinis bbacilliformis efaecalis cdiphtheriae abaumannii cdifficile fpsychrophilum arcobacter pputida cupsaliensis mcatarrhalis taylorella chelveticus mbovis pdamselae blicheniformis scanis edwardsiella abaumannii_2 plarvae mycobacteria leptospira_3 streptomyces koxytoca ecoli lsalivarius soralis hcinaedi mhyopneumoniae ppentosaceus mabscessus bhyodysenteriae xfastidiosa leptospira_2 bcc pmultocida_rirdc vcholerae2 bintermedia yruckeri suberis tenacibaculum bhenselae magalactiae brucella clanienae msynoviae vcholerae cinsulaenigrae dnodosus chlamydiales lmonocytogenes mflocculare miowae shominis szooepidemicus sgallolyticus hparasuis mcanis ranatipestifer campylobacter mcaseolyticus efaecium orhinotracheale brachyspira csepticum bwashoensis cfreundii paeruginosa sagalactiae rhodococcus achromobacter mplutonius liberibacter cfetus hsuis aphagocytophilum hpylori ureaplasma bordetella spseudintermedius

 

スキームを指定して実行。

mlst --scheme vibrio Vibrio*.gbk.gz

引用

Seemann T, mlst Github https://github.com/tseemann/mlst

 

"This publication made use of the PubMLST website (https://pubmlst.org/) developed by Keith Jolley (Jolley & Maiden 2010, BMC Bioinformatics, 11:595) and sited at the University of Oxford. The development of that website was funded by the Wellcome Trust".

 

関連


 

 

 

特定の領域由来のロングリードを高速選抜する selectION

 

 

SelectION: Identification of predefined genomic regions in large nanopore DNA

f:id:kazumaxneo:20200530135435p:plain

 London Calling 2017

 

インストール

ubuntu18.04LTSでテストした。

ビルド依存

requires gcc > 5 and the following libraries:

  • boost filesystem
  • boost program_options
  • boost system
  • libhdf5

Boost is available through standard package sources. Libhdf5 is downloaded and build by the install script. 

本体 Github

git clone https://github.com/PayGiesselmann/selectION
cd selectION
mkdir build
cd build
cmake -DCMAKE_BUILD_TYPE=Release ..
make
make install

> selection index-h

# selection index-h

Program: SelectION

Usage: selection <command> [options]

Commands: index : Build FM-Index for reference sequence

scan : Scan input for reads matching specified positions

selection scan -h

# selection scan -h

Usage: selection scan [options] <db.prefix> <input.fq> <outputDir>

Read selection options:

  -f [ --filter ] arg        Input selection filter

  -s [ --sam ] arg           Write pseudo alignment in sam format to file

  -t [ --threads ] arg (=1)  Number of threads

  -q [ --quality ] arg (=20) Quality threshold for filtered reads

  --scanPrefix arg (=1000)   Prefix of read to use for alignment

 

 

 実行方法

1、index

参照ゲノムのインデックスを作成する。

selection index -t 8 ref.fa

 

2、Scan

input.fq内のすべてのポジションを推定し、その結果をout.samに書き込む。

selection scan -t 8 ref.fa input.fq ./ --sam ./out.sam

フルのアラインメントを しないため動作は非常に高速。1GB程度のファイルなら10秒程度でsam出力する。

 

動画を見てもらえば分かりますがまだ開発中です。更新がありましたら追記します。ONT fast5からの直接選抜などが予定されているようです。 

引用

GitHub - giesselmann/selectION: Rapid linking of long reads to a reference genome

 

 

 

メタゲノムのビニング後の解析を行う自動化されたパイプライン MetaSanity

2020 5/29 構成を修正、タイトル変更

 

 マイクロバイオーム研究の重要性はますます一般的になっており、さまざまな生態系(例:海洋、構築、宿主関連など)を理解するために不可欠である。研究者は、微生物ゲノムの分析のため、高度に再現性と品質を担保して解析パイプラインを実行できなければならない。 MetaSanityは、既存の11の広く使用されているゲノム評価とアノテーションスイートの分析を単一の配布可能なワークフローに組み込み、柔軟で拡張可能なデータ分析パイプラインを可能にして、微生物研究者の作業負荷を軽減する。 MetaSanityは、個別の再現可能なワークフローを提供するように設計されており、研究者のニーズによって、(1)推定される系統発生アサインを提供しながら微生物ゲノムの全体的な品質を決定でき、(2)適合する特異性の程度が異なる構造的および機能的な遺伝子アノテーションを割り当てることができる。このソフトウェアスイートは、いくつかのツールの結果を組み合わせて、全体的な代謝機能とペプチダーゼと炭水化物活性酵素の推定細胞外局在化に関する幅広い洞察を提供する。重要なことは、このソフトウェアは、関連するすべての出力をSQLデータベースに保存することにより、「ビッグデータ」分析の組み込み最適化を提供し、ユーザーが研究に最も影響を与える要素のすべての結果を照会できるようにする事である。MetaSanityは、GNU General Public License v.3.0に基づいて提供され、https://github.com/cjneely10/MetaSanityからダウンロードできる。このアプリケーションはPython3 / CythonおよびC ++で実装されており、Dockerイメージとして配布される。 

 

f:id:kazumaxneo:20200222111255p:plain

MetaSanity pipeline schema.  Preprintより転載。

 

wiki

https://github.com/cjneely10/MetaSanity/wiki/2-Installation

 

2019

 

インストール

ubuntu18.04LTSでテストした。

本体 Github

conda create -n metasanity -y
conda activate metasanity
#インストールスクリプトのダウンロード
wget
https://raw.githubusercontent.com/cjneely10/MetaSanity/master/install.py -O install.py

#実行。docker image(link)がダウンロードされる。root権限でなければsudo実行する。数GBあるので注意。
python3 install.py
=> docker imagesができ、さらにMetaSanity/ができる

./MetaSanity.py -h

$ ./MetaSanity.py -h

usage: MetaSanity.py [-h] -d DIRECTORY -c CONFIG_FILE [-a]

                     [-o OUTPUT_DIRECTORY] [-b BIOMETADB_PROJECT]

                     [-t TYPE_FILE] [-p] [-z] [-v]

                     program

 

MetaSanity: Run meta/genomes evaluation and annotation pipelines

 

Available Programs:

 

FuncSanity: Runs gene callers and annotation programs on MAGs

(Flags:  --directory --config_file --cancel_autocommit --output_directory --biometadb_project --type_file --prokka --reevaluate_quality)

PhyloSanity: Evaluates completion, contamination, and redundancy of MAGs

(Flags:  --directory --config_file --cancel_autocommit --output_directory --biometadb_project)

 

positional arguments:

  program               Program to run

 

optional arguments:

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

  -d DIRECTORY, --directory DIRECTORY

                        Directory name containing genomes

  -c CONFIG_FILE, --config_file CONFIG_FILE

                        Config file name

  -a, --cancel_autocommit

                        Cancel commit to database

  -o OUTPUT_DIRECTORY, --output_directory OUTPUT_DIRECTORY

                        Output directory prefix, default out

  -b BIOMETADB_PROJECT, --biometadb_project BIOMETADB_PROJECT

                        BioMetaDB_project name (updates values of existing database)

  -t TYPE_FILE, --type_file TYPE_FILE

                        type_file name formatted as 'file_name.fna\t[Archaea/Bacteria]\t[gram+/gram-]\n'

  -p, --prokka          Use PROKKA gene calls instead of prodigal search

  -z, --remove_intermediates

                        Remove intermediary genome directories, default: False

  -v, --version         Print version of MetaSanity and exit

(metasanity) kazu@kazu:~/Documents/MetaSanity$ 

 

 

ライセンスの関係でSignalP4.1 and RNAmmer1.2は別に導入する必要がある。必須ではないため、ここでは説明しない。公式wikiを参照。続いてデータベースのダウンロードを行う。

 

データベースの準備

ダウンロードスクリプトが用意されている。自動作成されたMetaSanity/にてスクリプトを実行する。

cd MetaSanity/
python3 download-data.py

以下のデータベースがダウンロードされる。 データベースサイズは46GBほどあるので空き容量に注意する。

  • GTDB-Tk release 89.0
  • CheckM 2015_01_16
  • KoFamScan (most recent release as of Feb 2020)
  • dbCAN HMMdb-V8
  • virsorter-data-v2 

 

  

実行方法

大きく分けて2つのモードがある。PhyloSanityモードでは、CheckMを用いてMAGの完全性と汚染度を評価し、FastANIを用いて冗長性を評価する。オプションでGTDB-Tkを用いて系統を予測する。FuncSanityモードでは、PROKKA、KEGG、InterProScan、炭水化物酵素(CAZy)データベース、MEROPS データベースを含む機能アノテーション パイプラインをカスタマイズしてランできる。オプションで、PSORTb と SignalP4.1 を使用して、予測された外膜および細胞外タンパク質をターゲットに 、MEROPS の一致をフィルタリングできる。KEGG の結果は KEGG-Decoder で処理され、代謝機能とパスウェイの視覚的なヒートマップが提供される。いずれのランにもconfigファイルを使う。

 

1、準備

作業ディレクトリを作り、そこにマスターのコンフィグファイルをコピーする。

#1 作業ディレクトリの作成
mkdir ExampleRun && cd ExampleRun && mkdir genomes

#2 テストデータをgenomes/にコピーする。
cp <path>/<to>/MetaSanity/ExampleSet/TOBG-* genomes/

#3 configファイルをカレントにコピーする。PhyloSanity.iniとFuncSanity.iniを両方ランするので、2つともコピー。
cp <path>/<to>MetaSanity/Config/Docker/PhyloSanity.ini .
cp <path>/<to>MetaSanity/Config/Docker/FuncSanity.ini .

 PhyloSanity.ini(configファイル)

f:id:kazumaxneo:20200528225256p:plain

パスは仮想環境のものなので修正する必要はない。パラメータのみ必要であれば変更する。利用可能なすべてのフラグはComplete-PhyloSanity.iniを参照。

 

PhyloSanity.ini

f:id:kazumaxneo:20200528230515p:plain

f:id:kazumaxneo:20200528230519p:plain

f:id:kazumaxneo:20200528230518p:plain

利用可能なすべてのフラグはComplete-FuncSanity.iniをを参照。

 

 

2、ラン

PhyloSanity

PhyloSanityのコンフィグを指定して実行する。dockerのランにroot権限が必要ならsudo実行する。

MetaSanity -c PhyloSanity.ini -d genomes/ PhyloSanity

出力ディレクトリ(out/)

f:id:kazumaxneo:20200528225803p:plain

out/evaluation.tsv - GTDBtk、checkm、fastANIなどの結果が統合される。

f:id:kazumaxneo:20200528225809p:plain



FuncSanity

次はFuncSanityを動かす。FuncSanityのコンフィグを指定して実行する。

MetaSanity -c FuncSanity.ini -d genomes/ FuncSanity

グラム陰性のバクテリア向けになっている。それ以外のゲノムを使う場合はヘルパースクリプトを使う(wiki参照)。

 

out/に結果は追加される。サンプルそれぞれのTSVファイルと、全サンプル統合したTSVファイルが出力される。

f:id:kazumaxneo:20200529124051p:plain

 

KEGG

combined_results

f:id:kazumaxneo:20200529123545p:plain

biodata_results

f:id:kazumaxneo:20200529123656p:plain

 

function.TSV 

f:id:kazumaxneo:20200529123656p:plain

f:id:kazumaxneo:20200529123751p:plain

 

f:id:kazumaxneo:20200529123751p:plain

f:id:kazumaxneo:20200529123755p:plain

f:id:kazumaxneo:20200529123804p:plain

f:id:kazumaxneo:20200529123807p:plain

f:id:kazumaxneo:20200529123810p:plain

f:id:kazumaxneo:20200529123819p:plain

(横に長いので以下省略)

 ProdigalとPROKKAに続き、InterProScan、KEGG KoFamScanとKEGG-Decoder、proteinの細胞内局在、シグナルペプチド予測と切断部位予測、peptidasesとExtracellular Peptidase予測、catalytic and carbohydrate-binding module予測(炭水化物の合成、修飾、分解に関与する炭水化物活性酵素(CAZymes)と関連する炭水化物結合モジュールの予測)、virus予測が行われる。

 

他に結果を管理するリレーショナルデータベースBioMetaDBがある。

引用

MetaSanity: An integrated microbial genome evaluation and annotation pipeline
Christopher J Neely, Elaina D Graham, Benjamin J Tully
Bioinformatics, Published: 19 May 2020

 

MetaSanity: An integrated, customizable microbial genome evaluation and annotation pipeline.

Neely C., Graham E. D., Tully, B. J.

BioRxiv, Posted October 08, 2019.

 

関連


 

 

 

 

 

 

 

 

 

 

 

 

パンゲノム解析を行う GET_HOMOLOGUES

2020 5/28 追記

 

 GenBank のような公開データベースに登録されているゲノムの数が増え続けていることから、種の遺伝子レパートリーを比較するためのツールが開発されている。このような比較には、種分化イベントの後に共通の祖先から発散したと仮定され、パラログよりも生物間で機能が保存されている可能性が高いオルソログ遺伝子の同定が含まれる(ref.1)。このため、オルソログはゲノムアノテーションや進化研究の重要な要素となっている(ref.2,3)。生命の他のどの領域よりも高速に配列決定されている細菌間では(ref.4)、直交配列を検出するための一般的なヒューリスティックな方法は、単にBLASTの相互ヒットを探すことであり(ref.5, 6)、このタスクのために様々なソフトウェアが利用可能である(ref.7)。これらのツールと増加するゲノム配列を組み合わせることで、最近の研究では、細菌のゲノムは、関心のあるグループのすべての分離株で共有される遺伝子(コアゲノム)と、菌株特異的/部分的に共有される遺伝子を含むモザイクであることを示唆する証拠がいくつか提供されている(ref.8)。コアゲノムとグループ内の残りの遺伝子の合計がパンゲノムと定義される(ref.9)。

 ここでは、GNU General Public Licenseの下でリリースされているオープンソースのソフトウェアパッケージであるGET_HOMOLOGUESを紹介する。このソフトウェアは、Linux/Mac OS Xコンピュータシステム上で、系統的に異なる距離にある細菌株のパンゲノムおよび比較ゲノム解析のために特別に設計され、テストされている。このソフトウェアはいくつかの点でユニークである。ゲノムデータのダウンロード、ユーザーが選択した配列特徴の抽出、BLASTおよびHMMERジョブの実行、結果の索引付け、クラスタリング、および解析を含む、完全に自動化された高度にカスタマイズ可能な解析パイプラインを実装している。最新のマルチプロセッサアーキテクチャやコンピュータクラスタを利用して、時間のかかるBLASTやHMMERジョブを並列化することができる。Berkeley DBを使用して一時的なデータをディスクに書き込むか、双方向ベストヒット(BDBH)アルゴリズムヒューリスティックバージョンを呼び出すことで、適度なマシン(8GB RAM未満)で大規模なデータセット(例えば、101のEscherichia coliゲノムを解析している)を扱うことができる。サポートされている配列クラスタリングアルゴリズムの組み合わせによって復元されたコンセンサスクラスタの計算を含む、遺伝子ファミリーの解析と生成を容易にするための補助スクリプトが統合されている。また、R関数を呼び出して統計解析を行い、コアプロットやパンゲノムプロットなどの結果をグラフィカルに表示するためのスクリプトも用意されている。また、多様な比較ゲノム解析も可能である。最後に、インストール作業を簡単にするためのインストールスクリプトが用意されており、また、ハンズオンチュートリアルを含む非常に詳細なマニュアルも用意されているため、使い勝手の良いパッケージとなっている。

 ここでは、公開されている完全なゲノムの中からオルソログを同定するデータベースであるOMA(Orthologous Matrix)の最新版からダウンロードした50のストレプトコッカスゲノムのセットを解析することで、これらの機能の一部を紹介する(ref.10)。著者らはいくつかの理由からこの属を選んだ。この属は非常に高いレベルのゲノム可塑性を示している(ref.11)。最初のパンゲノミクス解析は、Tettelinらの先駆的な研究であるStreptococcus agalactiae (ref.12)で行われ、その後、ヒトの主要な病原体であるS. pyogenes (ref.13)やS. pneumoniae (ref.14)を含むこの属の多様な種について、非常に詳細な比較ゲノミクス研究が行われており、StreptococcusはGET_HOMOLOGUESソフトウェアの優れたテストケースとなっている。

 

 このソフトウェアは、BLAST+(ref.16)とOrthoMCL(バージョン1.4)のコードベース(ref.17)の上に構築されており、3つの一般的なシーケンスクラスタリングアルゴリズムをサポートしている。OrthoMCL(OMCL)、COGtriangles(ref.18)、および双方向ベストヒット(BDBH)アルゴリズムの著者ら自身の実装をサポートしている(論文補足資料の図S1参照)。これらのアプローチは、戦略は異なるがBLASTの双方向ベストヒットを証拠として直交配列を呼び(ref.19)、同一ゲノム内でのベストヒットを持つ遺伝子、すなわち最近のinparalogues(ref.20)を区別している。さらに、タンパク質ドメインのPfamアノテーションを容易にするためにHMMER(http://hmmer.org)が統合されており(ref.21)、オルソログの割り当てを混乱させる可能性のある異なるドメイン構造を持つ配列を含むクラスタをフィルタリングして除外することができる。入力ファイルが GenBank 形式の場合、ヌクレオチド配列とアミノ酸配列の両方のクラスタが生成され、必要に応じて、オルソログ遺伝子に挟まれたオルソログ遺伝子間領域を抽出することもできる。また、GenBankファイル内のゲノム座標は、少なくとも1つの保存された近傍領域を有する同調性クラスタの選択に利用することができる(論文補足資料の図S2参照)。最後に、GenBankファイルには様々なDNAフィーチャーが含まれているため、例えばtRNA遺伝子に着目して検索することも可能である。この場合、デフォルトのBLASTPジョブの代わりにBLASTN検索が実行される。
 メインのPerlスクリプトget_homologues.plに加えて、このソフトウェアには、その後の解析に役立ついくつかの補助スクリプトがバンドルされている。例えば、複数のアルゴリズム(COGtrianglesやOMCLなど)によって生成されたインターセクションクラスタは、compare_clusters.plで簡単に選択することができ、ユーザーはコンセンサスクラスタのみを扱うことができる。さらに、パンゲノムとコアゲノムをプロットするためのスクリプト plot_pancore_matrix.pl が提供されている。plot_pancore_matrix.pl により、Tettelinら(9)とWillenbrockら(ref.22)の指数モデルをフィッティングしてコアとパンゲノムのサイズを推定できる。本論文で図示されているように、PHYLIPパッケージ(redf.23)のPARSプログラムを用いて、パンゲノミックマトリクスを表形式とPHYLIP形式で提供しており、 parsimony criterionの下でパンゲノミックツリーを自動生成することができる。parse_pangenome_matrix.plスクリプトは、系統特異的な遺伝子ファミリーや拡張の同定に焦点を当てた比較ゲノミクスや、コアゲノム、クラウドゲノム、シェルゲノムのコンパートメントの計算やグラフ化に役立つ。コア(core)は考慮されているすべてのゲノム/taxaに含まれる遺伝子、ソフトコア(soft core)はKaasと共同研究者の研究(25)で記述されているように考慮されているゲノム/taxaの95%に含まれる遺伝子、クラウド(cloud)は少数のゲノム/Ttaxaにのみ存在する遺伝子(カットオフは、最もぽピュレーションの多い非コアクラスタクラスとそのすぐ隣のクラスとして定義されている)、シェル(shell)は残りの遺伝子で、いくつかのゲノム/taxaに存在している遺伝子と定義される。このスクリプトはまた、Snipenら(ref.26)の二項混合モデルの下でコアおよびパンゲノムサイズの推定値を計算する。ソースのゲノムデータがGenBank形式で提供されている場合、parse_pangenome_matrix.plは、その系統から選択されたリファレンスゲノムの線形化された遺伝地図上にクレード特異的な遺伝子をプロットするように動かすことができる。compare_clusters.plで計算されたパンゲノミックツリーは、parse_pangenome_matrix.plで比較するグループのメンバーを選択するのに便利である。

 

Manual

 

GET_HOMOLOGUES

 

インストール

リリースからmacosのビルドをダウンロードし、macos10.14でテストした。

本体 Github

リリースから、各プラットフォーム向けにGET_HOMOLOGUES-ESTとGET_HOMOLOGUESのバイナリがダウンロードできる。その後、データベースをダウンロードしてインストールする。

cd get_homologues-macosx-20200226/
./install.pl

./get_homologues.pl

$ ./get_homologues.pl

 

usage: ./get_homologues.pl [options]

 

-h this message

-v print version, credits and checks installation

-d directory with input FASTA files ( .faa / .fna ),           (overrides -i,

   GenBank files ( .gbk ), 1 per genome, or a subdirectory      use of pre-clustered sequences

   ( subdir.clusters / subdir_ ) with pre-clustered sequences   ignores -c, -g)

   ( .faa / .fna ); allows for new files to be added later;    

   creates output folder named 'directory_homologues'

-i input amino acid FASTA file with [taxon names] in headers,  (required unless -d is set)

   creates output folder named 'file_homologues'

 

Optional parameters:

-o only run BLAST/Pfam searches and exit                       (useful to pre-compute searches)

-c report genome composition analysis                          (follows order in -I file if enforced,

                                                                ignores -r,-t,-e)

-R set random seed for genome composition analysis             (optional, requires -c, example -R 1234,

                                                                required for mixing -c with -c -a runs)

-s save memory by using BerkeleyDB; default parsing stores

   sequence hits in RAM

-m runmode [local|cluster]                                     (default local)

-n nb of threads for BLAST/HMMER/MCL in 'local' runmode        (default=2)

-I file with .faa/.gbk files in -d to be included              (takes all by default, requires -d)

 

Algorithms instead of default bidirectional best-hits (BDBH):

-G use COGtriangle algorithm (COGS, PubMed=20439257)           (requires 3+ genomes|taxa)

-M use orthoMCL algorithm (OMCL, PubMed=12952885)

 

Options that control sequence similarity searches:

-X use diamond instead of blastp                               (optional, set threads with -n)

-C min %coverage in BLAST pairwise alignments                  (range [1-100],default=75)

-E max E-value                                                 (default=1e-05,max=0.01)

-D require equal Pfam domain composition                       (best with -m cluster or -n threads)

   when defining similarity-based orthology

-S min %sequence identity in BLAST query/subj pairs            (range [1-100],default=1 [BDBH|OMCL])

-N min BLAST neighborhood correlation PubMed=18475320          (range [0,1],default=0 [BDBH|OMCL])

-b compile core-genome with minimum BLAST searches             (ignores -c [BDBH])

 

Options that control clustering:

-t report sequence clusters including at least t taxa          (default t=numberOfTaxa,

                                                                t=0 reports all clusters [OMCL|COGS])

-a report clusters of sequence features in GenBank files       (requires -d and .gbk files,

   instead of default 'CDS' GenBank features                    example -a 'tRNA,rRNA',

                                                                NOTE: uses blastn instead of blastp,

                                                                ignores -g,-D)

-g report clusters of intergenic sequences flanked by ORFs     (requires -d and .gbk files)

   in addition to default 'CDS' clusters

-f filter by %length difference within clusters                (range [1-100], by default sequence

                                                                length is not checked)

-r reference proteome .faa/.gbk file                           (by default takes file with

                                                                least sequences; with BDBH sets

                                                                first taxa to start adding genes)

-e exclude clusters with inparalogues                          (by default inparalogues are

                                                                included)

-x allow sequences in multiple COG clusters                    (by default sequences are allocated

                                                                to single clusters [COGS])

-F orthoMCL inflation value                                    (range [1-5], default=1.5 [OMCL])

-A calculate average identity of clustered sequences,          (optional, creates tab-separated matrix,

by default uses blastp results but can use blastn with -a      recommended with -t 0 [OMCL|COGS])

-P calculate percentage of conserved proteins (POCP),          (optional, creates tab-separated matrix,

by default uses blastp results but can use blastn with -a      recommended with -t 0 [OMCL|COGS])

-z add soft-core to genome composition analysis                (optional, requires -c [OMCL|COGS])

 

This program uses BLAST (and optionally HMMER/Pfam) to define clusters of 'orthologous'

genomic sequences and pan/core-genome gene sets. Several algorithms are available

and search parameters are customizable. It is designed to process (in a HPC computer

cluster) files contained in a directory (-d), so that new .faa/.gbk files can be added

while conserving previous BLAST results. In general the program will try to re-use

previous results when run with the same input directory.

 

./compare_clusters.pl

$ ./compare_clusters.pl

 

[options]:

-h  this message

-d  comma-separated names of cluster directories                  (min 1 required, example -d dir1,dir2)

-o  output directory                                              (required, intersection cluster files are copied there)

-n  use nucleotide sequence .fna clusters                         (optional, uses .faa protein sequences by default)

-r  take first cluster dir as reference set, which might contain  (optional, by default cluster dirs are expected

    a single representative sequence per cluster                   to be derived from the same taxa; overrides -t,-I)

-s  use only clusters with syntenic genes                         (optional, parses neighbours in FASTA headers)

-t  use only clusters with single-copy orthologues from taxa >= t (optional, default takes all intersection clusters; example -t 10)

-I  produce clusters with single-copy seqs from ALL taxa in file  (optional, example -I include_list; overrides -t)

-m  produce intersection pangenome matrices                       (optional, ideally expects cluster directories generated

                                                                   with get_homologues.pl -t 0)

-x  produce cluster report in OrthoXML format                     (optional)

-T  produce parsimony-based pangenomic tree                       (optional, requires -m)

 

./plot_matrix_heatmap.sh

$ ./plot_matrix_heatmap.sh

# ERROR: no input tab file defined!

    

    USAGE synopsis for: [plot_matrix_heatmap.sh v.v1.0.4_31Jan18]:

       plot_matrix_heatmap.sh <string (tab-delimited matrix file)> [-t|-m|-v|-o|-p|-H|-W|-C|-r|-k|-N|-M]

    

    REQUIRED

       -i <string> *Avg_identity.tab file     

       

    OPTIONAL:

     * filter out excessive redundancy in the tab-delimited ANI matrix file

       -c <float> (maximum) identity cut-off value (e.g. 97.3) [def: 100]

       -d <int> maximum number of decimals in matrix display [0-2; def: 0]

 

     * tweak the graphical output:

       -a <integer> angle to rotate node angles           [def 45]

       -t <string> text for plot title                    [def input_tab_file_name]

       -m <integer> margins_horizontal                    [def 18]

       -v <integer> margins_vertical                      [def 18]

       -o <string> output file format                     [def svg]

       -p <integer> points for plotting device            [def 15]    

       -H <integer> output device height                  [def 10]    

       -W <integer> output device width                   [def 15]     

       -C <flag> do not reorder clusters                  [def reorder clusters and plot dendrogram]

       -r <flag> remove column names and cell contents    [def names and cell contents are printed]

       -k <string> text for scale X-axis                  [def "Value"]

       -b <integer> right border (margin) in dendrograms  [def 10]

       -X <float> character expansion factor              [def 1.0]

     * Goodness of clustering

       -K <integer> max. number of clusters               [def: ; NOTE: 2 <= K <= n-1 ]

 

    RUN NJ ANALYSIS using ANI matrix (average identity matrix) generated by get_homologues.pl -A -t 0 -M|G

       -N <flag> will compute a distance matrix (ANdist) from the input similarity matrix

                 and use the former to construct a NJ tree and save it as newick string,

as well as a complete linkage dendrogram for convenient cluster delimitations

at ANdist cutoffs of 6,5 and 4, which correspond to ANI values of 94%, 95% and 96%, respectively

Do not use with a binary presence-absence matrix

    

    SUBSET MATRIX WITH REGULAR EXPRESSIONS

       -x <string> regex, like: 'species1|species2'       [def ]

    

    MORE DETAILED HELP:

       -M <flag> prints gplot installation instructions and further usage information

       

    EXAMPLE:

      plot_matrix_heatmap.sh -i Avg_identity.tab -c 99.5 -d 1 -t "Genus X ANIb (OMCL all clusters)" -N -o pdf -m 22 -v 22 -p 20 -H 20 -W 30 -x 'sp1|sp2' -a 45 -b 12 -X 0.9 -K 20

 

    #------------------------------------------------------------------------------------------------------------------

    AIM: Plot ordered heatmaps with row and col. dendrogram, from squared numeric (distance or presence-absence) matrix,

         like the *Avg_identity.tab matrices produced by get_homologues.pl with the -A flag

          

    OUTPUT:  svg|pdf files

    

    DEPENDENCIES:

         R packages ape and gplots. Run plot_matrix_heatmap.sh -M for installation instructions.

    

    NOTES:

      1. To improve the display, shorten the genome/taxon names as much as possible!

      2. Use -m and -v to control horizontal and vertial margins to the plot

         to fit the taxon labels, and -X character_expansion_factor to control font size

    #------------------------------------------------------------------------------------------------------------------

 

 

# Run check_dependencies() ... looks good: R is installed.

>./plot_pancore_matrix.pl

$ ./plot_pancore_matrix.pl

 

[options]: 

-h  this message

-i  input .tab data file                   (required, generated by get_homologues.pl)

-f  type of fit [pan|core_Tettelin|core_Willenbrock|core_both]

                                           (default: core_Tettelin, PubMed:16172379,18088402)

-F  font scale for fitted formulas         (optional, default=0.8)

-a  save snapshots for animations in dir   (optional, example: -a animation)

 

 

実行方法

解析するゲノムのGenbankファイル(*1)またはアミノ酸fastaファイルを集め、そのディレクトリを指定する。"-d"で指定するgenbankファイルは、gzip圧縮した状態で提供してもよい。

#default settings
perl get_homologues.pl -d genbank_dir/ -n 12

#orthoMCL使用。inparalogue含めず、10genome中9ゲノムヒット
perl get_homologues.pl -d genbank_dir/ -M -e -t 9 -n 12
  • -n   number of threads for BLAST/HMMER/MCL in 'local' runmode (default=2)
  • -G   use COGtriangle algorithm (COGS, PubMed=20439257) (requires 3+ genomes|taxa)
  • -M   use orthoMCL algorithm (OMCL, PubMed=12952885) Options that control sequence similarity searches:
  • -X   use diamond instead of blastp  (optional, set threads with -n)
  • -e    exclude clusters with inparalogues (by default inparalogues are included)
  • -t    report sequence clusters including at least t taxa (default t=numberOfTaxa,
    t=0 reports all clusters [OMCL|COGS])

 

 

テストデータを用いた解析の流れ

#1 download
wget https://github.com/eead-csic-compbio/get_homologues/releases/download/v3.3.2/get_homologues-x86_64-20200226.tgz
tar -xvf get_homologues-x86_64-20200226.tgz && rm get_homologues-x86_64-20200226.tgz
cd get_homologues-x86_64-20200226/

#2-1 run get_homologues(BDBH algorithm)
#4genomeのアミノ酸faaファイルが収められているsample_buch_fastaを対象にget_homologuesをラン
./get_homologues.pl -n 12 -d sample_buch_fasta .
#run後、カレントにsample_buch_fasta_homologues/ができる。

#2-2 run get_homologues(COG triangles algorithm)
./get_homologues.pl -n 12 -d sample_buch_fasta -G

#2-3 run get_homologues(OMCL algorithm)
./get_homologues.pl -n 12 -d sample_buch_fasta -M

#3 venn diagram(#4のため-Tと-mも指定)
./compare_clusters.pl -o sample_intersection -T -m -d sample_buch_fasta_homologues_BDBH/BuchaphCc_f0_alltaxa_algBDBH_e0_,sample_buch_fasta_homologues_COG/BuchaphCc_f0_alltaxa_algCOG_e0_,sample_buch_fasta_homologues_OMCL/BuchaphCc_f0_alltaxa_algOMCL_e0

#3の出力

f:id:kazumaxneo:20200525182927p:plain

 

#4 heatmap (Rのライブラリが必要)
install.packages("gplots")
install.packages("dendextend")
install.packages("factoextra")
install.packages("ape")
#run
./plot_matrix_heatmap.sh -i sample_intersection/pangenome_matrix_t0.tab -o pdf \
-r -H 8 -W 14 -m 28 -t "sample pangenome (clusters=180)" -k "genes per cluster"

#4の出力(別のサンプルデータ)

f:id:kazumaxneo:20200527230612p:plain

#5 cloud, shell、 core遺伝子の計算
./parse_pangenome_matrix.pl -m sample_intersection/pangenome_matrix_t0.tab -s

#5の出力

f:id:kazumaxneo:20200527231231p:plain

f:id:kazumaxneo:20200527231256p:plain

 

 

#6 
./get_homologues.pl -d sample_buch_fasta -c -n 12

#7
./plot_pancore_matrix.pl -i sample_buch_fasta_homologues/core_genome_algBDBH.tab -f core_Tettelin
  • -c   report genome composition analysis (follows order in -I file if enforced,
    ignores -r,-t,-e)

f:id:kazumaxneo:20200527231654p:plain


#ANI推定

./get_homologues.pl -d sample_plasmids_gbk -a 'CDS' -A -t 0 -M

cd sample_plasmids_gbk_homologues
../plot_matrix_heatmap.sh -i EcoliST131plasmidpKC394_f0_0taxa_CDS_algOMCL_e0_Avg_identity.tab \
-d 2 -t "CDS ANI matrix with 2 decimals"
  • -A   calculate average identity of clustered sequences, (optional, creates tab-separated matrix, by default uses blastp results but can use blastn with -a      recommended with -t 0 [OMCL|COGS])
  • -a   report clusters of sequence features in GenBank files  (requires -d and .gbk files, instead of default 'CDS' GenBank features                                               example -a 'tRNA,rRNA',  NOTE: uses blastn instead of blastp, ignores -g,-D)

f:id:kazumaxneo:20200527232148p:plain

f:id:kazumaxneo:20200527232151p:plain

f:id:kazumaxneo:20200527232157p:plain




引用

GET_HOMOLOGUES, a Versatile Software Package for Scalable and Robust Microbial Pangenome Analysis
Bruno Contreras-Moreira, Pablo Vinuesa

Appl Environ Microbiol. 2013 Dec; 79(24): 7696–7701

 


*1

DFASTでアノテーション付けすると上手く認識したが、PROKKAでアノテーション付けしたGenBnakファイルを使うとエラーになった。

 

#自分用メモ 

#pan-genome解析でコア遺伝子を探索。50genomeとしたら50で全ゲノムに存在。refも別途指定しているなら(refもgenbank-dirに含めないと解析されないので注意)
./get_homologues.pl -d genbank-dir/ -M -t 50 -r ref.gbk -e -n 20
  • -r    reference proteome .faa/.gbk file (by default takes file with least sequences; with BDBH sets first taxa to start adding genes)

 

 

VCFのSVコールをロングリードでジェノタイピングする SVJedi

 

構造変異(SV)の研究は急速に拡大している。その結果、第三世代シークエンシング技術のおかげで、特にヒトゲノムにおいて発見されたSVの数が増加している。同時に、臨床診断のようないくつかのアプリケーションでは、新たにシーケンシングされた個体を、よく定義された特徴的なSVでジェノタイピングすることが重要になる。これまでに、ショートリードデータを対象としたSVジェノタイピングツールはいくつか開発されているが、Pacific BiosciencesやOxford Nanopore Technologiesのように、新たにロングリードシーケンスされたサンプルに既知のSVが存在するかどうかを評価するための専用のツールはなかった。
 ロングリードシーケンシングデータから既知のSVをジェノタイピングする新しい方法を提示する。この方法は、各構造変異の2つの対立遺伝子を表す代表的な対立遺伝子配列のセットを生成することに基づいている。ロングリードは、これらの対立遺伝子配列にアラインメントされる。アライメントを解析し、情報量の多いものだけを残すようにフィルタリングして、各SV対立遺伝子の存在と対立遺伝子頻度を定量化して推定する。長いリードを持つSVをジェノタイプ化するためのSVJedi法の実装を提供する。この手法は、シミュレーションと実データセットの両方に適用され、高いジェノタイピング精度を達成した。本研究では、SVJediが他の既存のロングリードジェノタイピングツールよりも優れた性能を有していることを示し、また、他のアプローチ、すなわちSV発見やショートリードSVジェノタイピングアプローチと比較して、SVジェノタイピングが大幅に改善されていることを示した。

 

インストール

依存

  • Python3
  • Minimap2
  • NumPy
  • Biopython

Github

#bioconda(link)
conda install -c bioconda -y svjedi

svjedi.py -h

$ svjedi.py -h

usage: svjedi.py [-h] [--version] -v <vcffile> [-r <refgenomefile>]

                 [-a <refallelefile>] [-i [<readfile> [<readfile> ...]]]

                 [-p <paffile>] [-t <nb_threads>] [-o <output>]

                 [-dover <dist_overlap>] [-dend <dist_end>] [-ms <minNbAln>]

                 [-d [<seq data type>]]

 

Structural variations genotyping using long reads

 

optional arguments:

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

  --version             show program's version number and exit

  -v <vcffile>, --vcf <vcffile>

                        vcf format

  -r <refgenomefile>, --ref <refgenomefile>

                        fasta format

  -a <refallelefile>, --allele <refallelefile>

                        fasta format

  -i [<readfile> [<readfile> ...]], --input [<readfile> [<readfile> ...]]

                        reads

  -p <paffile>, --paf <paffile>

                        PAF format

  -t <nb_threads>, --threads <nb_threads>

                        Number of threads

  -o <output>, --output <output>

                        genotype output file

  -dover <dist_overlap>

                        breakpoint distance overlap

  -dend <dist_end>      soft clipping length allowed for semi global

                        alingments

  -ms <minNbAln>, --minsupport <minNbAln>

                        Minimum number of alignments to genotype a SV

                        (default: 3>=)

  -d [<seq data type>], --data [<seq data type>]

 

テストラン

git clone https://github.com/llecompte/SVJedi.git
cd SVJedi/HG002_son/
svjedi.py -v test.vcf -r genome.fasta -i simulated-reads.fastq.gz

 入力(test.vcf)

f:id:kazumaxneo:20200526225511p:plain

出力(genotype_results.vcf)

f:id:kazumaxneo:20200526225507p:plain

 

引用

SVJedi: Genotyping structural variations with long reads
Lolita Lecompte, Pierre Peterlongo, Dominique Lavenier, Claire Lemaitre
Bioinformatics, Published: 21 May 2020

ロングリードのマルチプルシーケンスアラインメントを行う abPOA (POAのSIMD実装拡張)

 

 マルチプルシーケンスアラインメント(MSA)問題を解決するために、Leeら(2002)によって最初に導入されたのがPartial order alignment(POA)である。POAでは、MSAをdirected acyclic graph(DAG)(有向無閉路グラフ)として表現し、動的計画法wiki)(DP)を用いて配列をDAGに反復的にアラインメントする。その後、アラインメントグラフに最も重いバンドル化アルゴリズムを適用することで、複数のコンセンサス配列が生成される(Lee, 2003)。最近、Pacific Biosciences(PacBio)およびOxford Nanopore Technologies(ONT)プラットフォーム上でのロングリードシーケンシング技術の出現と人気の高まりに伴い、POAの再評価と関心が高まっており、このアルゴリズムは現在、エラーが発生しやすいロングリードのエラー訂正とアセンブリに広く使用されている(Loman et al、2015; Vaser et al、2017; Volden et al、2018; Gao et al、2019; Ruan and Li、2019)。
POAは古典的なMSAアルゴリズム(Lassmann and Sonnham-mer, 2002)よりもはるかに高速であるが、最先端のロングリードシーケンシング・プラットフォームによって生成される大規模なデータセットが大きな課題となっている。この課題に対処するため、オリジナルのPOAアルゴリズムを高速化するSIMD(Single Instruction Multiple Data)実装が使用された(Vaser et al.) POAのこのSIMDバージョンであるSPOAは、複数の要素を並列に処理する現代のプロセッサにおけるより広いSIMDレジスタを利用する。SIMDベクタは、DP行列の各行の複数の連続したセルのスコアを格納するために使用され、SIMD命令を使用して処理され、各ベクタに格納されたすべてのスコアのために並列更新される。SIMD並列化に加えて、もう一つのアクセラレーション戦略である「banded DP」もまた、配列対配列のアライメントツールにおいて広く用いられている(Chao et al. 具体的には、DP行列の各行または列において、特定の「バンド」内のセルのみを塗りつぶす必要がある。最近、この戦略のグラフバージョンがGraphAligner(Rautiainen and Marschall, 2019)で探求されたが、この場合、DPマトリックスの各行のバンドは、その行の最小スコアセルに基づいて動的に定義される。他のセルは、スコアがその最小スコアから一定の距離内にある場合にのみ、バンド内にあるとみなされる。GraphAlignerは他のグラフアライメントツールよりも高速だが、その戦略は、エラー率の高いロングリードのMSAやコンセンサスコールには適していないかもしれない。これは、GraphAlignerが一般的なスコアリング関数ではなく編集距離をアライメント指標として使用しているためであり、PacBioやONTのシークエンシングデータによく見られる挿入・欠失エラーの多い領域では、誤ったアライメントが発生する可能性がある。
 本研究では、SIMDを実装した適応型バンドDPを実行するPOAの拡張版であるabPOAを開発した。abPOAはスタンドアロンのMSAやコンセンサスコールツールとして機能するほか、任意のロングリードシーケンシングエラー訂正やアセンブリワークフローに簡単に統合することができる。

 

 

インストール

ubuntu18.04LTSとmacos10.14でテストした。

ビルド依存

  • gcc (>=6.4.0) and zlib

Github

#bioconda
conda install -c bioconda abpoa -y

#from source
git clone https://github.com/yangao07/abPOA.git
cd abPOA; make

#dot出力も行うならgraphvizも必要。
sudo apt-get install graphviz

abPOA

$ abPOA

 

abpoa: adaptive banded Partial Order Alignment

 

Usage: abpoa [options] <in.fa/fq> > cons.fa/msa.out

 

Options:

  Alignment:

    -m --aln-mode INT       alignment mode [0]

                              0: global, 1: local, 2: extension

    -M --match    INT       match score [2]

    -X --mismatch INT       mismatch penalty [4]

    -O --gap-open INT(,INT) gap opening penalty (O1,O2) [4,24]

    -E --gap-ext  INT(,INT) gap extension penalty (E1,E2) [2,1]

                            abPOA provides three gap penalty modes, cost of a g-long gap:

                            - convex (default): min{O1+g*E1, O2+g*E2}

                            - affine (set O2 as 0): O1+g*E1

                            - linear (set O1 as 0): g*E1

  Adaptive banded DP:

    -b --extra-b  INT       first adaptive banding parameter [10]

                            set b as < 0 to disable adaptive banded DP

    -f --extra-f  FLOAT     second adaptive banding parameter [0.01]

                            the number of extra bases added on both sites of the band is

                            b+f*L, where L is the length of the aligned sequence

  Input/Output:

    -l --in-list            input file is a list of sequence file names [False]

                            each line is one sequence file containing a set of sequences

                            which will be aligned by abPOA to generate a consensus sequence

    -o --output   FILE      ouput to FILE [stdout]

    -r --result   INT       output result mode [0]

                            - 0: consensus (FASTA format)

                            - 1: MSA (PIR format)

                            - 2: both 0 & 1

    -g --out-pog  FILE      dump final alignment graph to FILE (.pdf/.png) [Null]

 

    -h --help               print this help usage information

    -v --version            show version number

  Parameters under development:

    -a --cons-alg INT       algorithm to use for consensus calling [0]

                            - 0: heaviest bundling

                            - 1: heaviest in column

    -d --diploid            input data is diploid [False]

                            -a/--cons-alg will be set as 1 when input is diploid

                            and at most two consensus sequences will be generated

    -q --min-freq FLOAT     min frequency of each consensus for diploid input [0.30]

 

 

テストラン

git clone https://github.com/yangao07/abPOA.git
cd abPOA/test_data/
abpoa seq.fa > cons.fa

入力の配列

$ cat seq.fa 

>1

CGTCAATCTATCGAAGCATACGCGGGCAGAGCCGAAGACCTCGGCAATCCA

>2

CCACGTCAATCTATCGAAGCATACGCGGCAGCCGAACTCGACCTCGGCAATCAC

>3

CGTCAATCTATCGAAGCATACGCGGCAGAGCCCGGAAGACCTCGGCAATCAC

>4

CGTCAATGCTAGTCGAAGCAGCTGCGGCAGAGCCGAAGACCTCGGCAATCAC

>5

CGTCAATCTATCGAAGCATTCTACGCGGCAGAGCCGACCTCGGCAATCAC

>6

CGTCAATCTAGAAGCATACGCGGCAAGAGCCGAAGACCTCGGCCAATCAC

>7

CGTCAATCTATCGGTAAAGCATACGCTCTGTAGCCGAAGACCTCGGCAATCAC

>8

CGTCAATCTATCTTCAAGCATACGCGGCAGAGCCGAAGACCTCGGCAATC

>9

CGTCAATGGATCGAGTACGCGGCAGAGCCGAAGACCTCGGCAATCAC

>10

CGTCAATCTAATCGAAGCATACGCGGCAGAGCCGTCTACCTCGGCAATCACGT

上記のコマンドだとコンセンサス配列が出力される。

$ cat cons.fa 

>Consensus_sequence

CCACGTCAATCTATCGAAGCATACGCGGCAGAGCCGAAGACCTCGGCAATCAC

 

PIR fomatのマルチプルシークエンスアラインメント(多重整列)を出力。

abpoa seq.fa -r1 > cons.out
  • -r    output result mode [0]
     0: consensus (FASTA format)
     1: MSA (PIR format)
     2: both 0 & 1

出力

 cat cons.out 

>Multiple_sequence_alignment

---CGTCAAT-CTA-TC--G---AAGCA---TACG--CGGGC-AGAGCC--GAA---GACCTCGG-CAATCCA--

CCACGTCAAT-CTA-TC--G---AAGCA---TACG--C-GGC---AGCC--GAACTCGACCTCGG-CAATCAC--

---CGTCAAT-CTA-TC--G---AAGCA---TACG--C-GGC-AGAGCCCGGAA---GACCTCGG-CAATCAC--

---CGTCAATGCTAGTC--G---AAGCA---GCTG--C-GGC-AGAGCC--GAA---GACCTCGG-CAATCAC--

---CGTCAAT-CTA-TC--G---AAGCATTCTACG--C-GGC-AGAGCC--------GACCTCGG-CAATCAC--

---CGTCAAT-CTA-----G---AAGCA---TACG--C-GGCAAGAGCC--GAA---GACCTCGGCCAATCAC--

---CGTCAAT-CTA-TC--GGTAAAGCA---TACGCTC-TGT---AGCC--GAA---GACCTCGG-CAATCAC--

---CGTCAAT-CTA-TCTTC---AAGCA---TACG--C-GGC-AGAGCC--GAA---GACCTCGG-CAATC----

---CGTCAAT-GGA-TC--G----AG-----TACG--C-GGC-AGAGCC--GAA---GACCTCGG-CAATCAC--

---CGTCAAT-CTAATC--G---AAGCA---TACG--C-GGC-AGAGCC--G---TCTACCTCGG-CAATCACGT

 

abPOAはCのAPIsとしても機能する。詳細はGIthubとペーパーを確認して下さい。

 

引用

abPOA: an SIMD-based C library for fast partial order alignment using adaptive band
Yan Gao, Yongzhuang Liu, Yanmei Ma, Bo Liu, Yadong Wang, Yi Xing

bioRxiv, posted May 10, 2020

CRISPR/Cas9編集後のアンプリコンシークエンシングからindelのレポートを生成する CRISPR-DAV

 

 CRISPR/Cas9システムの簡便さと精度の高さは、遺伝子編集の新時代をもたらした。CRISPRを介在させたゲノム編集を用いた目的のクローンのスクリーニングは、その多重化により次世代シークエンシング(NGS)によって可能になった。ここでは、CRISPR NGSデータをハイスループットで解析するためのCRISPR-DAV(CRISPR Data Analysis and Visualization)パイプラインを紹介する。パイプラインでは、Burrows-Wheeler AlignerとAssembly Based ReAlignmentを用いて、small indelとlarge indelの検出を行い、結果は包括的なチャートとインタラクティブなアライメントビューのセットで表示される。CRISPR-DAVはGitHubとDocker Hubのリポジトリhttps://github.com/pinetree1/crispr-dav.githttps://hub.docker.com/r/pinetree1/crispr-dav/ で利用できる。

 CRISPR-DAVパイプラインは、アンプリコンベースのNGSから生成されたFASTQリードを処理する。低品質で非常に短いリードは、まずPRINSEQを用いて除去される(Schmieder et al、2011)。次に、BWAとABRAを使用して、高品質のリードをリファレンスゲノムまたはアンプリコン配列にアラインメントし、indelを検出する。著者らの経験では、ABRAはBWAによって検出されなかった大きなindelを検出した。アンプリコン内の全てのヌクレオチド位置でのindelの頻度を計算する。しかし、CRISPRからのindelは通常、すべてが同じヌクレオチド位置から始まるわけではない。そこで、CRISPRの効率を評価するために、リードレベルでの簡易的な測定値、%indelリードを定義する。まず、sgRNA(シングルガイドRNA)の標的配列にまたがるリードを総リードとしてカウントする。第二に、著者らの観察に基づいて、CRISPRによって誘発されたすべてのindelは、sgRNAの配列と重なっている。したがって、CRISPRによるindelリードとみなされるためには、このターゲット領域に少なくとも1つの塩基が挿入または欠失している必要がある。次に、リードの合計数に応じたindelリードの割合を計算する。HDRの効率を評価するために、リード中のHDRオリゴの所望の塩基変化を調べ、リードを4つのカテゴリに特徴付ける。(i)完全なHDR:すべての意図した塩基変化が起こり、それらの間にindelが存在しない、(ii)編集されたHDR:少なくとも1つの意図した塩基変化が起こるが、indelが存在する、指示された修復が起こった後にCRISPRによる再編集が原因である可能性が高い、(iii)部分的なHDR:一部の意図した塩基変化が起こるが、すべてではないが、indelが存在しない、(iv)非HDR:意図した塩基変化のいずれも示さないリード。データの理解を深めるためには、可視化が重要である。パイプラインは、以下を示すチャートを含むHTMLレポートを生成する:様々な段階でのリードカウント、アンプリコン内のリードの深さとindel頻度、indelリードのカウントとパーセンテージ、アレル、SNP、HDRの頻度、Canvas Xpress(http://canvasxpress.org)で有効化されたアライメントビュー、となる。

 

 

インストール

準備されているdockerイメージを使ってテストした。

Github

#dockerhub(link)
docker pull pinetree1/crispr-dav:latest

> /crispr-dav/crispr.pl

# /crispr-dav/crispr.pl

CRISPR data analysis and visualization pipeline.

 

Usage: ../../crispr.pl [options] 

 

    --conf <str> Configuration file. Required. See template /crispr-dav/conf.txt

        It has information about genome locations, tools, and parameters.  

 

    Specify a reference using --genome or --amp_fasta, but not both. 

    Use --genome for standard genome, such as hg19. Need to have paths of fasta file, 

    bwa index, and refGene coordinate file in the configuration file. To download the  

    coordinate file, go to UCSC Genome Browser, in TableBrowser, select group:Genes 

    and Gene Predictions, track:RefSeq Genes, table:refGene, region:genome,

    output format:all fields from selected table. The downloaded tab-delimited file 

    should have these columns:

    bin,name,chrom,strand,txStart,txEnd,cdsStart,cdsEnd,exonStarts,exonEnds,... 

    Use --amp_fasta when using a custom amplicon sequenece as reference. 

 

    --genome <str> Genome version (e.g. hg19) as specified in configuration file.

 

    --amp_fasta <str> Amplicon reference fasta file containing a single sequence. 

    --codon_start <int> Translation starting position in the amplicon reference sequence. 

        If the first codon starts at the first base, then the position is 1. No translation 

        will be performed if the option is omitted. No intron should be present in the 

        amplicon reference sequence if translation is needed. 

 

    --region <str> Required when --genome option is used. This is a bed file for amplicon region.

        The tab-separated fields are chr, start, end, genesym, refseqid, strand(+/-). 

        No header. All fields are required.

        The start and end are 0-based; start is inclusive and end is exclusive.

        Genesym is gene symbol. Refseqid is used to identify transcript coordinates in 

        UCSC refGene coordinate file. If refseqid is '-', no alignment view will be created. 

        Only one row is allowed this file. If an experiment has two amplicons, run the 

        pipeline separately for each amplicon. 

 

    --crispr <str> Required. A bed file containing one or more CRISPR sgRNA sites.

        Tab-delimited file. No header. Information for each site:

        The fields are: chr, start, end, CRISPR_name, sgRNA_sequence, strand, and 

        HDR mutations.  All fields except HDR mutations are required. The start and end 

        are 0-based; start is inclusive and end is exclusive. CRISPR names and sequences 

        must be unique. 

        HDR format: <Pos1><NewBase1>,<Pos2><NewBase2>,... The bases are desired new bases 

        on positive strand,e.g.101900208C,101900229G,101900232C,101900235A. No space. The 

        positions are 1-based and inclusive. 

 

    --fastqmap <str> Required. A tab-delimited file containing 2 or 3 columns. No header.

        The fields are sample name, read1 fastq file, and optionally read2 fastq file.

        Fastq files must be gizpped and and file names end with .gz.

 

    --sitemap <str> Required. A tab-delimited file that associates sample name with CRISPR 

        sites. No header. Each line starts with sample name, followed by one or more sgRNA

        guide sequences. This file controls what samples to be analyzed. 

 

    --merge Y or N. Default: Y. Merge paired-end reads before filtering and alignment.

    --sge Submit jobs to SGE queue. The system must already have been configured for SGE.

    --outdir <str> Output directory. Default: current directory.

    --help  Print this help message.

    --verbose Print some commands and information for debugging.

 

 

テストラン1

実行

docker run --rm -it -v $PWD:/Users/xyz/temp pinetree1/crispr-dav 
cd /crispr-dav/Examples/example1/
sh run.sh

以下のコマンドを実行している。

../..//crispr.pl --conf conf.txt --region amplicon.bed --crispr site.bed \
--sitemap sample.site --fastqmap fastq.list --genome genomex

conf.txtはコンフィグファイル。先頭のgenomeのパスだけ書き換えれば他のデータにも使える。

> cat conf.txt

f:id:kazumaxneo:20200524222547p:plain

(以下省略)

amplicon.bedは増幅領域を指定したbedファイル。 site.bedはCRISPR sgRNA sitesを指定したbedファイル。sample.siteはサンプル名とCRISPR sitesを指定したタブ区切りファイル。fastq.listはfastqのパスを指定したタブ区切りファイル。

f:id:kazumaxneo:20200524222930p:plain


ラン後、結果をホストと共有しているディレクトリに移してから終了。

mv deliverables /Users/xyz/temp
exit

index.html

f:id:kazumaxneo:20200524205917p:plain

f:id:kazumaxneo:20200524205918p:plain

GENEX_CR1_cx0.html

f:id:kazumaxneo:20200525000429p:plain

f:id:kazumaxneo:20200525000542p:plain

 

図はインタラクティブに操作できる。

f:id:kazumaxneo:20200525000511p:plain

 

テストラン2

準備。テスト1と違ってサンプルシートを用意し、ヘルパースクリプトを使ってランに必要なファイルを調整する。

docker run --rm -it -v $PWD:/Users/xyz/temp pinetree1/crispr-dav 
cd /crispr-dav/Examples/example1_a/
../../prepare_run.pl samplesheet.txt

ラン

cd amp_GENEX_chr4_35769_36280/
sh run.sh

 

引用
CRISPR-DAV: CRISPR NGS Data Analysis and Visualization Pipeline

Xuning Wang, Charles Tilford, Isaac Neuhaus, Gabe Mintier, Qi Guo , John N Feder, Stefan Kirov

Bioinformatics. 2017 Dec 1;33(23):3811-3812