macでインフォマティクス

macでインフォマティクス

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

再現性のあるメタゲノム解析を行うためのモジュール設計された自動パイプライン Sunbeam

2019 6/26 誤字修正 

 

 メタゲノミックショットガンシークエンシングは、関心のある微生物混合群からDNAを抽出し、無作為に抽出されたDNAをディープシーケンシングする。これは、特定の標的遺伝子領域が増幅およびシーケンシングされるマーカー遺伝子シーケンシング(例えば、細菌の16S rRNA遺伝子)とは対照的である。メタゲノムシーケンシングは、特にウイルスおよびバクテリオファージコミュニティの研究において、微生物学において重要な洞察を可能にしました[ref.1–9]。ただし、継続的な課題は、結果として得られる大規模データセットを標準的で信頼性の高い方法で分析および解釈することである[ref.16–23]。

 微生物のメタゲノムを調べる一般的な方法は、イルミナのシークエンシング技術を使用して、目的のサンプルから単離した断片化DNAから多数の短い(100〜250塩基対)リードを取得することである。シーケンスの取得後、シーケンスを使用して基礎となる生物学的洞察を得る前に、いくつかの後処理ステップを実行する必要がある[ref.21、24]。

 研究者はそれぞれの後処理ステップを実行するための多くのツールを自由に使えるようにしているが、その結果として、得られるアウトプットや下流の分析を変える可能性がある各ツールそれぞれに様々なパラメータが存在する。分析間でパラメータ、ツール、またはリファエンスデータベースのバージョンを変えると、異なるメタゲノムシーケンシング実験の結果を比較するのが難しくなる[ref.25]。逆に、研究間で一貫したワークフローを採​​用することで、実験は比較可能であり、下流分析は再現可能であることが保証される[ref.21]。使用されるソフトウェア、データベース、およびパラメータの文書化は、このプラクティスの重要な要素である。そうでなければ、一貫した再現可能なワークフローの利点が失われる。

 メタゲノムのポストプロセッシングワークフローは、その有用性と柔軟性を最大にするために以下の特性を持つべきである。それは共有コンピューティングシステムとクラウドにデプロイ可能であるべきであり、ソフトウェアパラメータとリファレンスデータベースの簡単な設定を許容する。中断の後に再開する能力、それは不必要なステップがスキップされるか無視されることができるようにモジュール式であるべきである、そしてそれは新しい手順がユーザによって加えられることを可能にするべきである。機関プラットフォームとクラウドプラットフォームの両方にワークフローを展開する機能により、さまざまなラボでさまざまなコンピューティング設定でワークフローを繰り返すことができ、研究者が機関とクラウドのコンピューティングリソースを柔軟に選択できるようになる。同様に、設定ファイルを使用して実行中のパラメータを記録する機能により、実験固有のソフトウェアパラメータを使用することができ、今後の参考資料として使用できる。

 いくつかの機能が効率的なデータ分析に貢献する。ワークフローのエラーや中断が最初からやり直す必要がない場合は、有益である。シーケンス実験では大量のデータが生成されるため、データ処理のステップを繰り返す必要があり、時間と費用がかかる可能性がある。さらに、ワークフローのすべてのステップがすべての実験に必要なわけではなく、実験によってはカスタム処理が必要な場合もある。実験を適切に処理するために、ワークフローは不要なステップをスキップし、必要に応じて後で実行する簡単な方法を提供するべきである。フレームワークを広く有用にするために、ユーザーは必要に応じてワークフローに新しいステップを簡単に追加して他のユーザーと共有できなければならない。これらの目標の多くを達成するパイプラインがいくつか開発されているが[ref.26-29]、メタゲノムデータセットの処理における柔軟性の向上と分析の長期的な再現性に対する我々(本著者ら)のニーズを満たしていなかった。

 ここでは、メタゲノムシーケンス実験から一貫した一連の後処理済みファイルを生成する、簡単に配置および構成可能なパイプラインであるSunbeamを紹介する。 Sunbeamは自己完結型で、管理者権限なしでGNU / Linuxシステムにインストールできる。 Snakemakeワークフロー言語での実装により、エラー処理、タスク再開、および並列コンピューティング機能を備えている[ref.30]。ほぼすべての手順は事前に指定された妥当なデフォルト値を使用して構成可能であるため、広範なパラメータ調整を行わなくても迅速に展開できる。 Sunbeamは、ローカルデータを直接使用することも、NCBIのシーケンスリードアーカイブ(SRA)からの外部データを使用することもできる[ref.31]。 Sunbeamは、新しい手順をユーザーが追加できるようにする単純なメカニズムを使用して拡張可能である。

 さらに、Sunbeamは、低品質またはホスト由来の配列を多く含む困難なサンプルのデータを処理できるカスタムソフトウェアを備えている。これらには、任意の数の宿主ゲノムまたは汚染ゲノム配列に対するカスタム調整された宿主由来のリード除去ステップ、およびKomplexity、問題のある低複雑性リードを下流の分析の前に迅速かつ正確に除去する新しい配列複雑性分析プログラムが含まれる。マイクロサテライトDNA配列は、ヒトゲノムのかなりの部分を占めており、個体間で非常に多様であり[ref.32-34]、単一のリファレンスゲノムに対するアラインメントによってそれらを除去することの難しさを増している。 Komplexityを開発したのは、ヌクレオチド配列の複雑さを分析するための既存のツール[ref.35-37]が、スピード、怪しいヒットの除去、およびfastqファイルのネイティブ処理の点で、本著者らのニーズを満たしていなかったためである。ここでは宿主関連の low-microbial biomass body sites [ref.19、38、39]、virome[ref.40]、およびマイクロバイオームの縦断サンプリング[ref.41、42]の公表され進行中の研究においてSunbeamを使用した。

 SunbeamはPythonBash、およびRustで実装されている。 GPLv3でライセンスされている。https://github.com/sunbeam-labs/sunbeamで自由に利用可能であり、ドキュメントはhttp://sunbeam.readthedocs.ioで利用可能である。

 Sunbeamは、「可用性と要件」に記載されているハードウェアとソフトウェアの初期要件を満たすGNU / Linuxディストリビューションにインストールできる。インストールは、リポジトリからソフトウェアをダウンロードして「install.sh」を実行することによって実行される。Sunbeamは、インストールまたは実行に管理者権限を必要としない。 Sunbeamの基本的な分析ワークフローをDebian 9、 CentOS 6と7、Ubuntu 14.04、16.04、18.04、および18.10、 Red Hat Enterprise 6および7、そしてSUSE Enterprise 12にインストールして実行できることを確認した。 

 SunbeamはCondaパッケージ管理システム[ref.43]を利用して、「可用性と要件」のセクションに記載されているもの以外にも追加のソフトウェア依存関係をインストールする。 Condaは、依存関係の自動解決を提供し、管理者以外のユーザーにソフトウェアのインストールを容易にし、追加のソフトウェアをパッケージ化してサードパーティのソフトウェアチャネルを管理するための標準化されたシステムを使用する。 Condaは、パイプラインの外側にある既存のソフトウェアとの競合を避けるために、追加のソフトウェア依存関係解決のための分離された環境も提供する。必要に応じて、CondaはSunbeamインストールスクリプトによってインストールされる。 Sunbeamが使用する分離ソフトウェア環境もインストールスクリプトによって作成される。

  Sunbeamパイプラインへの入力は、ローカルのファイルまたはSRAを通じて利用可能なサンプルのいずれかで、raw Illuminaシーケンスリードからなる。デフォルトでは、Sunbeamはリードに対して次の順序で予備操作を実行する。

クオリティ管理:オプションで、grabseq [ref.44]とsra-tools [ref.31]を使ってSRAからリードをダウンロードする。アダプター配列を除去し、そしてCutadapt [ref.45]およびTrimmomatic [ref.46]ソフトウェアを用いて塩基をクオリティフィルタリングする。クオリティフィルタリングを通過したリードペアは保持される。クオリティはFastQCを使用して評価され[ref.47]、別々のレポートにまとめられる。

低複雑度マスキング:各リードにおける配列複雑度は、後述する新規複雑度スコアリングアルゴリズムであるKomplexityを使用して評価される。ユーザーがカスタマイズ可能なシーケンス複雑度のしきい値を下回るリードは削除される。削除されたリード数のログは後で検査するために書き込まれる。

ホストリード除去:リードはbwaを使用してユーザー指定のホストまたは汚染物質シーケンスのセットに対してマッピングされる。特定の同一性および長さのしきい値内でこれらのシーケンスのいずれかにマップされるリードは削除される。削除されたリード数は後で検査するために記録される。

初期クオリティ管理プロセスの後、複数のオプションの下流工程を並行して実施することができる。分類ステップでは、汚染配列が除去されクオリティ管理されたリードがKraken [ref.49]を使用して分類的に分類され、タブ区切り形式とBIOM [ref.50]形式の両方で要約される。アセンブリステップでは、MEGAHIT [ref.51]を使用して、各サンプルからのリードを連続アセンブリする。事前に指定された長さを超えるコンティグには、あのテーションが付けられる。オープンリーディングフレーム(ORF)は、Prodigal [ref.52]を用いて予測される。次に、コンティグ全体(および関連するORF)を、コンティグ全体および推定ORFの両方を使用して、任意の数のユーザー指定のヌクレオチドまたはタンパク質BLAST [ref.53]データベースに対して検索する。結果は各サンプルのレポートにまとめられる。最後に、マッピングステップでは、bwa [ref.48]を使用してクオリティ管理されたリードを任意の数のユーザー指定のリファレンスゲノムまたは遺伝子セットにマッピングし、結果のBAMファイルをSAMtools [ref.54]を使用してソートおよびインデックス付けする。

Sunbeamは、出力ファイルを概念的に異なるフォルダにグループ化し、論理的な出力フォルダ構造を提供するよう設計されている。 Sunbeamからの標準出力には、クオリティ管理プロセスの各ステップからのfastqファイル、各リードの分類系統アサイン、各サンプルからアセンブリされたコンティグ、遺伝子予測、およびすべてのクオリティ管理されたリードのアラインメントファイルが含まれる。ほとんどのルールは、後の検査と要約のためにその操作のログを作成する。ユーザーは特定の出力を個別にまたはグループとして要求でき、パイプラインは必要なファイルを作成するのに必要なステップのみを実行する。これにより、ユーザーはモジュール式にパイプラインの任意の部分をスキップまたは再実行できる。

 

表1、他のツールとの機能比較

f:id:kazumaxneo:20190622233206p:plain

論文より転載

  

以下の作業が自動化されている (Githubより)

  • Quality control, including adaptor trimming, host read removal, and quality filtering;
  • Taxonomic assignment of reads to databases using Kraken;
  • Assembly of reads into contigs using Megahit;
  • Contig annotation using BLAST[n/p/x];
  • Mapping of reads to target genomes; and
  • ORF prediction using Prodigal.

 

Documentation

https://sunbeam.readthedocs.io/en/latest/

 

インストール

ubuntu16.0.4のminiconda3-4.3.30環境でテストした。

依存

  • python3.6 ~ 3.7

Gihutb

git clone -b stable https://github.com/sunbeam-labs/sunbeam sunbeam-stable
cd sunbeam-stable

#condaがない環境で実行するとcondaも含めて導入される
./install.sh
#condaをここで入れたなら
echo "export PATH=$PATH:/root/miniconda3/bin" >> ~/.bashrc && source ~/.bashrc
source activate sunbeam

テスト

#テストラン時に出た文字化けエラーはC.UTF-8をexportして応急処置

export LC_ALL=C.UTF-8
export LANG=C.UTF-8

> tests/run_tests.bash -e sunbeam

f:id:kazumaxneo:20190622180053j:plain

sunbeam

# sunbeam

usage: sunbeam [-h/--help,-v/--version] <subcommand>

 

subcommands:

  init         Create a new config file for a project using local or SRA data.

  run          Execute the pipeline.

  config       Modify or update config files.

  list_samples Make a list of samples from a directory.

 

optional arguments:

  -v, --version  show program's version number and exit

 

For more help, see the docs at http://sunbeam.readthedocs.io.

Unrecognized command.

> sunbeam init -h

# sunbeam init -h

usage: init [-h] [-f] [--output FILE] [--defaults FILE] [--template FILE]

            [--data_fp PATH] [--format STR] [--single_end]

            [--data_acc ACC [ACC ...]]

            project_fp

 

Initialize a new Sunbeam project in a given directory, creating a new config

file and (optionally) a sample list. The sample list source can be either a

folder of input files or a list of SRA accession numbers.

 

positional arguments:

  project_fp            project directory (will be created if it does not

                        exist)

 

optional arguments:

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

  -f, --force           overwrite files if they already exist

  --output FILE         name of config file (sunbeam_config.yml)

 

config file options:

  --defaults FILE       set of default values to use to populate config file

  --template FILE       custom config file template, in YAML format

 

sample list options:

  Options to automatically generate a sample list. --data_fp (for reading

  filenames from a specified folder) and --data_acc (for fetching SRA

  accession numbers to be downloaded during a run) are mutually exclusive.

  If --data_acc is used, all SRA Run (sample) entries corresponding to the

  given accession numbers (e.g. PRJNA###, SAMN###, SRR###) will be used to

  create the sample list. Note in this case project_fp should either be

  given before --data_acc or separated by a '--' to distinguish it from an

  accession number.

 

  --data_fp PATH        path to folder containing fastq.gz files

  --format STR          filename format for --data_fp (default: guessed)

  --single_end          fastq files are in single-end, not paired-end, format

                        for --data_fp

  --data_acc ACC [ACC ...]

                        list of SRA-compatible accession numbers

> sunbeam run -h

# sunbeam run -h

usage: sunbeam run [options] -- <snakemake options>

 

Executes the Sunbeam pipeline by calling Snakemake.

 

optional arguments:

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

  -s SUNBEAM_DIR, --sunbeam_dir SUNBEAM_DIR

                        Path to Sunbeam installation

 

You can pass further arguments to Snakemake after --, e.g:

    $ sunbeam run -- --cores 12

See http://snakemake.readthedocs.io for more information.

sunbeam config -h

# sunbeam config -h

usage: sunbeam config [-h] {update,modify} ...

 

positional arguments:

  {update,modify}

 

optional arguments:

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

> sunbeam list_samples -h

# sunbeam list_samples -h

usage: sunbeam list_samples [-h] [-s] [-f STR] data_fp

 

positional arguments:

  data_fp               Path to folder containing reads

 

optional arguments:

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

  -s, --single_end      Samples are single-end (not paired-end)

  -f STR, --format STR  Filename format (e.g. {sample}_R{rp}.fastq.gz).

                        (default: guessed)

 

最小限の手間でパイプラインにツールを追加することも可能になっている(未テスト)。

#sunbeam内にて
git clone https://github.com/sunbeam-labs/sbx_kaiju extensions/sbx_kaiju
cd extensions/sbx_kaiju/
conda install -c bioconda --file requirements.txt

#databaseをダウンロード後、ymlファイルにパスを追加。(kaiju github)
mkdir kaijudb
cd kaijudb
kaiju-makedb -s (database name)

#それからrun前にymlファイルに追加する
cat extensions/sbx_kaiju/config.yml >> my_config.yml

 

 

データベースの準備

1、kraken database(version.1.0)

ここではpre-build (2017)のMiniKraken DB_8GBデータベースを使う。

mkdir /data/kraken
cd /data/kraken/
wget https://ccb.jhu.edu/software/kraken/dl/minikraken_20171019_8GB.tgz
tar -zxvf minikraken_20171019_8GB.tgz && rm minikraken_20171019_8GB.tgz

ymlファイルのclassify:kraken_db_fp: にパスを記載する。上の通りなら

kraken_db_fp: '/data/kraken/minikraken_20171019_8GB'

 

2、filtering of contaminated sequences

コンタミネーション(ホストゲノムやテクニカルな汚染配列) の除去。例えばヒトゲノムの汚染を除く。リファレンスの拡張子はfastaにする必要がある。

mkdir /data/contamination_genome
cd /data/contamination_genome/
wget http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz
unpigz -c hg38.fa.gz -p 8 > hg38.fasta && rm hg38.fa.gz

ymlファイルのhost_fp: にパスを記載する。上の通りなら

host_fp: '/data/contamination_genome'

 

3、blast(n、p、x)

blastしてアノテーションをつけることができる。ここではNCBI NTを使う。抗生物質耐性遺伝子、virulance factor、NCBI nt、など目的に応じてデータベースを構築しておく。

mkdir /data/annotation_database/
cd /data/annotation_database/

#NCBI ntを入れる (容量に注意)
mkdir nt
cd nt
wget "ftp://ftp.ncbi.nlm.nih.gov/blast/db/nt.*.tar.gz"
for a in nt.*.tar.gz; do tar xzf $a; done


#CARDデータベース
mkdir /data/annotation_database/card
#調べてから追記します

#VFDBータベース

ymlファイルの blastdbs以下に記載する。ここではproteinは載せていないのでblastpは実行されない。

f:id:kazumaxneo:20190625194055j:plain

ダウンロードしたntは70以上に分かれているが、sunbeamはnt/ntで認識する。ここではblastデータベースは全て/data/annotation_database/にあり、その下位ディレクトリのnt/にNCBI ntデータベースファイルがあるという設定で進めている。

下の画像の通り、 nt/にはダウンロードし解凍したntのblastデータベースが保存されている。

f:id:kazumaxneo:20190625194446j:plain

 

 

実行方法

1、fastqの準備。

特定のディレクトリに解析対象のfastqをまとめる。 例えばこのようなサンプルを分析する。ファイル名は推測できるような名前になっていないといけない。下のファイル名だと問題なく認識した(sample1, sample2, ...)。

f:id:kazumaxneo:20190623155306j:plain

 

 

2、ランのconfigファイル作成。

ここではfastqを保存しているディレクトリを/sequencing/fastq、プロジェクトファイル名はmy_project1としている。

sunbeam init my_project1 --data_fp /sequencing/fastq

my_project1ディレクトリができ、その中にsamples.csvとsunbeam_config.ymlができる。

 

3、.ymlの編集。他のツールのconfiguration fileに相当するymlファイルができるので、ラン前に編集する。

vim my_project1/sunbeam_config.yml

sunbeam_config.yml

f:id:kazumaxneo:20190626131356j:plain

上から順番に見ていって、パラメータを変更していく。モジュラー設計のパイプラインなのでステップを丸ごと追加・削除することもできる。
 

 

 

4、実行する。-nでdry run

sunbeam run --configfile my_project1/sunbeam_config.yml

出力

f:id:kazumaxneo:20190625194650j:plain

 

 

SRAデータを使ったテストラン

例えばmetagenome 4サンプルを分析する。

mkdir ~/test4strain
sunbeam init /data/test4strain/ --data_acc ERR3277301 ERR3277295 ERR3277320 ERR3277274
#.ymlとsample.csvができる。

#~/testにできた.ymlを指定してラン
sunbeam run --configfile /data/test4strain/sunbeam_config.yml

ymlファイルで出力defaultは"subbeam_output"になっている。

ラン後、subbeam_outputを開く。

f:id:kazumaxneo:20190623130924j:plain

 

dowload、qc、assembly、classify、annotationに分かれている。

f:id:kazumaxneo:20190623130942j:plain

 

download/にはダウンロードされたfastqが収められている。

f:id:kazumaxneo:20190623131008j:plain

 

qc/にはcutadapt、trimmomaticによるアダプター除去/QT、low comlexility reads除去、ホストコンタミネーションリード除去後のそれぞれのfastqが収められている。したがってかなりサイズが大きくなる。

f:id:kazumaxneo:20190623130955j:plain

 

qc/reportsにはQC結果のレポート

f:id:kazumaxneo:20190623131812j:plain

classify/にはdefaultではkrakenによる分析結果が収められている。

f:id:kazumaxneo:20190623131025j:plain

全サンプル結果はBIOM形式とTSV形式でまとめられる。

all_samples.tsv

f:id:kazumaxneo:20190623132003j:plain

biomはqiimeやRのphyloseqなどで可視化できる。

 

assembly/にはmegahitのcontigsやカバレッジファイルが収められている。

f:id:kazumaxneo:20190623131040j:plain

全サンプルの全contigのカバレッジやサイズレポート

f:id:kazumaxneo:20190623131044j:plain



アノテーションにはユーザーが指定したblastデータベースによるアノテーション結果が収められる。

f:id:kazumaxneo:20190623131052j:plain

 

 

snakemakeのスケーラビリティにより、シングルノードサーバーから、クラスターの分散メモリ環境でも特別なチューニングなしにランできるようになっています。パラメータは

https://sunbeam.readthedocs.io/en/latest/usage.html#configuration

のCluster optionsで確認して下さい。

 

 

引用

Sunbeam: an extensible pipeline for analyzing metagenomic sequencing experiments

Clarke EL, Taylor LJ, Zhao C, Connell A, Lee JJ, Fett B, Bushman FD, Bittinger K

Microbiome. 2019 Mar 22;7(1):46.

 

 

 

複数のメタゲノムをその場で分析するための軽量で多機能なメタゲノム分析ツール SqueezeMeta(オフライン使用)

2020 11/19 condaインストール追記

 

  シーケンシング技術の改良によりメタゲノムシーケンシングが一般化し、メタゲノムシーケンシングがマイクロバイオームの構造および機能性を分析するための標準的な手順となった。メタゲノム実験によって生成された膨大な数のショートリード配列に対処するために、多くの新規なバイオインフォマティックツールおよびアプローチが開発されてきた。単に膨大な量のデータ以外に、メタゲノム解析は、結果に直接互換性がないことが多いさまざまなソフトウェアツールといくつかの標準化されていないステップを含む複雑な作業になっている。

 最近、携帯性の高いシーケンサー、特にナノポア技術に基づくシーケンサーの開発(Deamer et al、2016)は、迅速な結果を得る必要性が最も重要であるシナリオ、例えば疾病管理または伝染病の臨床シナリオにおいてin situシークエンシングを容易にした(Quick et al、2015、2016)。メタゲノムシークエンシングもまた、例えば南極氷の海洋調査がその場で行われており(Lim et al、2014; Johnson et al、2017)、サンプリングしてすぐにシークエンスを生成する能力が増大している。これにより、前日の結果に従って、今後のサンプリング実験の情報に基づく計画が可能になる。近い将来、このようなアプリケーションがますます使用されるようになると予測される。したがって、バイオインフォマティック分析は非常に短い時間(数時間)で実行されるべきであり、軽量コンピューティングインフラストラクチャに適しているはずである。

 標準的なメタゲノムパイプラインには、リードのキュレーション、アセンブリ、遺伝子予測、および結果として得られる遺伝子の機能および分類学アノテーションが含まれる。これらの分析の大部分を自動化するために、いくつかのパイプラインが作成されている(Li、2009; Arumugam et al、2010; Glass and Meyer、2011; Abubucker et al、2012; Eren et al、2015; Kim et al、2016) 。ただし、能力とアプローチの点では異なり、最も重要な違いの1つは、アセンブリ手順が必要かどうかになる。いくつかのプラットフォームは、アセンブリをスキップし、その結果、遺伝子予測をスキップし、代わりにrawリードのダイレクトアノテーションに頼っているが、rawリードを扱ういくつかの欠点がある:これは巨大なリファレンスデータベースに対する何百万ものシーケンスの相同性検索に基づいているので、それは通常非常に大きいCPU使用率を必要とする。特に分類学アサインについては、リファレンスデータベースはエラーを最小にするために可能な限り完全でなければならない(Pignatelli et al、2008)。さらに、正確なアサインを生成するにはシーケンス長が短すぎることが多い(Wommack et al、2008; Carr and Borenstein、2014)。

 多くの遺伝子を含むことが多いゲノムのより大きな断片を回収することができるので、アセンブリは賢明な手段である。完全な遺伝子配列およびそのコンテキストを有することは、その機能および分類学アサインをはるかに容易かつより信頼性の高いものにする。アセンブリの欠点は、異なるゲノムパーツを誤ってアセンブリことによるキメラ形成、およびいくつかのリード、特に存在量の少ない種からのリードをアセンブリすることができないことである。アセンブリされていないリードの割合は、いくつかの要因、特にシーケンシングデプスとマイクロバイオームの多様性に左右されるが、通常は少ない(多くの場合20%未満)。最近、最初の段階でアセンブリされなかったリードをリアセンブリし、このステップの性能を向上させるいくつかのツールが開発されている(Hitch and Creevey、2018)。結果のセクションで説明するように、関連するメタゲノムを共にアセンブリすることで、この問題を大幅に軽減することもできる。

 アセンブリはまた、ビニング方法を介して準完全なゲノム回収を容易にするので、賢明である。ゲノムの検索は、微生物と機能を結び付けることを可能にし、よってコミュニティの機能のより正確な生態学的記述に寄与するので、マイクロバイオームの研究における大きな前進につながる。例えば、(特に重要な機能に関与する)マイクロバイオームの重要なメンバーを決定すること、メンバー間の潜在的な相互作用を推測すること(例えば代謝の補完関係を探すこと)、および生態学的効果の理解を進めることが可能である。

 ビニングのための最良の戦略は、関連メタゲノムの共アセンブリである。異なる試料中のコンティグの存在量および組成を比較することによって、どのコンティグが同じ生物に属するかを決定することが可能である:これらのコンティグは類似のオリゴヌクレオチド組成、個々の試料中で類似の存在量、および異なる試料間では共変動パターンを有する。このようにして、マイクロバイオームの機能をより詳細に分析するための出発点として使用することができる、異なる完成レベルを有する数十または数百のゲノムビンを検索することが可能である。

 SqueezeMetaはメタゲノム/メタトランスクリプトミクスのための完全自動のパイプラインで、分析のすべてのステップを網羅している。これは、関連メタゲノムの共アセンブリおよびビニング手順による個々のゲノムの検索を可能にするマルチメタゲノムサポートを含む。

 SqueezeMetaと他のパイプラインの機能の比較を論文の表Table 1.1に示す。現在のほとんどのパイプラインでは、協調およびビニングのサポートは含まれていないが、外部ビニング結果をインポートして関連情報を表示することを許可するものもある。

SqueezeMetaには、既存のパイプラインとは異なる高度な特性がいくつかある。例えば、

1、各メタゲノムにおける個々の遺伝子存在量推定のためのリードマッピングと組み合わせた共アセンブリ手順。
2、個々のメタゲノムのマージを介して無制限の数のメタゲノムの処理を可能にする代替の共アセンブリ手法。
3、ナノポアロングリードのサポート。
4、個々のゲノムを回収するためのビニングおよびビンチェック。
5、コンティグとビンの分類学アノテーションの内部チェック。
6、リファレンスメタゲノムに対してcDNAリードをマッピングするメタトランスクリプトームのサポート。メタゲノムとメタトランスクリプトームの共アセンブリへのマッピングも利用できる。
7、結果を保存するためのMySQLデータベース。Webインターフェースを使用して簡単にエクスポート、共有および検査ができる。

 

我々(著者ら)は、サンプリングしたその場でのメタゲノムシーケンシング実験に利用できるよう、SqueezeMetaをコンピュータリソースの乏しい環境でも実行できるように設計した。すべてのパイプラインのコンポーネントを適切に設定することで、16GBのRAMを搭載したデスクトップコンピューターを使用して、個々のメタゲノムを完全に分析し、さらに関連するメタゲノムを統合することさえできた。本システムは完全に自動化されているため、技術的な知識やバイオインフォマティックな知識を必要としない非常に使いやすいものである。また、インターネット接続有無とは完全に無関係である。

SqueezeMetaはhttps://github.com/jtamames/SqueezeMetaからダウンロードできる。

 

 

f:id:kazumaxneo:20190215205638p:plain

Workflow of the three modes of operation of SqueezeMeta.

論文より転載

 

SqueezeMeta uses a combination of custom scripts and external software packages for the different steps of the analysis: (Githubより)

  1. Assembly
  2. RNA prediction and classification
  3. ORF (CDS) prediction
  4. Homology searching against taxonomic and functional databases
  5. Hmmer searching against Pfam database
  6. Taxonomic assignment of genes
  7. Functional assignment of genes
  8. Taxonomic assignment of contigs, and check for taxonomic disparities
  9. Coverage and abundance estimation for genes and contigs
  10. Estimation of taxa abundances
  11. Estimation of function abundances
  12. Merging of previous results to obtain the ORF table
  13. Binning with MaxBin
  14. Binning with MetaBAT
  15. Binning integration with DAS tool
  16. Taxonomic assignment of bins, and check for taxonomic disparities
  17. Checking of bins with CheckM
  18. Merging of previous results to obtain the bin table
  19. Merging of previous results to obtain the contig table
  20. Prediction of kegg and metacyc patwhays for each bin
  21. Final statistics for the run

 

 

インストール

ubuntu16.04でテストした(ホストOS macos10.14)。

本体 Github

GithubレポジトリのINSTALL-UBUNTU14+(リンク)の手順にしたがってインストールした。

git clone https://github.com/jtamames/SqueezeMeta.git
cd SqueezeMeta/

#以下のインストールスクリプトをdockerのrootで実行した。インストールスクリプトのsudoとスクリプト最後の200GBのデータセットダウンロードを消してから実行した。
bash INSTALL-UBUNTU14+

#blastp、pysamが入っていなかったので手動で導入。usearch(v11)はlinux 32bit版をオーサーのHPからダウンロード

2020 11/19追記(python3.7でテストした)
conda create -n SqueezeMeta -c bioconda -c conda-forge -c fpusan squeezemeta
conda activate SqueezeMeta

> perl scripts/SqueezeMeta.pl

# perl scripts/SqueezeMeta.pl

 

SqueezeMeta v0.4.2, Jan 2019 - (c) J. Tamames, F. Puente-Sánchez CNB-CSIC, Madrid, SPAIN

 

Please cite: Tamames & Puente-Sanchez, bioRxiv 347559 (2018); doi: https://doi.org/10.1101/347559; Tamames & Puente-Sanchez, Frontiers in Microbiology (2019), in press

 

Usage: SqueezeMeta.pl -m <mode> -p <projectname> -s <equivfile> -f <raw fastq dir> <options>

 

Arguments:

 

 Mandatory parameters:

   -m: Mode (sequential, coassembly, merged) (REQUIRED)

   -s|-samples: Samples file (REQUIRED)

   -f|-seq: Fastq read files' directory (REQUIRED)

   -p: Project name (REQUIRED in coassembly and merged modes)

   

 Filtering: 

   --cleaning: Filters with Trimmomatic (Default: No)

   -cleaning_options: Options for Trimmomatic (Default:LEADING:8 TRAILING:8 SLIDINGWINDOW:10:15 MINLEN:30)

   

 Assembly: 

   -a: assembler [megahit,spades,canu] (Default: megahit)

   -assembly_options: Options for required assembler

   -c|-contiglen: Minimum length of contigs (Default: 200)

   

 Mapping: 

   -map: mapping software [bowtie,bwa,minimap2-ont,minimap2-pb,minimap2-sr] (Default: bowtie) 

 

 Annotation:  

   --nocog: Skip COG assignment (Default: no)

   --nokegg: Skip KEGG assignment (Default: no)

   --nopfam: Skip Pfam assignment  (Default: no)

   -b|-block-size: block size for diamond against the nr database (Default: 8)

   -e|-evalue: max evalue for discarding hits diamond run  (Default: 1e-03)

   -miniden: identity perc for discarding hits in diamond run  (Default: 50)

   

 Binning:

   --nobins: Skip all binning  (Default: no)

   --nomaxbin: Skip MaxBin binning  (Default: no)

   --nometabat: Skip MetaBat2 binning  (Default: no)

   

 Performance:

   -t: Number of threads (Default: 12)

   -canumem: memory for canu in Gb (Default: 32)

   --lowmem: run on less than 16Gb of memory (Default:no)

 

 Other:

   --minion: Run on MinION reads (use canu and minimap2) (Default: no)

   

 Information:

   -v: Version number  

   -h: This help 

     

 

データベースの準備

オフライン環境になる前にダウンロードしておく。

200GBほどスペースが必要。

mkdir data_dir
perl scripts/preparing_databases/download_databases.pl data_dir

#condaで導入した場合
download_databases.pl data_dir

 ダウンロードされ、data_dir/db/に解凍される。

> ls -alth *

db:

total 120G

drwxr-xr-x 4 root root 4.0K Mar 20 06:47 ..

drwxr-xr-x 10 1000 1000 4.0K Feb 26 07:55 .

-rw-r--r-- 1 1000 1000 56 Feb 26 07:55 DB_BUILD_DATE

drwxr-xr-x 2 1000 1000 4.0K Feb 26 07:55 LCA_tax

-rw-r--r-- 1 1000 1000 1.4G Feb 26 06:18 Pfam-A.hmm

-rw-r--r-- 1 1000 1000 3.3G Feb 26 06:17 eggnog.dmnd

-rw-r--r-- 1 1000 1000 112G Feb 26 06:13 nr.dmnd

-rw-r--r-- 1 1000 1000 2.7G Feb 26 00:17 keggdb.dmnd

-rw-r--r-- 1 1000 1000 130M Feb 24 21:18 taxdb.btd

-rw-r--r-- 1 1000 1000 14M Feb 24 21:18 taxdb.bti

-rw-r--r-- 1 1000 1000 1.1K Feb 22 20:51 nr.pal

-rw-r--r-- 1 1000 1000 65M Nov 16 23:54 arc.all.faa

-rw-r--r-- 1 1000 1000 1.3M Nov 16 23:54 arc.scg.faa

-rw-r--r-- 1 1000 1000 198K Nov 16 23:54 arc.scg.lookup

-rw-r--r-- 1 1000 1000 16M Nov 16 23:54 bac.all.faa

-rw-r--r-- 1 1000 1000 244K Nov 16 23:54 bac.scg.faa

-rw-r--r-- 1 1000 1000 36K Nov 16 23:54 bac.scg.lookup

-rw-r--r-- 1 1000 1000 335 May 21 2018 ReadMe

-rw-r--r-- 1 1000 1000 820K Apr 27 2018 arc.hmm

-rw-r--r-- 1 1000 1000 790K Apr 27 2018 bac.hmm

-rw-r--r-- 1 1000 1000 1015K Apr 27 2018 euk.hmm

-rw-r--r-- 1 1000 1000 446K Apr 27 2018 mito.hmm

-rw-r--r-- 1 1000 1000 3.6M Apr 27 2018 bacar_marker.hmm

-rw-r--r-- 1 1000 1000 16M Apr 27 2018 marker.hmm

-rw-r--r-- 1 1000 1000 5.1K Apr 27 2018 .dmanifest

drwxr-xr-x 4 1000 1000 4.0K Jan 14 2015 genome_tree

drwxr-xr-x 2 1000 1000 4.0K Jan 14 2015 img

drwxr-xr-x 2 1000 1000 4.0K Jan 14 2015 pfam

-rw-r--r-- 1 1000 1000 80K Jan 14 2015 selected_marker_sets.tsv

-rw-r--r-- 1 1000 1000 21M Jan 14 2015 taxon_marker_sets.tsv

drwxr-xr-x 2 1000 1000 4.0K Jan 14 2015 test_data

drwxr-xr-x 2 1000 1000 4.0K Jan 14 2015 hmms

drwxr-xr-x 2 1000 1000 4.0K Jan 14 2015 hmms_ssu

drwxr-xr-x 2 1000 1000 4.0K Jan 14 2015 distributions

test:

total 16K

drwxr-xr-x 4 root root 4.0K Mar 20 06:47 ..

drwxrwxr-x 3 1007 1008 4.0K May 21 2018 .

drwxrwxr-x 2 1007 1008 4.0K May 21 2018 raw

-rw-r--r-- 1 1007 1008 156 May 21 2018 test.samples

 

 

 

実行方法

 3つのモードがある。テストデータを使い説明する。

data_dir/test/raw/

f:id:kazumaxneo:20190610171133j:plain

2セットのfastqが用意されている。

 

1、Sequential mode

すべてのサンプルを個別に処理し、順次分析する。 このモードはビニングを含まない。

perl scripts/SqueezeMeta.pl -m sequential -t 40 \
-f data_dir/test/raw/ -s data_dir/test/test.samples
  • -m  <sequential, coassembly, merged>: Mode (REQUIRED)
  • -s    Samples file (REQUIRED)
  • -f     Fastq read files' directory (REQUIRED)
  • -t     Number of threads (Default: 12)

 

2、Coassembly mode

すべてのサンプルからのリードがプールされ、まとめてアセンブリが実行される。 次に、個々のサンプルからのリードをCoアセンブリマッピングして、各サンプル量を取得し、 ビニングしてゲノムビンを得る。メモリ要件が高すぎてクラッシュする場合、3のmerged modeを使う。プロジェクト名が必要。

perl scripts/SqueezeMeta.pl -m coassembly -p projectname -t 40 \
-f data_dir/test/raw/ -s data_dir/test/test.samples
  • -p <string>: Project name (REQUIRED in coassembly and merged modes)

 

3、Merged mode

このモードでは無制限数のサンプルのCoアセンブリが可能なっている。 まずサンプルを個々にアセンブリし、得られたコンティグを単一の共アセンブリにマージする。 その後、Coassemblyモードの場合と同様にビニングが行われる。 キメラコンティグを作成する可能性がより高いので、これは推奨される手順ではない(可能であればcoassembly mode使用)。プロジェクト名が必要。

perl scripts/SqueezeMeta.pl -m merged -p projectname -t 40 \
-f data_dir/test/raw/ -s data_dir/test/test.samples

 

 

 

dockerイメージもあげましたが、ランの最後にクラッシュします。非推奨です。

docker pull kazumax/squeezemeta 

#ホストのカレントディレクトリとイメージの/dataをシェアして起動(pullを飛ばして以下を実行してもOK)
docker run -itv $PWD:/data/ kazumax/squeezemeta

> source ~/.profile
> cd ~/SqueezeMeta

引用

SqueezeMeta, A Highly Portable, Fully Automatic Metagenomic Analysis Pipeline
Javier Tamames, Fernando Puente-Sánchez

Front Microbiol. 2018; 9: 3349

 

ショートリードのマッピングを行う Whisper

 

 リファレンスゲノムへのリードのマッピングは、シークエンシングデータ解析パイプラインの最初のステップである。シーケンシングコストが削減していることから、合理的な時間内に増大する量の生成データを処理することができるアルゴリズムに対する必要性が高まっている。
 著者らはリードをソートし、それからリファレンスゲノムとそのreverse complementのsuffix arraysに対してそれらをマッピングするという考えに基づいて、正確で高性能なマッピングツールであるWhisperを紹介する。 タスクとデータの並列処理を採用し、一時データをディスクに格納すると、妥当なメモリ要件で優れた時間効率が得られる。 Whisperは、大規模なNGSリードコレクション、特にIlluminaの一般的なWGSカバレッジのデータ処理で優れた性能を示す。 実データによる実験では、よく知られているBWA-MEMおよびBowtie 2ツールで必要とされる時間の約15%で、同等の精度でこのソリューションが機能することが示されている。

 

 

本体 Github

リリースよりダウンロードする

https://github.com/refresh-bio/Whisper/releases

#linux v1.1
wget https://github.com/refresh-bio/Whisper/releases/download/v1.1/whisper
chmod +x whisper
sudo mv whisper /usr/local/bin/

#linux
wget https://github.com/refresh-bio/Whisper/releases/download/v1.1/whisper-index
chmod +x whisper-index
sudo mv whisper-index /usr/local/bin/

whisper

$ whisper

Whisper v. 1.1 (2018-07-10)

Usage:

   whisper [options] <index_name> @<files> 

   whisper [options] <index_name> file_se 

   whisper [options] <index_name> file_pe1 file_pe2

Parameters:

  index_name   - name of the index (as created by asm_pp)

  files        - name of the file containing list of FASTQ files with seq. reads

  file_se      - FASTQ file (single-end)

  file_pe[1|2] - FASTQ files (paired-end)

Options:

  -b <value> - no. of temporary files (minimum: 100, default: 384)

  -d[fr/ff/rf] - mapping orientation (default: -dfr (forward - reverse)

  -dist_paired <value> - max. distance for paired read (default: 1000)

  -e <value> - max. no of errors (default: 4, max: 5% of read length)

  -e-paired <value> - max. fraction of errors in paired read (default: 0.06)

  -enable-boundary-clipping <value> - enable clipping at boundaries when a lot of mismatches appears (default: 0)

  -filter <value> - store only mappings for given chromosome (default: )

  -gap-open <value> - score for gap open (default: -6)

  -gap-extend <value> - score for gap extend (default: -1)

  -gzipped-SAM-level <value> - gzip compression level of SAM/BAM, 0 - no compression (default: 0)

  -hit-merging-threshold <value> - minimal distance between different mappings (default: 12)

  -high-confidence-sigmas <value> - (default: 4)

  -hit-merging-wrt-first <value> - calculate distance in marged group w.r.t. first (default: 1)

  -m[f/s/a] - mode: first stratum/second stratum/all strata (default: first stratum)

  -mask-lqb <value> - mask bases of quality lower than value (default: 0)

  -out <name> - name of the output file (default: whisper)

  -penalty-saturation <value> - no. of sigmas for max. penalty in matching pairs (default: 7)

  -rg <read_group> - complete read group header line, (example: '@RG\tID:foo\tSM:bar')

                     '\t' character will be converted to a TAB in the output SAM/BAM,

                      while the read group ID will be attached to every read 

  -r[s|p] - single or paired-end reads <value> (default: single)

  -score-discretization-threshold (default: 0.5)

  -score-match <value> - score for matching symbol (default: 1)

  -score-clipping <value> score for clipping (default: -10)

  -score-mismatch <value> - score for mismatching symbol (default: -5)

  -sens <value> - turn on/off sensitive mode (default: 1)

  -sens-factor <value> - sensitivity factor (default: 3)

  -stdout - use stdout to store the output

  -store-BAM - turn on saving in BAM

  -t <value> - no. of threads (0-adjust to hardware) (default: 0)

  -temp <name> - prefix for temporary files (default: ./whisper_temp_)

  -x - load complete suffix arrays in main memory (default: 0)

Examples:

  whisper human @files

  whisper -temp temp/ human reads1.fq reads2.fq

  whisper -out result.sam -temp temp/ -t 12 human reads1.fq reads2.fq

kazu@kazu:~$ 

 

> whisper-index

$ whisper-index 

Whisper index construction v. 1.1 (2018-07-10)

Usage: 

   whisper-index <index_name> <ref_seq_file_name> <dest_dir> <temp_dir>

   whisper-index <index_name> <@ref_seq_files_name> <dest_dir> <temp_dir>

 

 

実行方法

たとえば、圧縮されていないFASTQ形式のサンプル読み取りのサイズが100GBで、処理が16個のCPUスレッドによって行われる場合、Whisperは600のファイルを開く。そのため、同時に開けるファイル数を上げる必要がある。linuxなら

ulimit -n 600

fastqサイズや使用スレッド数が増えると増加する。 必要な数は

num_files = num_bins (384 by default) + num_threads + 2 * total_size_of_FASTQ_in_GB

で算出できる。

 

1、indexing

出力と作業ディレクトリを指定する。ここでは作業ディレクトリは/tmpか/var/tmpあたりにする。

mdkir output_dir
whisper-index genome_index ref.fa output_dir /tmp

 

2、mapping

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

whisper -out output -t 8 output_dir/genome_index \
pair_1.fq.gz pair_2.fq.gz

output.samが出力される。

 

引用

Whisper: read sorting allows robust mapping of DNA sequencing data
Sebastian Deorowicz Agnieszka Debudaj-Grabysz Adam Gudyś Szymon Grabowski
Bioinformatics, Volume 35, Issue 12, June 2019, Pages 2043–2050

 

 

 

メタバーコディングのデータベース配列キュレーションなどを行うツールキット MetaCurator

 

 配列ベースの生物学的コミュニティの特徴付けの過程において、配列の教師ありのtaxonomic classification は重要な目標である。多数の配列分類ソフトウェアプログラムは、配列類似性を測り、そして配列類似性と分類学的所属との間の関係をモデル化することによってこれを達成する。(一部略)

 配列組成に基づく分類の間、関心のあるマーカーに隣接する無関係の配列領域は、SINTAX、RDPナイーブベイズ分類器、VSEARCHなどの多くのツールで使用されるk-mer分布を偏らせる(Wang et al、2007; Edgar 2016; Rognes et al、 2017)。上位N個の配列の距離分布を分析することによって分類学的境界を推定する分類器(Huson et al、2011; Bengtsson-Palme et al、2015; Gao et al、2017)では、配列重複の保持は距離分布の特徴付けに偏りがある。ただし、2つの配列が同じ分類群からのものであることを保証せずに分類群的に単純な配列重複の除去を行うと、 2つの分類群間の変異に関してデータベースから重要な情報が削除されることになる。
 これらの問題を解決するために、MetaCuratorツールキットを開発した。メインツールであるIterRazorは、例えば全ゲノムアセンブリを含む任意の入手可能な供給源から標的マーカーリファレンス配列を同定および抽出する。与えられたマーカーのリファレンスを抽出した後、DerepByTaxonomyは分類学を意識したアプローチを使ってリファレンス配列を逆複製するために使われる。追加のPythonツールおよびUnixシェルスクリプトは、階層的キュレーションのための分類系統のフォーマット化およびNCBIヌクレオチドデータベース内に一般的に見られる分類系統アーチファクトの除去を容易にする(NCBI Resource Coordinators、2018)。

 

 

インストール

ubuntu16.04にて、conda createで2.7.14の仮想環境を作ってテストした。

依存

MetaCurator is a command-line only toolkit which runs on typical Unix/Linux environments. It requires the following software to be installed and globally accessible.

  • Python 2.7.5 or greater
  • Perl 5.16.0 or compatible
  • VSEARCH 2.8.1 or compatible
  • HMMER3 3.1b2 or compatible
  • MAFFT 7.270 or compatible
#仮想環境を作って依存を導入 
conda create -n metacurator -c bioconda python=2.7.14 vsearch hmmer mafft

 本体 Github

wget https://github.com/RTRichar/MetaCurator/archive/MetaCurator_v1.0beta.1.tar.gz
tar xzf MetaCurator_v1.0beta.1.tar.gz
cd MetaCurator-MetaCurator_v1.0beta.1/
export PATH=$PATH:$PWD

IterRazor.py -h

# IterRazor.py -h

usage: IterRazor is a tool which uses hidden Markov models to identify and extract precise DNA barcode sequences of interest from the full diversity of reference sequence data available through sources such as NCBI (e.g. whole genome sequences, whole chloroplast sequences or Sanger seuqnces representative of variable portions of the locus of interest)

       [-h] -r REFERENCESFILE -i INPUTFILE -o OUTPUTFILE [-t THREADS]

       [--SaveTemp SAVETEMP] [-e HMMEVALUE] [-is ITERATIONSERIES]

       [-cs COVERAGESERIES]

 

optional arguments:

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

  -r REFERENCESFILE, --ReferencesFile REFERENCESFILE

                        Fasta file of sequences trimmed to exact region of

                        interest (header format should conform to --InputFile

                        header format). This should be 5 to 10 sequences

                        spanning a diversity of the phylogenetic tree of

                        interest, so as to build a relatively broadly

                        inclusive HMM from the first iteration (REQUIRED)

  -i INPUTFILE, --InputFile INPUTFILE

                        Input file of reference sequences to be extracted

                        (e.g. whole chloroplast genome sequences from NCBI).

                        Sequence headers should include a unique sequence

                        identifier, NCBI accessions are recomended, preceded

                        only by '>' (REQUIRED)

  -o OUTPUTFILE, --OutPutFile OUTPUTFILE

                        Output file for extracted sequences (REQUIRED)

  -t THREADS, --threads THREADS

                        Number of processors to use

  --SaveTemp SAVETEMP   Option for saving alignments and HMM files produced

                        during extraction

  -e HMMEVALUE, --HmmEvalue HMMEVALUE

                        Evalue threshold for hmm search

  -is ITERATIONSERIES, --IterationSeries ITERATIONSERIES

                        The IterationSeries and CoverageSeries arguments

                        dictate the number of iterative searches to run and

                        the minimum HMM coverege required for each round of

                        extracting. This argument consists of a comma-

                        separated list specifying the number of searches to

                        run during each round of extraction. The list can vary

                        in length, changing the total number of extraction

                        rounds conducted. However, if a search iteration fails

                        to yield any new reference sequence amplicons,

                        IterRazor will break out of that round of searching

                        and move to the next. By default, the software

                        conducts four rounds of extraction with 20, 10, 5 and

                        5 search iterations per round (-is 20,10,5,5).

  -cs COVERAGESERIES, --CoverageSeries COVERAGESERIES

                        This argument consists of a comma-separated list

                        specifying the minimum HMM coverage required per round

                        for an extracted amplicon reference sequence to be

                        added to the reference sequence fasta. The list can

                        vary in length but must be equal in length to the

                        IterationSeries list. By default, the software

                        conducts four rounds of extraction with coverage

                        limits of 1.0, 0.95, 0.9 and 0.65 for each respective

                        search round (-cs 1.0,0.95,0.9,0.65).

> MetaCurator.py -h

# MetaCurator.py -h

usage: MetaCurator.py [-h] -r REFERENCESFILE -i INPUTFILE -it INTAX -ot OUTTAX

                      -of OUTFASTA [-t THREADS] [--SaveTemp SAVETEMP]

                      [-e HMMEVALUE] [-tf TAXONOMIZRFORMAT] [-ct CLEANTAX]

                      [-is ITERATIONSERIES] [-cs COVERAGESERIES] [-mh MAXHITS]

                      [-di ITERATIONS]

 

optional arguments:

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

 

required arguments:

  -r REFERENCESFILE, --ReferencesFile REFERENCESFILE

                        Fasta file of sequences trimmed to exact region of

                        interest (header format should conform to inputfile

                        header format). This should be approximately 10

                        sequences spanning a diversity of the phylogenetic

                        tree of interest, so as to build an appropriately

                        inclusive HMM from the first iteration.

  -i INPUTFILE, --InputFile INPUTFILE

                        Input file of reference sequences to be extracted

                        (sequence headers should include a unique sequence

                        identifier, NCBI accessions are recomended, preceded

                        only by '>')

  -it INTAX, --InTax INTAX

                        Input taxonomic lineages (a tab-delimited file with

                        the first column containing unique sequence

                        identifiers compatible with the files provided under

                        '-r' and '-i'). Normally, this would be in

                        Metaxa2-compatible format, however, if lineages come

                        directly from Taxonomizr, you can skip reformatting

                        and declare '-tf True.' See below for more details.

  -ot OUTTAX, --OutTax OUTTAX

                        Filename for output taxonomic lineages

  -of OUTFASTA, --OutFasta OUTFASTA

                        Filename for output fasta formatted sequences

 

universal optional arguments:

  -t THREADS, --threads THREADS

                        Number of processors to use

  --SaveTemp SAVETEMP   Option for saving alignments and HMM files produced

                        during extraction and dereplication (True or False)

  -tf TAXONOMIZRFORMAT, --TaxonomizrFormat TAXONOMIZRFORMAT

                        Specify if tax lineages are in Taxonomizr format and

                        must first be converted to Metaxa2-compatible tab-

                        delimited format (True or False)

  -ct CLEANTAX, --CleanTax CLEANTAX

                        Specify if tax lineages are to be cleaned of common

                        NCBI artifacts and revised at unresolved midpoints

                        (True or False)

 

IterRazor optional arguments:

  -e HMMEVALUE, --HmmEvalue HMMEVALUE

                        Evalue threshold for HMM search

  -is ITERATIONSERIES, --IterationSeries ITERATIONSERIES

                        The IterationSeries and CoverageSeries arguments

                        dictate the number of iterative searches to run and

                        the minimum HMM coverege required for each round of

                        extracting. This argument consists of a comma-

                        separated list specifying the number of searches to

                        run during each round of extraction. The list can vary

                        in length, changing the total number of extraction

                        rounds conducted. However, if a search iteration fails

                        to yield any new reference sequence amplicons,

                        IterRazor will break out of that round of searching

                        and move to the next. By default, the software

                        conducts four rounds of extraction with 20, 10, 5 and

                        5 search iterations per round (i.e. '-is 20,10,5,5').

  -cs COVERAGESERIES, --CoverageSeries COVERAGESERIES

                        This argument consists of a comma-separated list

                        specifying the minimum HMM coverage required per round

                        for an extracted amplicon reference sequence to be

                        added to the reference sequence fasta. The list can

                        vary in length but must be equal in length to the

                        IterationSeries list. By default, the software

                        conducts four rounds of extraction with coverage

                        limits of 1.0, 0.95, 0.9 and 0.65 for each respective

                        search round (i.e. '-cs 1.0,0.95,0.9,0.65').

 

DerepByTaxonomy-specific optional arguments:

  -mh MAXHITS, --MaxHits MAXHITS

                        Max hits to find before ending search for any given

                        query using VSEARCH (default: 1000)

  -di ITERATIONS, --Iterations ITERATIONS

                        Number of iterative VSEARCH searches to run (default:

                        10). This is important in case the number of

                        replicates for a given sequence greatly exceeds '--

                        MaxHits'

 

 

テストラン

1、メタバーコーディングする領域を指定するため手動でトリミングして集めてきた8つのrbcL塩基配列"rbcL_Reps.fa"を入力とし、NCBIでrbcLとアノテートされている数千のrbcL塩基配列"rbcL_sample.fa"をキュレーションする。

cd TestMetaCurator/
MetaCurator.py -r rbcL_Reps.fa -i rbcL_sample.fa -it rbcL_sample.tax -tf True -ct True -of Test.fa -ot Test.tax

元のデータベースrbcL_sample.fa平均サイズは850.9-bpで2500配列あったが、キュレーション後は1691配列、平均476.6-bpに減少した。

> seqkit stats rbcL_sample.fa Test.fa

# seqkit stats rbcL_sample.fa Test.fa 

file            format  type  num_seqs    sum_len  min_len  avg_len  max_len

rbcL_sample.fa  FASTA   DNA      2,500  2,127,215       94    850.9    8,413

Test.fa         FASTA   DNA      1,691    805,958      319    476.6      494

 Taxnomy情報を表すTest.tax も出力される。

 

上ではNCBIに登録されているrbcL配列群rbcL_sample.fa を-iで指定しているが、実際のデータ解析(メタバーコーディングのプライマー作成、増幅領域のマルチプルアラインメントなど)で使用するときは、ゲノム、メタゲノムのアセンブリなどを-iで指定する。

引用

MetaCurator: A hidden Markov model-based toolkit for extracting and curating sequences from taxonomically-informative genetic markers

Rodney T. Richardson, Douglas B. Sponsler, Harper McMinn-Sauder, Reed M. Johnson

bioRxiv preprint first posted online Jun. 18, 2019

 

 

ヒトのガン原遺伝子/腫瘍抑制遺伝子の変異を視覚化するwebツール Mutplot

 

 シーケンシング技術開発はガン研究に革命をもたらした。約20年に及ぶ発展後、次世代シーケンシング(NGS)は速くて手頃な価格になっている。それは精密医療を臨床の現実にした。 NSGは、臨床現場での治療法を個別化し、研究情報を広げるための包括的なビッグデータを提供している。この技術的進歩は治療および研究のためのより多くの機会を生み出したが、それらは非常に大きく詳細であるため、結果として得られたデータを効率的に合成および要約するという問題も生じた。ビッグデータを手動でフィルタリングすると、エラーの可能性が高まり、整理するのに時間がかかる。ビッグデータも効果的に提示するのは困難である。ソフトウェアはこれらの問題をすべて回避する。この目的のためにいくつかのツールが利用可能である。ただし、ほとんどはプログラミングの背景を持つユーザー向けに設計されている。これは病院やそのような訓練を受けていない施設利用者の大部分を締め出す。 MutplotはWebブラウザで機能する機能し、簡単なカスタマイズのための柔軟性を提供する。臨床医や研究者が自分で使用するために特別に設計された。抽象的なビッグデータを視覚的な結果に変換する。さらに、Mutplotはすべてのプラットフォームで動作するオープンソースのツールであり、セキュリティ目的のため、ファイアウォールの内側で簡単に統合できる。

 Mutplotを、MutationMapper [ref.1]、Lollipops [ref.2]、Muts-needle-plot [ref.3]、Pfam [ref.4]、Plot Protein [ref.5]、trackViewer [ref.6]など、他の6つの最も一般的な突然変異プロットツールと比較した。それらのどれも技術的でないユーザーのためのすべての要件を満たさない(詳細は論文表1に示されている)。 MutationMapperを除くそれらすべては、プログラミングトレーニングを必要とするコマンドラインユーザーインターフェイスを使用する。 Muts-needle-plot、Plot Protein、TrackViewerでは、手動のドメイン入力が必要である。Lollipopsは、サンプル頻度が類似している突然変異とクラスタ化突然変異を区別することはできない。その上、手作業でデータを入力することは人的ミスを犯しがちであり、そしてそれは突然変異ハイライト機能を持たない。 Pfamは公開フォーマットではないJSONファイルを出力する。 MutationMapperはWebベースのユーザーインターフェースを使用しているので最良の選択と思われるが、それ自身の欠点もある。それは最も高い再発性の突然変異(アミノ酸の変化)を示すだけで、これは低い頻度であるがドライバー遺伝子の突然変異を排除するだろう[ref.7]。事実、腫瘍間の遺伝的異質性のために、多くのドライバー遺伝子が非常に低いvariant allele frequenciesで生じる。同じ遺伝子に複数の突然変異が発生した場合、MutationMapperは、ガン研究と個別化医療の進歩に不可欠な低発生突然変異を簡単に無視することがあり得る。さらに、2つのバリアントがMutationMapper内で非常に近くにあると、それらのうちの1つからの情報が重複する。 MutationMapperのもう1つの落とし穴は、スペースが限られている場合にドメイン名が自動的に切り捨てられることである。これらの欠点により、MutationMapperはNGS分析にはあまり理想的ではない。

 Mutplotには、さまざまなタンパク質の変異を可視化するための完全なワークフローが含まれている(論文図1)。バリアント情報(必要な4つの列はHugo_Symbol、Sample_ID、Protein_Change、およびMutaiton_Type。表S1参照)でファイル(タブ区切り形式またはカンマ区切り形式であること)を入力すると、Mutplotは自動的にUniProtから最新のタンパク質情報に接続する[ref.8]。ドロップダウンメニューを使用して、合計409のガン遺伝子および腫瘍抑制遺伝子が組み込まれている(論文表S2)。 Mutplotは選択された遺伝子のドメイン情報を取得する。アミノ酸頻度閾値のハイライトオプションは1、2、3、4、5、10、15、20、25、30に設定されている。両方の遺伝子とハイライト閾値オプションは、ソースコードをカスタマイズするだけで拡張できる。指示はGitHubにdepositされている:https://github.com/VivianBailey/Mutplot。
この情報を使って、Mutplotはドメイン情報、アミノ酸の位置、変異の頻度、アミノ酸の変化、変異の種類、そして強調表示された変異を含むタンパク質の図を作成する。アミノ酸の位置は正確な比率のために遺伝子の長さに合わせて調整されている。強調表示された変異は詳細なアミノ酸改変情報を有する。突然変異の種類と説明は、視覚化と区別を容易にするために色分けされている。複数の変異が密集している場合、Mutplotはスマートに他の変異を妨害することなく変異を標識する方法を見つけ出す。 Mutplotは出力オプションに関しても高い柔軟性を与える。 JEPG、PDF、PNG、画像ダウンロード用のSVGをサポートしている。また、入力データと更新されたUniprotデータベースから取得した対応するドメイン情報から、ダイアグラム用に選択したデータをエクスポートすることもできる。

このソースコードGitHub: https://github.com/VivianBailey/Mutplotで利用かのであり、非営利目的で使用することができる。そして簡単にアクセスしたり、修正したり、他のパイプラインやソフトウェアに統合することができる。ソースコードを修正すると、MutplotはWebアプリケーションからファイアウォールの内側にあるパーソナルコンピュータまたはサーバに移行する可能性がある。これは、厳格なセキュリティ規制に従う機関にとって大きな選択肢となる。さらに、GitHubにはMutplotの完全なドキュメント、ソースコードのカスタマイズ方法の説明、そして将来のリリースもGitHubに説明が載っている。

WebアプリはRプログラミング言語で開発された。 shinyパッケージ、ggplot2パッケージ、plyrパッケージ、httrパッケージ、drawProteinsパッケージ、およびggrepelパッケージが使用されている。

 

Github

  

 使い方

https://bioinformaticstools.shinyapps.io/lollipop/ にアクセスする。

f:id:kazumaxneo:20190619160133j:plain

 

入力ファイルはタブ/カンマ区切りの4列からなるフォーマットで用意する必要がある。それぞれの列には、Hugo_Symbol、Sample_ID、Protein_Change、およびMutaiton_Typeが記載されている必要がある。

入力ファイル例(direct link)。

f:id:kazumaxneo:20190619213507p:plain

 

Browseからファイルをアップロードする。

f:id:kazumaxneo:20190619213710p:plain

ヘッダーがある場合は「ヘッダー」を選択。入力ファイルに合わせセパレータを指定する。

 

視覚化する遺伝子を409のリスト中から選択。ここではテストデータを使っているので、TP53遺伝子か、BRCA1を選ぶ。

f:id:kazumaxneo:20190619160137j:plain

 利用できる遺伝子リスト(direct link)。

 

 

ハイライト表示するアミノ酸頻度のしきい値を設定する。

f:id:kazumaxneo:20190619221044p:plain

defaultでは全変異のハイライト表示は無しになっているが、数値を選ぶことで、多くのサンプルに共通する変異に絞ってハイライト表示できる。上に載せた例では、R175Hの変異が5サンプルで共通しており、他はそれぞれ1サンプルのみから由来。

 

ここでは1を選択、出力フォーマットにはpdfを指定。

f:id:kazumaxneo:20190619221615p:plain

Download Plotで図がダウンロードされる。

 

f:id:kazumaxneo:20190619221443p:plain

 

図の右端にレジェンドが表示される。変異のタイプの他、モチーフ位置も表現されている。

f:id:kazumaxneo:20190619221832p:plain

 

ハイライト表示するアミノ酸頻度のしきい値を5に変更

f:id:kazumaxneo:20190619221428p:plain

5サンプル以上で共通する変異のみハイライト表示された。

 

 

 

引用

Mutplot: An easy-to-use online tool for plotting complex mutation data with flexibility
Zhang W, Wang C, Zhang X

PLoS One. 2019 May 15;14(5)

 

 

 

 

バクテリア、アーキア、プラスミドの複製起点(ori)データベース DoriC

2019 6/21 誤字修正、コマンド修正

2023/10/19 URL修正

 

 すべての生物において、DNA複製は複製機構の構築段階で正確に制御されている(ref.1)。複製起点は特定のゲノム遺伝子座であり、そこでは二本鎖DNAがほどけて一本鎖DNA鋳型を形成して新しい鎖の合成を開始する。大部分の細菌において、複製起点(oriC)は、主要なイニシエータータンパク質DnAによって認識されるいくつかのDnaAボックスモチーフを含む。また複製起点(oriC)にはAT含有量の高いDNA unwinding element (DUE)を含む領域、ここでは一本鎖DNAもDnaAにより認識される(ref.2〜5)、も含まれる 。ATリッチなDNA unwinding elementは古細菌複製起点にも不可欠であることがわかっており、これはorigin recognition proteins (ref.6,7)の結合部位として働くorigin recognition boxes (ORBs) に隣接している。多数のプラスミドにおいて、origin of vegetative replication (oriV) はしばしばダイレクトリピートまたはiteron DNA配列(wiki)からなり、これらはRepタンパク質と相互作用して複製開始の過程で初期複合体を形成する(ref.8)。oriVのiteronの位置の近くにATに富む領域もあり、これはDNA巻き戻し要素として働く(ref.9)。

 複製起点が通常、dnaA、orc1 / cdc6およびrep遺伝子などの複製関連遺伝子の隣にあることは興味深い。chromosomeおよびプラスミド上の原核生物複製起点の類似構造は、同じ枠組みに基づいて起源予測のためのアルゴリズムを設計する機会を提供する。当初、Ori-Finderは細菌のchromosome上のoriC領域を同定するために開発された(ref.10)。

 3つのドメインの独立したドメインの1つとして、ほとんどの古細菌は地球上の様々な極端な環境に存在しており、特定のhabitsが実験的方法による複製起源の同定を困難にしている(ref.11)。したがって、ウェブベースのツールOri-Finder2は、古細菌ゲノムのoriC領域をin silicoで予測するために開発されたものであり、予測された結果は実験室での古細菌起源の同定に役立つ可能性がある(ref.12)。

 プラスミドは、chromosome外の自己複製する遺伝要素であり、細菌、古細菌酵母、そしていくつかの高等真核細胞に広く見られる(ref.8)。プラスミドはしばしば抗生物質耐性やtoxin–antitoxinシステムのような宿主細胞にいくつかの特別な特徴をもたらす遺伝子を持っている(ref.13)。したがって、プラスミドの自律的なDNA複製は細胞の生存にとって重要である。Vegetative replicationの起源はプラスミドの最も重要な要素の一つである。これまで、oriVの位置と特徴は、RK2、F、P1、R6K、pPS10プラスミドなどの広範囲のプラスミドでよく理解されていた(ref.14–18)。しかしながら、シーケンシングされた多数のプラスミド上でoriVを自動的に同定するためにバイオインフォマティクスツールが緊急に必要とされている。

 複製起点に関する関連研究を容易にするために、Ori-Finderシステムの予測がオンラインデータベースにまとめられた(ref.19–21)。 2007年に、oriC領域のデータベースであるDoriCが最初に公に利用可能になり、2013年に、DoriC 5.0は細菌ゲノムと古細菌ゲノムの両方の複製起点を含めた(ref.22、23)。過去6年間で、次世代シーケンシング技術の急速な進歩と様々な微生物ゲノムプロジェクトからのシーケンシングされたゲノムの蓄積がDoriCの拡大を促進しており、この拡大されたデータベースはDoriC 10.0としてここに提示される。はじめてのプラスミドDoriCデータベースとOri-Finderシステムは、複製起点の構造と機能のより良い理解を確実にし、そしてDNA複製における開始段階の調節メカニズムへの新しい洞察を提供する。これまでのところ、DoriCに保存されている予測の多くは現在実験室で検証されており、過去のDoriCデータベースとOri-Finderシステムに基づくより多くのアプリケーションが我々の最近のarticleで見直された(ref.24, pubmed)。

 今回のリリースでは、DoriCの内容は次のようにバージョン5.0と比較して大幅に改善されている。(i)細菌chromosome上のoriCは4倍の1633から7580に増加した。 (ii)古細菌のchromosome上のoriCは86から226に増加した。 (iii)NCBI recordから検索された348の注釈付き起点および修正Ori-Finderシステムによって861の予測起点を含む1209のプラスミド複製起点が初めて提示された。 (iv)originの機能に重要な新規リピートトリヌクレオチドモチーフである DnaA-trio要素を含む、細菌複製起点中のより多くの配列要素が組み込まれる。 DnaA-trioは、DoriCによるバイオインフォマティクス分析によって細菌界全体で高度に保存されている、origin の巻き戻しおよびDNAヘリカーゼローディングにおいて役割を果たす(ref.20)。 DnaA-trio様配列をDoriCデータベースで検索した後、その情報を対応するoriCレコードに追加した。さらに、データベースのユーザーインターフェースをより便利で直感的にわかるように再設計した(論文図1)。(以下略)

 

 

使い方

DoriCにアクセスする。

version12はこちら

 

以下、古いバージョンでの説明

f:id:kazumaxneo:20190620140653j:plain

 Browseからbacteira、Archaea、plamidのいずれかを選択する。

f:id:kazumaxneo:20190620215431p:plain

 

登録されている配列が表示される。

f:id:kazumaxneo:20190620215623p:plain

 

DoriC accession numberをクリックすることで、ポジションなどの詳細を表示できる。

f:id:kazumaxneo:20190620215740p:plain

 

Z-curves (bacteira、Archaeaのみ)

f:id:kazumaxneo:20190620221938p:plain

 

Browse in NCBI

f:id:kazumaxneo:20190620220620p:plain

 

BLAST検索はテスト時動作しなかった。

 

downloadからはRefseq IDや配列を含むCSVファイルをダウンロードできる。

http://tubic.org/doric/public/index.php/download

f:id:kazumaxneo:20190620222748p:plain

解凍。bacteira、Archaea、plamid3つのCSVがある。

f:id:kazumaxneo:20190620222731p:plain

excelで開いた。

f:id:kazumaxneo:20190620223326p:plain

 

おまけ

ori配列を使ってアセンブルを支援する。

awkCSVファイルからfasta形式に変換する。変換後、EMBOSSパッケージのseqret (紹介) を使い、エラーを除きつつ適当な文字数で改行。

#ubuntuで実行
#bacteira
awk 'BEGIN{FS=","}{ OFS = "" }{print ">", $3, "\n", $NF}' tubic_bacteria.csv \
| sed '1,2d' - > bacteira_oriC.fasta

#archaea
awk 'BEGIN{FS=","}{ OFS = "" }{print ">", $3, "\n", $NF}' tubic_archaea.csv \
| sed '1,2d' - > archaea_oriC.fasta

#plasmid
awk 'BEGIN{FS=","}{ OFS = "" }{print ">", $3, "\n", $NF}' tubic_plasmid.csv \
| sed '1,2d' - > plasmid_ori.fasta

#emboss seqret
seqret archaea_oriC.fasta archaea_oriC_corrected.fasta
seqret bacteira_oriC.fasta bacteira_oriC_corrected.fasta
seqret plasmid_ori.fasta plasmid_corrected.fasta

*$NFで最後のカラムのみ出力。OFS = ""でスペースを排除。sedに渡して先頭のコメント行(awk処理で2行になっている)を排除。EMBOSSのseqretはゲノム登録用に開発されており、潜在的なエラーをあらかた排除できる。fold -w "INT" やperl/awkよりオススメ。

 

De novoアセンブル支援に使う。

BandageにGFA/fastgを読み込みblast検索する。

f:id:kazumaxneo:20190621010852p:plain

完全ではないが、メタゲノムのアセンブリから、アーキアバクテリア、プラスミドを見分けるのに役立つと思われる。

引用

DoriC 10.0: an updated database of replication origins in prokaryotic genomes including chromosomes and plasmids
Hao Luo, Feng Gao
Nucleic Acids Research, Volume 47, Issue D1, 08 January 2019, Pages D74–D77

 
DoriC 5.0: an updated database of oriC regions in both bacterial and archaeal genomes

Gao F, Luo H, Zhang CT

Nucleic Acids Res. 2013 Jan;41(Database issue):D90-3


DoriC: a database of oriC regions in bacterial genomes

Gao F, Zhang CT

Bioinformatics. 2007 Jul 15;23(14):1866-7. Epub 2007 May 12

 

2023/10/19

DoriC 12.0: an updated database of replication origins in both complete and draft prokaryotic genomes
Mei-Jing Dong, Hao Luo, Feng Gao

Nucleic Acids Res. 2023 Jan 6;51(D1):D117-D120. 

 

参考

bioinformatics - Identifying the origin of replication of an unannotated *E. coli* plasmid - Biology Stack Exchange

ショートリードによるpolishingも行う高速なロングリードアセンブラ Raven (旧名 Ra)

2020 5/23 タイトル補足、ravenインストール追記

2020 8/11 引用にpreprint追記

2021 5/24 論文引用

2022/06/08 help更新

 

 Ra(現在はRaven)は、第3世代シーケンシングによって生成されたrawシーケンシングリードの高速で使いやすいアセンブラである。 以下の図に示すように、RaはMinimap2、Rala、およびRaconで構成されている。

 Raは入力としてFASTA / FASTQフォーマット(gzipで圧縮可能)のrawシーケンシングリードを含む単一のファイルを取り、高精度の一連のコンティグをstdoutにFASTAフォーマットで出力する。 さらに、FASTA / FASTQ形式の第2世代シーケンスリードファイル(gzip圧縮対応)を第2引数として受け取り、ショートリードでpolishingすることで最終的なアセンブリを完成させることができる。

 

 

f:id:kazumaxneo:20190618010005p:plain

Ra flow chart. Githubより

参考スライド ( Rayan Chikhi, CNRS, Univ Lille BiG talk, Lund University)

http://rayan.chikhi.name/pdf/big18_large_genome_assembly.pdf

http://rayan.chikhi.name/pdf/big18_large_genome_assembly.pdf

とても面白い内容です。ワクワクしますね。 Raについては最後の方で触れられています。

 

2021/11/29

 

インストール

ubuntu16.04のminiconda3.4.0.5環境でテストした(docker使用、ホストOS ubuntu18.04)。

依存

  • gcc 4.8+ or clang 3.4+
  • cmake 3.2+

Github

 


#--recursiveをつけて依存するminimap2、racon、そしてralaモジュールも含めてクローンする。
git clone --recursive https://github.com/rvaser/ra.git ra
cd ra
mkdir build
cd build
cmake -DCMAKE_BUILD_TYPE=Release ..
make -j 8
cd bin/

> ./ra --help

# ./ra --help

/usr/bin/env: 'python': No such file or directory

root@ebaf20ee92c1:~/ra/build/bin# source ~/.profile 

root@ebaf20ee92c1:~/ra/build/bin# ./ra --help

usage: ra [-h] [-u] [-t THREADS] [--version] -x {ont,pb}

          sequences [ngs_sequences]

 

positional arguments:

  sequences             input file in FASTA/FASTQ format (can be compressed

                        with gzip) containing third generation sequences for

                        assembly

  ngs_sequences         input file in FASTA/FASTQ format (can be compressed

                        with gzip) containing next generation sequences for

                        polishing (default: None)

 

optional arguments:

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

  -u, --include-unused  output unassembled and unpolished sequences (default:

                        False)

  -t THREADS, --threads THREADS

                        number of threads (default: 1)

  --version             show program's version number and exit

 

required arguments:

  -x {ont,pb}           sequencing technology of input sequences (default:

                        None)

Raven

#bioconda
mamba create -n raven -y
conda activate raven
mamba install -c bioconda raven-assembler -y

#from source
git clone --recursive https://github.com/lbcb-sci/raven.git raven
cd raven && mkdir build && cd build
cmake -DCMAKE_BUILD_TYPE=Release .. && make
./bin/raven

> ./raven 

usage: raven [options ...] <sequences> [<sequences> ...]

 

  # default output is to stdout in FASTA format

  <sequences>

    input file in FASTA/FASTQ format (can be compressed with gzip)

 

  options:

    -k, --kmer-len <int>

      default: 15

      length of minimizers used to find overlaps

    -w, --window-len <int>

      default: 5

      length of sliding window from which minimizers are sampled

    -f, --frequency <double>

      default: 0.001

      threshold for ignoring most frequent minimizers

    -p, --polishing-rounds <int>

      default: 2

      number of times racon is invoked

    -m, --match <int>

      default: 3

      score for matching bases

    -n, --mismatch <int>

      default: -5

      score for mismatching bases

    -g, --gap <int>

      default: -4

      gap penalty (must be negative)

    --graphical-fragment-assembly <string>

      prints the assembly graph in GFA format

    --resume

      resume previous run from last checkpoint

    --disable-checkpoints

      disable checkpoint file creation

    -t, --threads <int>

      default: 1

      number of threads

    --version

      prints the version number

    -h, --help

      prints the usage

 

 

 

実行方法

ONTのアセンブリ

ra -t 20 -x ont long_reads.fq > output.fa

GFAファイルとアセンブリされた配列output.faが出力される。

 

Pacbioのアセンブリ

ra -t 20 -x pb long_reads.fq short_reads.fq > output.fa

 

ONTのリードとショートリードのハイブリッドアセンブリ。ペアエンドショートリードデータはあらかじめ結合しておく。

ra -t 20 -x ont long_reads.fq short_reads.fq > output.fa

 

引用

Raven: a de novo genome assembler for long reads

Robert Vaser, Mile Sikic

bioRxiv, Posted August 10, 2020

 

2021 5/24

Time- and memory-efficient genome assembly with Raven | Nature Computational Science

 

論文中のパフォーマンス比較結果をみると多くのmetricsについてバランスの取れた結果を出しており、Ravenが優秀なアセンブラであることが分かりますね。

 

 

 

関連


 

補足 

同名のショートリード用アセンブリツールが別にあります。