macでインフォマティクス

macでインフォマティクス

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

原核生物ゲノムのアセンブリ結果およびアノテーションを改善するためのWebプラットフォーム ReNoteWeb

 

 DNA塩基配列の解読にかかる費用と時間が短縮されたことにより、NCBI(National Center for Biotechnology Information)のような公開データベースへの生物情報の寄託が大幅に増加した。1回の実行で大量のデータが生成されるため、この新しい特徴を持つデータを扱い、原核生物ゲノムのアセンブリアノテーションなどの解析を支援することができるアルゴリズムを開発する必要性に迫られることになった。長年にわたり、この作業を自動化するパイプラインや計算ツールが開発され、その結果、与えられた生物、特に非モデル生物の遺伝的内容を知るための総時間を短縮し、バイオテクノロジーに応用可能なターゲットの同定と連携することができるようになった。自動アノテーションツールの場合、その結果の正確さは文献で広く認められているが、自動プロセスで推測された情報をキュレーターが検証し、充実させるという手動キュレーションプロセスを排除するものではない。この作業には、研究対象生物の遺伝子産物の数に正比例する時間が必要である。この作業を支援するために、シンプルで直感的なインターフェースを備えたウェブツールReNoteWebを提供し、オリジナルのゲノム配列に欠けている産物を特定する可能性を持つアセンブリ強化作業を行う。さらに、ReNoteWebは、高精度な外部データベースから取得した情報をもとに、すべてのプロダクトに対してアノテーション処理を行うことができる。データ処理を担うエンジンはJAVAで開発され、ウェブプラットフォームはYiiフレームワークのリソースを利用している。このプラットフォームで作成されるアノテーションは、手作業によるキュレーションプロセスの全体的な時間を短縮することを目的としている。本ツールの検証には23の生物を使用した。NCBIデータベースで利用可能なこれらの同じ生物のアノテーションとRASTプラットフォームで実行されたアノテーションを比較することにより、効率が検証された。本ツールは、http://biod.ufpa.br/renoteweb/ で利用できる。

 

Guide

http://biod.ufpa.br/renoteweb/files/Guide_ReNoteWeb.pdf

 

webサービス

http://biod.ufpa.br/renoteweb/にアクセスする。

プロジェクトを閲覧したり、新しくプロジェクトを作成することができる。

 

1つプロジェクトを作ってみる。Create Projectタブを選択する。

プロジェクト名、メールアドレス、実行したい内容を選択する。ここではAnnotationを選択。

 

結果はSearchタブでkeyを入力することでダウンロードできる。そのkeyは記載したメールアドレスにジョブを投げたタイミングで送られている。

improved assembly

アセンブリfastaファイルを指定し、次の画面でfastqを指定する(fastq, fastq.tar.xz, fq, fq.tar.xz)。tar.xzに圧縮するなら"tar -Jcf R1.fq.tar.xz R1.fastq"など。

試したfastqデータでは、tar.xzで固めるとtar.gzの2/3くらいのファイルサイズになった。

 

追記

混雑しているのか1日経っても結果が返ってこないので、また結果が出てきたら追記します。

 

2022/10/07追記

自分の環境では正しくランすることができませんでした。

引用

ReNoteWeb - Web platform for the improvement of assembly result and annotation of prokaryotic genomes

Gislenne da Silva Moia, Antônio Sérgio Cruz Gaia, Mônica Silva de Oliveira, Victória Cardoso Dos Santosa, Jorianne Thyeska Castro Alves, Pablo Henrique Caracciolo Gomes de Sá, Adonney Allan de Oliveira Veras

Gene. 2022 Nov 30;844:146819

 

参考


全自動のトランスポーザブル・エレメントのアノテーションと解析のパイプライン Earl Grey

 

 トランスポーザブル・エレメント(TE)は、ほぼ全ての真核生物ゲノムに存在し、様々な進化過程に関与している。TEに関する研究は非常に盛んだが、そのアノテーションと特性解析は、特に非専門家にとって依然として困難である。(i)断片的で重複するTEアノテーションは、TE数およびカバレッジの誤った推定につながる可能性がある。(ii)リピートモデルは、5'および3'領域の捕捉が不十分で、全長の小さな割合しか表さないことがある。また、既存のパイプラインは、インストール、実行、データの抽出が困難な場合がある。これらの問題を解決するために、本著者らはEarl Greyを開発した。このパイプラインは、真核生物のゲノムアセンブリに含まれるTEのキュレーションとアノテーションをユーザーフレンドリーな形で行うために設計された全自動トランスポーザブル・エレメント・アノテーション・パイプラインである。

 シミュレーションゲノム、3つのモデルゲノムアセンブリ、3つの非モデルゲノムアセンブリを用いた結果、Earl Greyは、現在広く用いられているTEアノテーション手法よりも優れており、非冗長化TEライブラリでより長いTEコンセンサス配列を生成し、それを用いて重複のない断片的なTEアノテーションを生成することにより、上記の問題を改善することができた。Earl Greyは、TEアノテーションベンチマーク(MCC:0.99)および分類(97%の正解率)において、既存のソフトウェアと比較して高いスコアを獲得した。
 Earl Greyは、包括的かつ完全に自動化されたTEアノテーションツールキットであり、他のバイオインフォマティクスツールと互換性のある標準フォーマットで、研究者に論文用の要約図とアウトプットを提供する。Earl Greyはモジュール式であり、今後のリリースにおいて、品質管理面や独自の解析に特化したモジュールを追加することが可能である。

 

インストール

ubuntu18に導入した。

Github

#Repeatmaskerなどの依存するツールを導入(YAML)
git clone https://github.com/TobyBaril/EarlGrey
cd ./EarlGrey
mamba env create -f earlGrey.yml
conda activate earlGrey
./configure #*1
chmod +x ./earlGrey

> ./earlGrey 

____________________

< Checking Parameters >

 --------------------

        \   ^__^

         \  (oo)\_______

            (__)\       )\/                

                ||----w |

                ||     ||

    #############################

    earlGrey

    Required Parameters:

        -g == genome.fasta

        -s == species name

        -o == output directory

        -r == RepeatMasker search term (e.g arthropoda/eukarya)

 

    Optional Parameters:

        -t == Number of Threads (DO NOT specify more than are available)

        -l == Repbase species subset library (FASTA format)

        -i == Number of Iterations to BLAST, Extract, Extend - Only enter a number between 5 and 10 (Iteration outputs are 0-based!) (Default: 5)

        -f == Number flanking basepairs to extract (Default: 1000)

        -d == Maximum non-TE distance between annotations to consider a pair to be from the same cluster (100 is recommended)

    -h == Show help

 

    Example Usage:

 

    earlGrey -g bombyxMori.fasta -s bombyxMori -o /home/toby/bombyxMori/repeatAnnotation/ -r arthropoda -t 16

 

 

    Prerequisites - These must be configured prior to using Earl Grey:

        - RepeatMasker (Version 4.1.2)                                                                                                                                                                                                               - RepeatModeler (version 2.0.2)                                                                                                                                                                                                                                                                                                                                                                                                                                                       Ensure you have run the ./configure script from the earlGrey installation directory

    Ensure RepeatMasker has been configured with the desired repeat libraries (RepBase and Dfam 3.2 are recommended)

 

    Queries can be sent to:

    tb529@exeter.ac.uk

 

    Please make use of the GitHub Issues and Discussion Tabs at: https://github.com/TobyBaril/EarlGrey

    #############################

 

 

実行方法

ゲノム配列を指定する。

earlGrey -g genome.fasta -s species_name -o repeatAnnotation/ -r arthropoda -t 16
  • -s    species name
  • -r    RepeatMasker search term (e.g arthropoda/eukarya)
  • -g   genome.fasta
  • -o   output directory

 

  • 最も重要な結果はsummaryFiles/に保存される。中間結果も手動でキュレーションや調査を行う事に備えて保持される(各ファイルについてはレポジトリで説明されている)。

引用

Earl Grey: a fully automated user-friendly transposable element annotation and analysis pipeline
Tobias Baril,  Ryan M. Imrie,  Alex Hayward

bioRxiv, Posted July 02, 2022

 

関連


*1

mamba install -c bioconda repeatmasker -y

mamba install -c bioconda repeatmodeler -y

 

BUSCO遺伝子を使った系統解析のsnakemakeワークフロー

 

2022/09/13 追記、誤字修正

2024/01/16 タイトル修正、わかりにくい説明を修正

 

Current Protocols in BioinformaticsにBUSCOの使い方に関する論文が出ています。(引用1)。その論文のパートの1つ;”Support Protocol 3: BUILDING PHYLOGENOMIC TREES”では、BUSCOがシングルコピーの系統マーカー遺伝子を取り出すことを利用して、ゲノム配列から遺伝子セットを取し出し、複数の遺伝子から系統推定するプロトコルについて説明されています。BUSCOの標準のコマンドではサポートしていない方法のため、この論文では、snakemakeベースでいくつかの外部ツールとスクリプトを組み合わせてこれを行っています。手順を確認してみます。

 

インストール

依存

  • snakemake
  • conda

この系統推定パイプラインでは、busco_protocol/support_protocol3/envsにある3つのenvbiroment.yamlの環境がアクティベートされる(リンク)。これらの環境のツールは、初回ランで自動導入される。

git clone https://gitlab.com/ezlab/busco_protocol
cd busco_protocol/support_protocol3/

Snakefileとconfig.yamlが確認できる。

 

config.yamlファイル

config.yamlにはgenomeのfastaファイルのパス、GFF形式のアノテーションファイルのパス、出力、BUSCO5.2.2に対応したBUSCOデータセットなどを書く。3つの環境のYAMLファイルへの相対パスも書かれている。これらの環境は、BUSCOのラン、Rのスクリプトのラン、コア遺伝子の抽出とMSA&トリミングなどで使い分けられる(リンク)。

genomesがgenomeのfastaファイルを含むディレクトリへのパス、GFFsがGFF形式のアノテーションファイルを含むディレクトリへのパス、resultsが出力ディレクトリ、geno_suff: ".fna" ならfastaファイルの拡張子は.fna、insecta_odb10がBUSCOデータセット(参考;v4リンク)。BUSCOデータセットは使うゲノムの分類によって変える必要がある。分類がよく分からないゲノムなら、マーカー遺伝子数は少ないがfungi_odb10のような大きな分類を使う手もある。percent_missing: 0だとすべてのゲノムに存在。keep_duplicates: "No"ならduplicates BUSCOに相当するマーカー遺伝子は使われない(デフォルトでNo)。

 

準備
上ではgenomeとなっているので、genome/に系統推定したいゲノムのfastaファイル(.fna)、GFFs/にアノテーションのgffファイルを配置する(対応するfastaとgffのファイル名は拡張子以外一致している必要がある)。

さらに、カレントにrulesディレクトリ、scriptsディレクトリ、envsディレクトリ、ルールの中にもenvディレクトリを配置する必要がある。本テストデータではこのままで問題ない。

 

 

Snakefileファイル

configfileではカレントパスにあるconfig.yamlが指定されている。Snakefileの中でincludeして呼び出されるSnakefileファイルは相対パスでrules下に含まれている。

 

最後に呼び出されてiqtreeを実行するSnakefile;rules/run_iqtree.smkが以下。

iqtree2のランで、アウトグループを-o Zootermopsis_nevadensis”としているが、これは論文中でgenomes/に配置されたfastaファイルの1つ;Zootermopsis_nevadensis.fnaをアウトグループにしているため。自分のデータでは修正するか、このオプション自体を消す必要がある。他のパラメータについても同様にチェックしておく。

 

 

ラン

準備ができたら実行する。30スレッド指定。マーカー遺伝子のトリミングやMSAで標準エラー出力が膨大に出てくるので、ファイルに保存するか捨てる。

snakemake --cores 30 --use-conda 2> error.log

 

それぞれのマーカータンパク質配列はdata/に、連結タンパク質のMSAや系統推定結果はresults/直下に保存される。

 

results/

phylogeny.treefileが最終出力。6つ見えるサブディレクトリはゲノム毎のBUSCO5.2.2出力。

(論文で使用されている6つの昆虫ゲノムの結果)。

 

data/extracted/にはゲノム毎にディレクトリができる。中には、それぞれのマーカータンパク質配列が保存されている。

 

smakemakeで実行されるフローはシンプルなため、解析手順に興味がある方は、rulesに含まれている8つのファイルを開いて確認してください。修正したいパラメータ設定があればこのファイルを書き換えるだけです。簡単ですね(その場合、コマンドだけでなく、どこを変えたかのも記録も残しておきましょう)。

 

この論文はオープンアクセスです。詳しくは論文を読んで下さい。

 

2024/01/15追記 ~tidyverseの導入失敗に伴うエラー~

SP3.yamlからcondaの仮想環境が作られるが、SP3-2.yamlのR3.6のtidyverseの導入に失敗するため、ラン途中でエラーが発生するようになってしまっている。対策として、SP3-2.yamlを開いてRのバージョンを3.6から4.2に書き換えて実行したがそれでもsnakemakeによる自動インストールで失敗していた。そこで環境をconda activateでアクティベート後(コマンド実行時のエラーログに長い数値からなる環境名をアクティベートしたというメッセージが残っている。おそらく最初に発生したエラーの少し前のlogあたり)、"mamba install r::r-tidyverse -y"でtidyverseを導入してdeactivate後、それから実行することで回避した。

引用
1

BUSCO: Assessing Genomic Data Quality and Beyond
Mosè Manni, Matthew R Berkeley, Mathieu Seppey, Evgeny M Zdobnov

Curr Protoc. 2021 Dec;1(12):e323

 

関連


参考

 

 

 

(ヒト)バリアントの機能的アノテーションリソース FAVOR

 

 大規模な全ゲノムシークエンシング(WGS)研究とバイオバンクにより、多数のコーディングおよびノンコーディングバリアントが急速に生成されている。これらは、ヒト疾患の遺伝的基盤を明らかにするための前例のないリソースを提供する。バリアント機能アノテーションは、WGS解析、結果の解釈、疾患や形質に関連する原因バリアントの優先順位付けにおいて重要な役割を担っている。既存の機能アノテーションデータベースは、オンラインクエリを実行する範囲が限られていたり、大規模なWGS研究およびバイオバンクの遺伝子型データを下流の解析のために機能的にアノテーションすることができなかったりすることがある。本著者らは、このような緊急のニーズに応えるために、Functional Annotation of Variants Online Resources (FAVOR)を開発した。FAVORは、ゲノム上に存在する90億個の一塩基バリアント(SNV)の要約と可視化、およびバリアント、遺伝子、地域レベルのオンラインクエリーを可能にする多面的な総合オンラインポータルを提供する。また、複数のソースからのバリアント機能情報を統合してバリアントの機能特性を記述し、ヒトの表現型に影響を与える原因となりうるバリアントの優先順位付けを容易にする。さらに、大規模なシーケンス研究の遺伝子型とバリアントの機能アノテーションデータを機能的にアノテーションし、アノテーションされたGDSファイル形式で効率的に保存し、下流の解析を容易にする拡張性のあるアノテーションツール、FAVORannotatorが提供されている。FAVORおよびFAVORannotatorは、https://favor.genohub.org で入手できる。

 

Github

(HPより)FAVORannotatorは、FAVORデータベースを用いて全ゲノム/全エキソームシーケンス(WGS/WES)研究のオフライン機能アノテーションを行うためのオープンソースRプログラム( GitHub )です。機能アノテーションデータと入力された遺伝子型データ(VCF)を結合し、オールインワンのaGDSファイルを作成します。aGDSフォーマットでは、遺伝子型と機能アノテーションの両方のデータを1つのファイルに格納することができます。(補足;FAVORannotatorのデータベースはhttps://favor.genohub.org/のFAVORannotatorタブからダウンロードできる。)

 

webサービス

https://favor.genohub.org/にアクセスする。

バリアントを入力する。ポジション、rsID、遺伝子名、領域などに対応している。

 

出力例

Summary

結果の説明とAnnotationリソースについてはDocumentationタブに記載されています。

 

Full tables

 

Figures

 

 

Batch annotationタブでは、入力ファイルをアップロードすることで、より大規模なクエリを実行することができる。

 

sign upタブからユーザー登録しておくとアップデート情報を受け取れるようです。

引用

FAVOR: Functional Annotation of Variants Online Resource and Annotator for Variation across the Human Genome
Hufeng Zhou, Theodore Arapoglou, Xihao Li, Zilin Li, Xiuwen Zheng, Jill Moore, Abhijith Asok, Sushant Kumar, Elizabeth E. Blue, Steven Buyske, Nancy Cox, Adam Felsenfeld, Mark Gerstein, Eimear Kenny, Bingshan Li, Tara Matise, Anthony Philippakis, Heidi Rehm, Heidi J. Sofia, Grace Snyder, NHGRI Genome Sequencing Program Variant Functional Annotation Working Group, Zhiping Weng, Benjamin Neale, Shamil R. Sunyaev, Xihong Lin

bioRxiv, Posted August 29, 2022

 

関連


バリアントコーリングを自動化する柔軟でスケーラブルなパイプライン grenepipe

 

 本著者らは、個体や集団のハイスループットな生シーケンスデータから遺伝子型バリアントコールまでのデータ処理を効率化するオールインワンSnakemakeワークフローであるgrenepipeを開発した。このパイプラインは、一般的なソフトウェアツールを単一の設定ファイル内で提供し、ソフトウェアの依存関係を自動的にインストールし、クラスタ環境でのスケーラビリティに高度に最適化され、単一のコマンドで実行できる。grenepipeはGPLv3の下で公開されており、github.com/moiexpositoalonsolab/grenepipeで自由に利用することができる。

 

Wiki

https://github.com/moiexpositoalonsolab/grenepipe/wiki

Quick Start and Full Example

https://github.com/moiexpositoalonsolab/grenepipe/wiki/Quick-Start-and-Full-Example

 

インストール

依存

  • conda(mamba含む)
  • snakemake == 6.0.5
  • python == 3.7.10
  • pandas == 1.3.1
  • numpy == 1.21.2

Github

git clone https://github.com/moiexpositoalonsolab/grenepipe.git
cd grenepipe/
mamba env create -f envs/grenepipe.yaml -n grenepipe
conda activate grenepipe

 

 

実行方法

1、fastqの場所とサンプル名の関係を示したsample.tsvを用意する。フォーマットは、1列名はsample名、2列名は1回のシークエンシングデータなら1(複数ある場合は、下の画像のサンプルBのように2つ目のfastqには2をつける)、3列目はプラットフォーム名(GATKで影響)、4列名はfastqのパス(ペアエンドならforward fastq)、5列目はペアエンドのみでR2 (reverse) のfastqのパス。

https://github.com/moiexpositoalonsolab/grenepipe/wiki/Setup-and-Usage#samples-table

ここではsample.tsvとして保存する。

補足;tools/generate-table.pyを使うと、このサンプルTSVを自動で作成できる。fastqを含むディレクトリを指定する。

tools/generate-table.py <path>/<to>/fastq_dir

 

config.yamlの中でsample.tsvの場所を指定する。config.yamlのテンプレートはレポジトリの直下に用意されている。line13で作成したsample.tsvのパスを指定する。

https://github.com/moiexpositoalonsolab/grenepipe/blob/master/config.yaml

 

2、このconfig.yamlを結果を保存したい空のディレクトリに配置する。そのディレクトリのパスを指定して実行する。

snakemake --cores 20 --use-conda --directory /path/to/output

ディレクトリのYAMLファイルが認識され、パイプラインがスタートする。その後、envs/にある環境がインストールされるが、このインストールにかなり時間がかかる可能性がある。試したときはqualimap.yamlの環境導入でハングアップしていた(conda側の問題)。

 

引用

grenepipe: A flexible, scalable, and reproducible pipeline to automate variant calling from sequence reads 
Lucas Czech,  Moises Exposito-Alonso
Bioinformatics, Published: 02 September 2022

 

海洋環境ゲノムをマイニングするためのオンラインサービス The Ocean Gene Atlas v2.0

 

 Tara Oceansの海洋メタゲノムやメタトランスクリプトームのような大規模データリソースを用いて遺伝子の生物地理に関する仮説を検証するには、多大なハードウェアリソースとプログラミングスキルが必要になる。今回リリースされた「Ocean Gene Atlas」(OGA2)は、大規模かつ複雑な海洋環境ゲノムデータベースをマイニングするための、自由に利用できる直感的なオンラインサービスである。OGA2のデータセットが拡張され、Tara Oceansのポートフォリオから、(i) 真核生物のMetagenome-Assembled-Genome (MAGs) と Single-cell Assembled Genomes (SAGs) (10.2E+6 coding genes), (ii) Ocean Microbial Reference Gene Catalogue version 2 (46.0E+6 non-redundant genes) 、(iii) 924 MetaGenomic Transcriptomes (7E+6 unigenes), (iv) 530 MAGs from an Arctic MAG catalogue (1E+6 genes) and (v) 1888 Bacterial and Archaeal Genomes (4.5E+6 genes)が追加され、さらにMalaspina 2010世界一周航海からのデータセット: (vi) 317 Malaspina Deep Metagenome Assembled Genomes (0.9E+6 genes)も収録している。OGA2が可能にした新しい分析には、海洋環境データセットとRefSeqデータベースの両方からの配列の相同性という文脈でユーザーのクエリを可視化する系統樹推論が含まれている。アプリケーションプログラミングインターフェース (API) により、ユーザーはコマンドラインツールを使用して OGA2 に問い合わせることができるようになり、ローカルワークフローの統合が可能になった。最後に、遺伝子発現量は、利用可能な環境変数のいずれかを使用して、地図表示上で直接対話的にフィルタリングすることができる。Ocean Gene Atlas v2.0は、https://tara-oceans.mio.osupytheas.fr/ocean-gene-atlas/ で自由に利用できる。

 

help

https://tara-oceans.mio.osupytheas.fr/ocean-gene-atlas/build/pdf/Ocean-Gene-Atlas_User_Manual.pdf

 

webサービス

https://tara-oceans.mio.osupytheas.fr/ocean-gene-atlas/にアクセスする。

 

データベースを選択する。

クエリの配列は、fasta形式の塩基配列(10kb以下)、アミノ酸配列、作成済みのHMMプロファイル、PfamのID、遺伝子名、以前の結果の視覚化(メールアドレスを記入してランするとTSVファイルが入手できる。それを指定する。)に対応している。

E-valueでの閾値はデフォルトで1E-10となっている。HMMプロファイルでのランは相当な時間かかるので、高感度さが重要では無いなら、塩基配列アミノ酸配列による検索が推奨される。

 

ここでは 上のメニューからTry an Example with OM-RGC dataset (prokaryotes)を選択した。

Submitをクリック。

 

出力例

mapが複数あるのは、それぞれのmapに異なる変数を表示させて比較しやすくするため。

 

バブルプロットの変数は図の左上のボタンから変更できる。

 

helpより

  • mapの用語について; DCM: deep chlorophyll maximum layer; SRF: upper layer zone; MES: mesopelagic zone; MIX: marine epipelagic mixed layer, FSW: filtered sea water, ZZZ: marine water layer、をそれぞれ表す。
  • バブルプロットのサンプルの円の大きさは、クエリのホモログの存在量に比例している。円は選択された分数によって色分けされている。Y軸は環境パラメータ値を表す。アルカリ度、アンモニウム5m*、炭素総量、CDOM*、クロロフィル_A、CO2、CO3、密度、深度、沿岸距離、HCO3、鉄5m*、硝酸5m*、亜硝酸5m*、NO2、NO3、NO3_NO2、NPP_C*、O2、PAR、pH、PIC*、PO4、POC*、塩分、Si、温度(submitする時にBubble plotsの数を増やすと出力される)。
  • 海洋モデルから推定された値には星印が付けられている。

 

引用
The Ocean Gene Atlas v2.0: online exploration of the biogeography and phylogeny of plankton genes
Caroline Vernette, Julien Lecubin, Pablo Sánchez, Tara Oceans Coordinators; Shinichi Sunagawa, Tom O Delmont, Silvia G Acinas, Eric Pelletier, Pascal Hingamp, Magali Lescot 

Nucleic Acids Res. 2022 Jun 10;50(W1):W516-W526

 

関連


メタゲノムの株レベルプロファイリングを行う PStrain

2022/09/07, 9/8  追記

 

 微生物群は、人間の病気や生理活動に不可欠な役割を担っている。ゲノム配列の株レベルの違いにより、微生物の機能が異なることがある。ショットガン・メタゲノムシーケンスを用いると、微生物群集の菌株を実用的にプロファイリングすることができる。しかし、菌株間の配列が非常に類似しているため、現在の手法は未発達である。本著者らは、遺伝子型頻度から、同じ一塩基変異(SNV)遺伝子座における菌株の遺伝子型を推定できることを確認した。また、同じリードがカバーする異なる遺伝子座のバリアントは、同一株上に存在する証拠となる可能性がある。
 PStrainは、MetaPhlAn2マーカー遺伝子のSNVに基づき、遺伝子型頻度と複数のSNV座をカバーするリードを利用して、反復的に株のプロファイリングを行う最適化手法である。PStrainは、従来の手法と比較して、菌株の存在量と遺伝子型の推定性能をそれぞれ平均87.75%と59.45%向上させることに成功した。PStrainを大腸ガンの2つのコホートのデータセットに適用したところ、Bacteroides coprocola株の配列がCRCと対照サンプルで有意に異なることがわかった。これは、CRCの腸内細菌叢におけるB.coprocolaの役割の可能性を初めて報告するものである。

 

インストール

依存

Third party pipelines:

  • SamTools-1.3.1
  • bowtie2-2.3.1
  • metaphlan2/3
  • GATK-3.5 (java-1.8.0_241)
  • picard-tools-2.1.0

Github

mamba create -n PStrain python=3 -y
conda activate PStrain
pip install  numpy scipy pandas pulp pysam pickle5
#third party (GATkは登録が必要*1)
mamba install -c bioconda gatk==3.5 -y
mamba install -c bioconda samtools -y
mamba install -c bioconda bowtie2 -y
mamba install -c bioconda picard -y

#本体
git clone https://github.com/wshuai294/PStrain.git
cd PStrain/scripts/

> python3 PStrain.py -h

usage: python3 PStrain.py [options] -c/--conf <config file> -o/--outdir <output directory>

 

The config file should follow the format:

 

//

sample: [sample1_ID]

fq1: [forward reads fastq]

fq2: [reverse/mate reads fastq]

//

sample: [sample2_ID]

fq1: [forward reads fastq]

fq2: [reverse/mate reads fastq]

...

 

Help information can be found by python3 PStrain.py -h/--help, config file format for single end reads , and additional information can be found in README.MD or https://github.com/wshuai294/PStrain.

 

PStrain: profile strains in shotgun metagenomic sequencing reads.

 

required arguments:

  -c , --conf       The configuration file of the input samples.

  -o , --outdir     The output directory.

 

options:

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

  -p , --proc       The number of process to use for parallelizing the whole pipeline (default is 1), run a sample in each process.

  -n , --nproc      The number of CPUs to use for parallelizing the mapping with bowtie2(default is 1).

  -w , --weight     The weight of genotype frequencies while computing loss, then the weight of linked read type frequencies is 1-w. The value is between 0~1. (default is 0.0)

  --lambda1         The weight of prior knowledge while rectifying genotype frequencies. The value is between 0~1. (default is 0.0)

  --lambda2         The weight of prior estimation while rectifying second order genotype frequencies. The value is between 0~1. (default is 0.0)

  --species_dp      The minimum depth of species to be considered in strain profiling step (default is 5).

  --snp_ratio       The SNVs of which the depth are less than the specific ratio of the species's mean depth would be removed. (default is 0.45).

  --qual            The minimum quality score of SNVs to be considered in strain profiling step (default is 20).

  --similarity      The similarity cutoff of hierachical clustering in merge step (default is 0.95).

  --elbow           The cutoff of elbow method while identifying strains number. If the loss reduction ratio is less than the cutoff, then the strains number is determined. Default is 0.24.

  --bowtie2         Path to bowtie2 binary. If your platform is not Linux, you need to specify your own bowtie2 binary.

  --bowtie2-build   Path to bowtie2 binary. If your platform is not Linux, you need to specify your own bowtie2-build binary.

  --samtools        Path to samtools binary (default version is in the packages path). If your platform is not Linux, you need to specify your own samtools binary.

  --metaphlan2      Path to metaphlan2 script (default version is in the packages path).

  --gatk            Path to gatk binary (default version is in the packages path).

  --picard          Path to picard binary (default version is in the packages path).

  --dbdir           Path to marker gene database (default is in the ../db path).

  --prior           The file of prior knowledge of genotype frequencies in the population.

> python3 PStrain_V30.py -h

usage: python3 PStrain.py [options] -c/--conf <config file> -o/--outdir <output directory>

 

The config file should follow the format:

 

//

sample: [sample1_ID]

fq1: [forward reads fastq]

fq2: [reverse/mate reads fastq]

//

sample: [sample2_ID]

fq1: [forward reads fastq]

fq2: [reverse/mate reads fastq]

...

 

Help information can be found by python3 PStrain.py -h/--help, config file format for single end reads , and additional information can be found in README.MD or https://github.com/wshuai294/PStrain.

 

PStrain: profile strains in shotgun metagenomic sequencing reads.

 

required arguments:

  -c , --conf           The configuration file of the input samples.

  -o , --outdir         The output directory.

 

optional arguments:

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

  -p , --proc           The number of process to use for parallelizing the

                        whole pipeline (default is 1), run a sample in each

                        process.

  -n , --nproc          The number of CPUs to use for parallelizing the

                        mapping with bowtie2(default is 1).

  -w , --weight         The weight of genotype frequencies while computing

                        loss, then the weight of linked read type frequencies

                        is 1-w. The value is between 0~1. (default is 0.0)

  --lambda1             The weight of prior knowledge while rectifying

                        genotype frequencies. The value is between 0~1.

                        (default is 0.0)

  --lambda2             The weight of prior estimation while rectifying second

                        order genotype frequencies. The value is between 0~1.

                        (default is 0.0)

  --species_dp          The minimum depth of species to be considered in

                        strain profiling step (default is 5).

  --snp_ratio           The SNVs of which the depth are less than the specific

                        ratio of the species's mean depth would be removed.

                        (default is 0.45).

  --qual                The minimum quality score of SNVs to be considered in

                        strain profiling step (default is 20).

  --similarity          The similarity cutoff of hierachical clustering in

                        merge step (default is 0.95).

  --elbow               The cutoff of elbow method while identifying strains

                        number. If the loss reduction ratio is less than the

                        cutoff, then the strains number is determined. Default

                        is 0.24.

  --bowtie2             Path to bowtie2 binary. If your platform is not Linux,

                        you need to specify your own bowtie2 binary.

  --bowtie2-build       Path to bowtie2 binary. If your platform is not Linux,

                        you need to specify your own bowtie2-build binary.

  --samtools            Path to samtools binary (default version is in the

                        packages path). If your platform is not Linux, you

                        need to specify your own samtools binary.

  --metaphlan3          Path to metaphlan3 (default version is in the packages

                        path).

  --metaphlan3_output_files 

                        If you have run metaphlan3 already, give metaphlan3

                        result file in each line, the order should be the same

                        with config file. In particular, '--tax_lev s' should

                        be added while running metaphlan3.

  --gatk                Path to gatk binary (default version is in the

                        packages path).

  --picard              Path to picard binary (default version is in the

                        packages path).

  --dbdir_V30           Path to marker gene database (default is in the

                        ../db_V30 path).

  --prior               The file of prior knowledge of genotype frequencies in

                        the population. Not providing this is also ok.

 

 

実行方法

ランするにはfastqのパスとサンプル名の関係を示したconfigファイルが必要。

https://github.com/wshuai294/PStrain/blob/master/test/config.txt

 

マーカー遺伝子のDBはオーサーによって準備されている(リンク)。ダウンロードしてディレクトリに配置する。

DB_dir/

 

configファイル、CPU数、出力ディレクトリ、マーカーDB、サードパーティのツールのパスを指定して実行する。

python3 PStrain_V30.py -c config.txt -n 30 -o out  --metaphlan3 <path>/<to>/metaphlan --bowtie2 <path>/<to>/bowtie2 --bowtie2-build <path>/<to>/bowtie2-build --samtools <path>/<to>/samtools  --picard <path>/<to>/picard.jar --gatk <path>/<to>/GenomeAnalysisTK.jar --dbdir_V30  <path>/<to>/<DB_dir>
  • -n     The number of CPUs to use for parallelizing the mapping with bowtie2(default is 1).
  • -o     The output directory
  • -c      The configuration file of the input samples.

 

出力

"-o"で指定した出力ディレクトリには、configファイルで指定されたsample IDそれぞれのサブディレクトリ、全サンプルをマージしたmergeというサブディレクトリができる。

sample IDのサブディレクトリの例

サブディレクトリ中のresultに結果が保存される。

strain_RA.txtにはどの種にどれだけの株が見つかったのか示されている。

strain_RA.txtの例(Githubより)

2列目のSpecies_RA カラムは、サンプル中の種の相対的な存在量、3列目のStrain_ID 列は、サンプル中の各生物種の菌株ID、4列目のCluster_ID カラムは、全サンプル中の菌株の対応するクラスタを意味する。

 

  • mergeというサブディレクトリには、どの種にどれだけの株が見つかったのかの個数がまとめられてたファイル、全サンプル中の各菌株クラスターのコンセンサス配列、全サンプル中の各生物種の株番号が示されたファイル、その他が出力される。
  • 自分で簡単に種固有のパイプラインを構築できる(Gihtub参照)。
  • PStrain-tracerというPStrainの下流解析のパイプラインも開発されている。PStrain-tracerは、全ゲノムショットガンメタゲノミクスデータから、ゲノム配列や株の性質を推定するツールで、Fecal microbiota transplantation (FMT) の分析に使用されている(論文)。
  • PStrain-tracerにはテストスクリプトが付属しており、簡単に動作確認できる(依存関係を揃えて試したところ、問題なく全結果が出力された)。

 

引用

PStrain: an iterative microbial strains profiling algorithm for shotgun metagenomic sequencing data 
Shuai Wang, Yiqi Jiang, Shuaicheng Li Author Notes
Bioinformatics, Volume 36, Issue 22-23, 1 December 2020, Pages 5499–5506

 

関連


*1

GATK3のダウンロード

https://console.cloud.google.com/storage/browser/gatk-software/package-archive/gatk;tab=objects?pli=1&prefix=&forceOnObjectsSortingFiltering=false

ダウンロード後、gatk-registerをバイナリを指定して叩くとコピーされパスが通る。

#初回はgatk-registerを走らせて実行ファイルのパスを叩く
gatk-register <path>/<to>/<your>/GenomeAnalysisTK.jar

その後、PStrainを実行するときはJava Archive(.JARファイル)の方を指定する。Picardも同様。