macでインフォマティクス

macでインフォマティクス

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

全工程が自動化された高速なRNA seq解析webサービス RaNA-Seq(60以上のモデル生物に対応)

2020 4/10 タイトル修正、説明と図追加

 

RaNA-Seqは、RNA-Seqデータを迅速に解析・可視化するためのクラウドプラットフォームである。FASTQファイルの定量、品質管理指標の計算、発現変動遺伝子の解析の実行、機能解析による結果の説明を可能にすることで、数分で完全な解析を実行する。本解析パイプラインは、一般的に受け入れられている再現性の高いプロトコルを適用しており、Webインターフェースは2つの簡単なステップで適用することができる。解析結果は、解釈と公開の準備ができているインタラクティブなグラフィックスとレポートとして表示される。RaNA-Seq ウェブサービスhttps://ranaseq.eu で自由に利用できる。

 我々(著者ら)は、完全かつ信頼性の高いRNA-Seq解析を行うための最適化された計算パイプラインを設計した。論文図1に我々のパイプラインのフローチャートを示す。FASTQファイルはFastpツール(Chen et al., 2018)で前処理され、Salmon(Patro et al., 2017)で発現定量化が行われる。定量化に基づいて、RaNA-Seqはインタラクティブなグラフを含む品質管理レポートを生成する。これらのグラフは、パッケージRJSplot(BarriosおよびPrieto、2018)を用いて生成され、グラフ間の距離はSERE(Schulzeら、2012)を用いて計算される。発現変動遺伝子の解析は、DESeq2、EdgeRおよびlimma(Law et al、2014; Love et al、2014; Robinson et al、2010)を用いて行うことができ、結果として得られたデータは、インタラクティブなレポートとしてD3で表現される。
 これに加えて、本パイプラインでは、各発現変動結果に対して、over-representation analysis と gene set enrichment analysis (GSEA)を行う。エンリッチメント解析は、RパッケージGOseq(Young et al., 2010)を用いて実行される。この方法は、発現変動結果に含まれる長い転写物や高発現の転写物が過剰に検出されることによる、RNA-Seqデータの結果の偏りを回避するものである。GSEAについては、fgsea Rパッケージを統合した(Sergushichev, 2016)。これは、 cumulative statistical calculationsを用いて、事前にランク付けされたGSEAを高速実行する。機能アノテーションデータベースについては、NCBI BioSystemsデータベース(Geer et al., 2009)からGene Ontologyとパスウェイアノテーションを統合したアノテーションデータベースを作成した。

 我々(著者ら)は、バイオインフォマティクス手法の経験のないユーザーにも自信を持ってRNA-Seqデータを解析できるようにすることを目的として、RaNA-Seqを設計した。理想的なRNA-Seq解析インターフェースが持つべきであるPoplawskiら(2016)によって記述された主な推奨事項に従った。補足表S1は、他のRNA-Seq解析プラットフォームと比較して、我々のプラットフォームが最も完成度の高い解析アプリケーションの一つであることを示している。

 パフォーマンスの高さはRaNA-Seqのもう一つの重要な特徴である。完全なRNA-Seq解析は、2つの簡単なステップで数分で実行される(論文補足図S1および補足表S2)。本解析には、FASTQファイルの定量、リードトリミング、サンプルの品質管理、発現変動、 functional enrichment解析、gene set enrichment analysis(GSEA)が含まれている。現在、RaNA-Seqは、10のヒトRNA-Seqサンプルのプロジェクトに対して、平均30分でこれらすべての処理を実行することができ、これは著者らの経験上、最速のRNA-Seq解析ウェブプラットフォームである(論文補足図S2の実行時間のガントチャートを参照)。

 

 

f:id:kazumaxneo:20200408094448p:plain

RaNA-SEQ Workflow.  HPより転載

 

RaNA-SEQ Tutorial 

 

 

helpからマニュアルPDFをダウンロードできます。結果の解釈の仕方についても簡単に説明されています。読んで下さい。

https://ranaseq.eu/help

 

webサービス

https://ranaseq.eu/ にアクセスする。

f:id:kazumaxneo:20200408094255p:plain

 

Analysisタブに移動。アノニマスモードで使用する事もできるが、登録が推奨されている。

f:id:kazumaxneo:20200408094341p:plainその場でアカウントは作成され、すぐにログインできるようになる。

 

1、サンプルのアップロード

解析の流れを確認していく。

RaNA-Seqは、Reference genomeを指定し、シーケンシングリード(fastq)をアップロードするだけでRNA-seq解析を行うことができる。

f:id:kazumaxneo:20200407213325p:plain

 

まずリファレンスゲノムを指定する。インクリメンタルサーチに対応しているので、タイプして絞り込む。

f:id:kazumaxneo:20200407213522p:plain

 

リファレンスゲノムがない場合、問い合わせると登録してもらえる。

f:id:kazumaxneo:20200407213251p:plain

(?をクリックすると右下にhelpが表示されます)

 

 

 

ペアエンドがシングルエンドかを指定して、gzip形式の圧縮fastqをアップロードする(拡張子は.fastq.gzになっている事)。まずペアエンドかシングルエンドか選択し、

f:id:kazumaxneo:20200408094053p:plain

Advancedをクリックしてライブラリタイプを選択する。

f:id:kazumaxneo:20200410095248p:plain

 

ファイルを指定するか、ウィンドウ内にドラッグ&ドロップすると、

f:id:kazumaxneo:20200407213830p:plain

 

リストに登録されていく。全ファイルを一度にドラッグ&ドロップしてO.K。ペアエンドは2つのファイルを一度にドラッグ&ドロップしてもペアで登録される。

f:id:kazumaxneo:20200408093919p:plain

時々認識されないことがある。リストに出ない場合はドラッグ&ドロップをやり直す。

 

 アップロードのゲージが進行していくので、全サンプルの アップロードが終わるまで待つ。

f:id:kazumaxneo:20200408095006p:plain

 アップロードに失敗しているサンプルはfastqの右にxが付く。この状態だと左端のチェックマーク✔︎をONにできない。右端のボタンをクリックして消し、もう一度アップロードする。

f:id:kazumaxneo:20200410094332p:plain

 

select allを押して全選択した。

f:id:kazumaxneo:20200408101734p:plain

補足

アップロードが完了すると、そのままマッピングが開始される。リスト右端にあるウオッチマークをクリックすると、サンプルの品質を確認できる。

f:id:kazumaxneo:20200408101825p:plain

 

補足

手持ちのfastq以外に、SRA/ENAのプロジェクトのダイレクトアップロードにも対応している。

f:id:kazumaxneo:20200410100246p:plain

プロジェクトのidentifierを打ち込む。

 

サンプルのペアを確認して解析をスタートさせる。

 

f:id:kazumaxneo:20200410103920p:plain

いくつかのレビューでは、それぞれの研究デザインに合わせて解析を行う必要があると結論づけている。解析ツールを変更して発現解析を繰り返すことを推奨する(マニュアルより)。

f:id:kazumaxneo:20200410104447p:plain

 

 

 

 

2、結果

終わった項目から閲覧できるようになる。

f:id:kazumaxneo:20200409222751p:plain

 

 

Quantification

Salmonを使ったリード定量値が出力される。

f:id:kazumaxneo:20200409224024p:plain

定量結果はTSVやexcel形式でダウンロードできる。

f:id:kazumaxneo:20200410111217p:plain

 

QC

f:id:kazumaxneo:20200409223822p:plain

 f:id:kazumaxneo:20200409223830p:plain

f:id:kazumaxneo:20200409223827p:plain

f:id:kazumaxneo:20200409223836p:plain


DE results

f:id:kazumaxneo:20200410111403p:plain

f:id:kazumaxneo:20200410111412p:plain

 

f:id:kazumaxneo:20200410111453p:plain

 

Functional enrichment

f:id:kazumaxneo:20200409222648p:plain

 

Gene OntologyとPathwayが選べる。

f:id:kazumaxneo:20200409224140p:plain

GOはAmiGO2ともリンクしている。

 

pathwayはKEGGとBIOCYCを利用している。

f:id:kazumaxneo:20200409222633p:plain
KEGGのデータはKEGG pathwayにリンクしている。

 

Barplot

f:id:kazumaxneo:20200410091956p:plain


Bubbleplot

f:id:kazumaxneo:20200410091918p:plain

 

Network

f:id:kazumaxneo:20200410091637p:plain

Symmetric Heatmap

f:id:kazumaxneo:20200410091633p:plain

TriangleとSquareから選べる。

 

GSEA (GOとpathway)

f:id:kazumaxneo:20200409224355p:plain

 

RUG plots 

f:id:kazumaxneo:20200410091310p:plain

Network

f:id:kazumaxneo:20200410091436p:plain



 

Symmetric Heatmap

f:id:kazumaxneo:20200410091320p:plain

TriangleとSquareから選べる。

感想

fastqからスタートするため、IDの互換性なども気にせず結果を出せる素晴らしいサービスです。ただ、fastqのアップロードがやや不安定で、アップロードが完了していてもペアエンドの片方しか認識しないなどの現象が発生しました。私だけかもしれませんが、アップロードは注意して下さい。使用したウェブブラウザかデータに問題があったのかもしれません(macos10.14.6のsafari最新版使用)。

引用

RaNA-Seq: interactive RNA-Seq analysis from FASTQ files to functional analysis
Carlos Prieto, David Barrios
Bioinformatics, Volume 36, Issue 6, 15 March 2020, Pages 1955–1956

 

関連

関連ツール

植物遺伝子のファミリー分類とエンリッチメント解析を行うwebサービス GenFam

 

 ハイスループットシーケンシング(HTS)技術を用いたゲノムスケールの研究では、様々な実験条件で発現が異なる遺伝子のリストが作成されている。これらの遺伝子リストは、下流の機能的遺伝学的解析を導くために、生物学的に関連する遺伝子と関連する機能を絞り込むために、さらにマイニングされる必要がある。そのためには、遺伝子オントロジー(GO)の語彙に基づいた遺伝子の機能アノテーションに頼るエンリッチメント解析ツールを用いて、ユーザーが定義したリストの中で統計的に過剰発現している遺伝子を決定するのが一般的なアプローチである。ここでは、遺伝子ファミリーに基づいて遺伝子のアノテーション、分類、エンリッチメントを可能にする新しい計算アプローチであるGenFamを提案する。遺伝学的には、遺伝子の種類ごとに、その遺伝子ファミリーを分類して解析することが可能である。その結果、既存の機能的エンリッチメントツールと比較して、その堅牢性と網羅性の高さが実証された。植物生物学者が利用しやすいように、GenFamはウェブベースのアプリケーションとして提供されており、ユーザーは遺伝子IDを入力し、表形式とグラフ形式の両方でエンリッチメント結果をエクスポートすることができる。また、様々な統計的エンリッチメント検定や多重検定補正法から選択して解析パラメータをカスタマイズすることができる。また、ウェブベースのアプリケーション、ソースコード、データベースは自由にダウンロードして使用することができる。

 GenFamは、よくアノテーションされたシロイヌナズナ(Berardiniら、2015)およびイネ(Oryza sativa)(Kawahara et al、2013)ゲノム、文献検索、およびPfamタンパク質ファミリーデータベース(El-Gebali et al、2019)に基づいて、遺伝子を384の代表的でユニークな遺伝子ファミリーに分類しており、我々(著者ら)の知る限りでは最大のコレクションである。Pfamの共通保存ドメインを同定し、相同遺伝子配列間のドメイン構成を利用して遺伝子ファミリーを割り付けた。これらの高度に保存されたドメインは、タンパク質の機能を定義し、タンパク質をコードする遺伝子を遺伝子ファミリーに分類する。保存されたシグネチャータンパク質ドメインは、配列ベースの類似性解析ツール[例えば、BLAST(Altschul et al、1997)]では困難な発散または遠縁なホモログを検出する能力を有する。したがって、ドメインベースの検索手法は、BLASTベースの相同性検索よりも、遺伝子ファミリーに属する遺伝子をより多く同定することになる。

 GenFamの実装にあたっては、Phytozome(v12)で公開されている植物ゲノムリソースを活用し、背景となるキュレーションされたデータベースを開発した。ユーザー定義の入力リスト内のすべての遺伝子IDをこのリファレンスデータベースにマッピングして遺伝子をファミリーに割り当て、その後、バックグラウンドデータベースと比較することで、入力リスト内の過剰発現遺伝子ファミリーを計算する。また、このデータベースを作成するために、60種類の植物ゲノムのタンパク質配列を用いて、保存されているタンパク質ドメインを同定し、既知遺伝子、未分類遺伝子、新規遺伝子にファミリーを割り当てた。それぞれのタンパク質ドメインは、タンパク質ファミリー隠れマルコフモデル(HMM)プロファイル(Pfamリリース32.0)を用いて、HMMER(v3.1b2)によって予測された(El-Gebali et al、2019)。シグネチャーが保存されたタンパク質ドメインの存在に基づいて遺伝子を分類し、遺伝子ファミリーに割り当てるためのルールを確立した(論文表S1)。このアプローチにより、アノテーションが欠落しているオーファンな遺伝子、アノテーションが正しくない遺伝子、およびそれぞれのゲノムデータベースの間に存在する新規遺伝子を含めて、分類を最大化することができた。また、バックグラウンドデータベースは、ファミリー間での遺伝子メンバーの重複や重複を除去するためにキュレーションを行った。植物ゲノムから384の代表的な遺伝子ファミリーとそれに対応する遺伝子(平均で約41%)をデータベースに統合することができた(論文表S2)。これは、他の既存のデータベースと比較すると、60種の植物にまたがる遺伝子ファミリーの包括的なコレクションである。例えば、最近公開されたポプラの遺伝子ファミリーデータベース(GFDP)では、6,551のポプラ遺伝子をシロイヌナズナゲノム由来の145の遺伝子ファミリーに分類している(Wang et al. PlantTFDB(v4.0)。PlnTFDB(v3.0)では、58および84の転写因子遺伝子ファミリーに分類した(Jin et al、2017; Perez-Rodriguez et al、2010)。同様に、別のデータベースと解析ツールキットであるPlantGSEAは、イネ(118遺伝子ファミリー)やトウモロコシ(81遺伝子ファミリー)のような、主によくアノテーションされたゲノムから遺伝子ファミリーを輸入している13種の植物種の遺伝子ファミリー解析をサポートしている。

 すべての遺伝子ファミリーデータは、PostgreSQLデータベースを用いてフォーマット化し、様々な統計的エンリッチメント手法を用いて分類とエンリッチメント解析を行った。タンパク質ドメインの完全なアノテーションと遺伝子ファミリーの分類を含むGenFamデータベースは、GenFamのウェブサイト(http://mandadilab.webfactional.com/home/dload/)からダウンロードできる。各遺伝子ファミリーに割り当てられた遺伝子数とバックグラウンド遺伝子の総数の詳細な統計は論文表S2に示されている。

 GenFamのWebサーバは、Python3 (https://www.python.org/)、Django 1.11.7 (https://www.djangoproject.com/)、PostgreSQL (https://www.postgresql.org/)のデータベースを使用して実装している。データフォーマットや統計解析のコードはすべてPythonスクリプト言語を使用して実装した。Python は、統計解析、グラフィックス、Web アプリとの統合のためのパッケージがよく開発されている本格的なプログラミング言語である。そのため、GenFamの開発にはRなどの他の言語よりもPythonを選択した。高レベルのPythonウェブフレームワークDjangoを用いて構築した。Python ウェブフレームワークは WebFaction (https://www.webfaction.com/) を使用してホストされた。ウェブベースのテンプレートは、Bootstrap、HTML、CSSを用いて設計した。GenFamは、Internet ExplorerMicrosoft EdgeGoogle ChromeMozillaSafariを含むすべての主要なブラウザと互換性がある。また、分析されたデータは、matplotlib(Droettboom et al、2016)Pythonプロットライブラリを用いて可視化された。

 

Documentation

http://mandadilab.webfactional.com/home/doc/

Github

 

 

 

webサービス

http://mandadilab.webfactional.com/home/にアクセスする。

f:id:kazumaxneo:20200408193407p:plain

 

GeneIDをペーストする。Phytozome ID(Phytozome locus、Phytozome transcript、Phytozome pacID)に対応している。

f:id:kazumaxneo:20200408201203p:plain

Gossypium raimondiiのexample データを貼り付けた。

 

植物種を選択する。

f:id:kazumaxneo:20200408193914p:plain

現在60種類の植物ゲノムをサポートしている。

 

IDタイプ、有意差検定を選択する。Fisher exact test、Hypergeometric distribution、Binomial distribution(二項分布)、Chi-squared distributionの4つの検定をサポートしている。1000 ID未満の小規模なデータセットではFisher exact test,かchi-square test, hypergeometric distributionが推奨され、大規模なデータセットには二項分布を使用することが推奨されている。

f:id:kazumaxneo:20200408201311p:plain

多重比較検定はBonferroni、Bonferroni-Holm、Benjamini-Hochberg検定をサポートしている(参考HP)。

 

Runをクリックして実行する。

 

出力

f:id:kazumaxneo:20200408193314p:plain

Enriched gene familes

P値<0.05でエンリッチされた遺伝子ファミリー

f:id:kazumaxneo:20200408203508p:plain

 

All gene families

ユーザーデータから分類された全ての遺伝子ファミリーの結果

f:id:kazumaxneo:20200408203555p:plain

 

Get Figures

f:id:kazumaxneo:20200408193317p:plain

 

引用.
GenFam: A web application and database for gene family-based classification and functional enrichment analysis

Bedre R, Mandadi K

Plant Direct. 2019 Dec 4;3(12):e00191

 

関連


 

参考

Bonferroni法、Holm法、False Discovery Rate | 大阪大学腎臓内科

tumorサンプルのテロペアリピート数を推定する telomerehunter

 

 テロメアは、真核生物の染色体の末端にある核タンパク質の複合体である。ヒトでは、テロメアDNAは主にノンコーディングのt型(TTAGGG)リピートで構成されているが、c型(TCAGGG)、g型(TGAGGG)、j型(TTGGGG)リピートで構成されています。しかし、c- (TCAGGG)、g- (TGAGGG)、j- (TTGGGG)のテロメア可変リピート(TVR)や他の6量体配列のバリエーションも存在する[ref.1,2,3]。テロメア細胞分裂のたびに短くなり [ref.4]、臨界テロメア長に達すると、DNA損傷反応が誘発され、細胞の老化またはアポトーシスを引き起こす [ref.5, 6]。

 細胞分裂数の限定を回避するため、腫瘍はテロメア維持機構(TMM)としてテロメラーゼの活性化[ref.7]またはalternative lengthening of telomeres(ALT)[ref.8](参考)を採用している。テロメラーゼは、染色体末端にt型リピートを付加する酵素である[ref.9]。対照的に、ALTはテロメア領域の組換えに基づいており、不均一な長さのテロメア[ref.8]および配列構成[ref.3、10]を含むいくつかの特徴をもたらす。

 これらのTMMは腫瘍形成に極めて重要であり、がん治療のための貴重な創薬標的となっている[ref.11]。しかし、様々なタイプの腫瘍におけるこれらのメカニズムを正確に同定し、それを阻害するためには、異なるテロメア構造についてのより多くの洞察が必要である。過去数十年の間に、テロメアの長さと ALT の状態を評価するためのいくつかの実験的方法、例えばテロメア qPCR、末端制限断片(TRF)分析、C サークルアッセイなどが確立されてきた [ref.12, 13]。

 ハイスループットシーケンシングの進歩に伴い、テロメア含量を測定する代替方法が登場した。いくつかの研究では、全ゲノムシークエンシング(WGS)データのテロメアリピートを含むショートリード数がテロメア含量を推定するために使用できることが示され、確立された実験的方法と同等の結果が得られている[ref.10, 14, 15, 16, 17, 18]。この種の解析は、最近発表されたいくつかのがん研究[19,20,21]で説明されているように、がんデータにおけるテロメアの特徴についての貴重な洞察をもたらす。ここでは、一致した腫瘍とコントロールペアのために特別に設計されたテロメア含有量を決定するための新しい計算ツール、TelomereHunterを紹介する。既存のツールとは対照的に、TelomereHunterはアラインメント情報を考慮に入れ、テロメア配列中のバリアントリピートの豊富さを報告する。この論文では、TelomereHunterの主な機能を紹介し、ALT陽性およびALT陰性の腫瘍サンプルにおける例示的な結果の解釈について議論し、テロメア含有量推定のための生物学的アッセイと比較してツールを特徴づけ、異なるシーケンスプロトコルテロメア含有量の定量化に与える影響を評価する。

 

HP

https://www.dkfz.de/en/applied-bioinformatics/telomerehunter/telomerehunter.html

 

インストール

ubuntu18.04のpython2.7環境でテストした(docker使用、ホストmacos10.14)。

依存

operating system: Linux

for telomere read extraction and calculation of telomere content

  • python 2.7.9 (does not work for python 3!)
  • pysam 0.9.0
  • PyPDF2 1.26.0
  • samtools 1.3.1

for visualization

  • R 3.3.0
  • ggplot2 2.1.0
  • reshape2 1.4.1
  • gridExtra 2.2.1
  • RColorBrewer 1.1-2
  • cowplot 0.9.2
  • svglite 1.2.1
#pip (PyPI) 視覚化のためのRのパッケージは別途導入する必要あり。
pip install telomerehunter

> telomerehunter -h

$ telomerehunter -h

 

 

TelomereHunter  Copyright (C) 2015  Lina Sieverling, Philip Ginsbach, Lars Feuerbach

This program comes with ABSOLUTELY NO WARRANTY.

This is free software, and you are welcome to redistribute it

under certain conditions. For details see the GNU General Public License

in the license copy received with TelomereHunter or <http://www.gnu.org/licenses/>.

 

 

TelomereHunter 1.1.0

 

 

usage: telomerehunter [-h] [-ibt TUMOR_BAM] [-ibc CONTROL_BAM] -o OUTPUT_DIR

                      -p PID [-b BANDING_FILE] [-rt REPEAT_THRESHOLD_SET]

                      [-rl] [-mqt MAPQ_THRESHOLD] [-d]

                      [-r REPEATS [REPEATS ...]] [-con] [-gc1 LOWERGC]

                      [-gc2 UPPERGC] [-nf]

                      [-rc TVRS_FOR_CONTEXT [TVRS_FOR_CONTEXT ...]]

                      [-bp BP_CONTEXT] [-pl] [-pff {pdf,png,svg,all}] [-p1]

                      [-p2] [-p3] [-p4] [-p5] [-p6] [-p7] [-p8] [-prc]

 

Estimation of telomere content from WGS data of a tumor and/or a control

sample.

 

optional arguments:

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

  -ibt TUMOR_BAM, --inputBamTumor TUMOR_BAM

                        Path to the indexed input BAM file of the tumor

                        sample.

  -ibc CONTROL_BAM, --inputBamControl CONTROL_BAM

                        Path to the indexed input BAM file of the control

                        sample.

  -o OUTPUT_DIR, --outPath OUTPUT_DIR

                        Path to the output directory into which all results

                        are written.

  -p PID, --pid PID     Sample name used in output files and diagrams

                        (required).

  -b BANDING_FILE, --bandingFile BANDING_FILE

                        Path to a tab-separated file with information on

                        chromosome banding. The first four columns of the

                        table have to contain the chromosome name, the start

                        and end position and the band name. The table should

                        not have a header. If no banding file is specified,

                        the banding information of hg19 will be used.

  -rt REPEAT_THRESHOLD_SET, --repeatThreshold REPEAT_THRESHOLD_SET

                        The number of repeats needed for a read to be

                        classified as telomeric. If no repeat threshold is

                        defined, TelomereHunter will calculate the

                        repeat_threshold depending on the read length with the

                        following formula: repeat_threshold =

                        floor(read_length * 6/100)

  -rl, --perReadLength  Repeat threshold is set per 100 bp read length. The

                        used repeat threshold will be: floor(read_length *

                        repeat_threshold/100) E.g. Setting -rt 8 -rl means

                        that 8 telomere repeats are required per 100 bp read

                        length. If the read length is 50 bp, the threshold is

                        set to 4.

  -mqt MAPQ_THRESHOLD, --mappingQualityThreshold MAPQ_THRESHOLD

                        The mapping quality needed for a read to be considered

                        as mapped (default = 8).

  -d, --removeDuplicates

                        Reads marked as duplicates in the input bam file(s)

                        are removed in the filtering step.

  -r REPEATS [REPEATS ...], --repeats REPEATS [REPEATS ...]

                        List of telomere repeat types to search for. Reverse

                        complements are automatically generated and do not

                        need to be specified! By default, TelomereHunter

                        searches for t-, g-, c- and j-type repeats (TTAGGG

                        TGAGGG TCAGGG TTGGGG).

  -con, --consecutive   Search for consecutive repeats.

  -gc1 LOWERGC, --lowerGC LOWERGC

                        Lower limit used for GC correction of telomere

                        content. The value must be an integer between 0 and

                        100 (default = 48).

  -gc2 UPPERGC, --upperGC UPPERGC

                        Upper limit used for GC correction of telomere

                        content. The value must be an integer between 0 and

                        100 (default = 52).

  -nf, --noFiltering    If the filtering step of TelomereHunter has already

                        been run previously, skip this step.

  -rc TVRS_FOR_CONTEXT [TVRS_FOR_CONTEXT ...], --repeatsContext TVRS_FOR_CONTEXT [TVRS_FOR_CONTEXT ...]

                        List of telomere variant repeats for which to analyze

                        the sequence context. Reverse complements are

                        automatically generated and do not need to be

                        specified! Counts for these telomere variant repeats

                        (arbitrary and singleton context) will be added to the

                        summary table. Default repeats: TCAGGG TGAGGG TTGGGG

                        TTCGGG TTTGGG ATAGGG CATGGG CTAGGG GTAGGG TAAGGG).

  -bp BP_CONTEXT, --bpContext BP_CONTEXT

                        Number of base pairs on either side of the telomere

                        variant repeat to investigate. Please use a number

                        that is divisible by 6.

  -pl, --parallel       The filtering, sorting and estimating steps of the

                        tumor and control sample are run in parallel. This

                        will speed up the computation time of TelomereHunter.

  -pff {pdf,png,svg,all}, --plotFileFormat {pdf,png,svg,all}

                        File format of output diagrams. Choose from pdf

                        (default), png, svg or all (pdf, png and svg).

  -p1, --plotChr        Make diagrams with telomeric reads mapping to each

                        chromosome. If none of the options p1/p2/p3/p4/p5/p6

                        are chosen, all diagrams will be created.

  -p2, --plotFractions  Make a diagram with telomeric reads in each fraction

                        (intrachromosomal, subtelomeric, junction spanning,

                        intratelomeric). If none of the options

                        p1/p2/p3/p4/p5/p6 are chosen, all diagrams will be

                        created.

  -p3, --plotTelContent

                        Make a diagram with the gc corrected telomere content

                        in the analyzed samples. If none of the options

                        p1/p2/p3/p4/p5/p6 are chosen, all diagrams will be

                        created.

  -p4, --plotGC         Make a diagram with GC content distributions in all

                        reads and in intratelomeric reads. If none of the

                        options p1/p2/p3/p4/p5/p6 are chosen, all diagrams

                        will be created.

  -p5, --plotRepeatFreq

                        Make histograms of the repeat frequencies per

                        intratelomeric read. If none of the options

                        p1/p2/p3/p4/p5/p6 are chosen, all diagrams will be

                        created.

  -p6, --plotTVR        Make plots for telomere variant repeats.

  -p7, --plotSingleton  Make plots for singleton telomere variant repeats.

  -p8, --plotNone       Do not make any diagrams. If none of the options

                        p1/p2/p3/p4/p5/p6/p7/p8 are chosen, all diagrams will

                        be created.

  -prc, --plotRevCompl  Distinguish between forward and reverse complement

                        telomere repeats in diagrams.

 

Contact Lina Sieverling (l.sieverling@dkfz-heidelberg.de) for questions and

support.

 

dockerイメージのビルド

git clone https://github.com/cancerit/telomerehunter-docker.git
cd telomerehunter-docker/
docker build -t telomerehunter .

 

 

実行方法

tumorとcontrol(任意)のbamを指定する。

telomerehunter -p prefix -o outdir -ibt tumor.bam -ibc control.bam
  • -ibt    Path to the indexed input BAM file of the tumor  sample.
  • -ibc   Path to the indexed input BAM file of the control  sample.
  • -o     Path to the output directory into which all results are written.
  • -p     Sample name used in output files and diagrams (required).

 

 

引用

TelomereHunter – in silico estimation of telomere content and composition from cancer genomes

Lars Feuerbach, Lina Sieverling, Katharina I. Deeg, Philip Ginsbach, Barbara Hutter, Ivo Buchhalter, Paul A. Northcott, Sadaf S. Mughal, Priya Chudasama, Hanno Glimm, Claudia Scholl, Peter Lichter, Stefan Fröhling, Stefan M. Pfister, David T. W. Jones, Karsten Rippe, Benedikt Brors
BMC Bioinformatics volume 20, Article number: 272 (2019)

 

 

可変数のタンデムリピート(VNTR)をジェノタイピングする adVNTR

 

 全ゲノムシークエンシングは、臨床パイプラインでメンデルバリアントを同定するために使用されることが多くなってきている。これらのパイプラインでは、より複雑な繰り返し配列のバリアントを無視して、一塩基変異(SNV)や構造変異に焦点を当てている。ここでは、短い(6-100bp)繰り返し単位の不正確なタンデム重複からなる可変数タンデムリピート(VNTR)の遺伝子型決定の問題を検討する。VNTRはヒトゲノムの3%を占め、コード領域に多く存在し、複数のメンデル病の原因となっている。既存のツールではVNTRを含む配列を認識することができるが、全ゲノム配列からVNTRのジェノタイピング(繰り返し単位数と配列変異の判定)を行うことは困難である。本研究では、隠れマルコフモデルを用いて各VNTRをモデル化し、繰り返し単位数をカウントし、配列変異を検出する手法であるadVNTRについて述べる。

 

Documentation

Quick Start — adVNTR 1.1.1 documentation

 

インストール

mac10.14でpython2.7の仮想環境を作ってテストした。

本体 Github

#bioconda (link)
conda creatte -n advntr-env python=2.7
conda activate advntr-env
conda install -c bioconda -y advntr

advntr -h

$ advntr -h

=======================================================

adVNTR 1.3.3: Genopyting tool for VNTRs

=======================================================

Source code: https://github.com/mehrdadbakhtiari/adVNTR

Instructions: http://advntr.readthedocs.io

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

 

usage: advntr <command> [options]

 

Command: genotype find RU counts and mutations in VNTRs

         viewmodel view existing models in database

         addmodel add custom VNTR to the database

         delmodel remove a model from database

 

advntr: error: too few arguments

advntr genotype -h

 advntr genotype -h

usage: advntr genotype [options]

 

Input/output options:

  -a/--alignment_file <file>      alignment file in SAM/BAM/CRAM format

  -r/--reference_filename <file>  path to a FASTA-formatted reference file for CRAM files. It overrides

                                  filename specified in header, which is normally used to find the reference

  -f/--fasta <file>               Fasta file containing raw reads

  -p/--pacbio                     set this flag if input file contains PacBio reads instead of Illumina reads

  -n/--nanopore                   set this flag if input file contains Nanopore MinION reads instead of

                                  Illumina

  -o/--outfile <file>             file to write results. adVNTR writes output to stdout if oufile is not

                                  specified.

  -of/--outfmt <format>           output format. Allowed values are {text, bed} [text]

 

Algorithm options:

  -fs/--frameshift                set this flag to search for frameshifts in VNTR instead of copy number.

                                  Supported VNTR IDs: [25561, 519759]

  -e/--expansion                  set this flag to determine long expansion from PCR-free data

  -c/--coverage <float>           average sequencing coverage in PCR-free sequencing

  --haploid                       set this flag if the organism is haploid

  -naive/--naive                  use naive approach for PacBio reads

 

Other options:

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

  --working_directory <path>      working directory for creating temporary files needed for computation

  -m/--models <file>              VNTR models file [vntr_data/hg19_selected_VNTRs_Illumina.db]

  -t/--threads <int>              number of threads [1]

  -u/--update                     set this flag to iteratively update the model

  -vid/--vntr_id <text>           comma-separated list of VNTR IDs

 

 

テストラン

ランにはトレーニング済みのモデルか、ユーザーがトレーニングしたモデルを用意する必要がある。

#pre-trained model
wget https://cseweb.ucsd.edu/~mbakhtia/adVNTR/vntr_data_recommended_loci.zip

> ls -l vntr_data/

$ ls -l vntr_data/ 

total 44280

-rw-r--r--  1 kazuma  staff   9207808  3 22  2019 hg19_selected_VNTRs_Illumina.db

-rw-r--r--  1 kazuma  staff  13463552  3 22  2019 hg19_selected_VNTRs_Pacbio.db

#bam
https://cseweb.ucsd.edu/~mbakhtia/adVNTR/quickstart/

 

bamを指定して実行する。解凍したモデルのvntr_data/は自動で認識される。

mkdir log_dir
advntr genotype --vntr_id 301645 --alignment_file CSTB_2_5_testdata.bam --working_directory ./log_dir/

引用

Targeted genotyping of variable number tandem repeats with adVNTR
Mehrdad Bakhtiari, Sharona Shleizer-Burko, Melissa Gymrek, Vikas Bansal, Vineet Bafna

Genome Research, 28(11), pp.1709-1719

 

 

オルガネラゲノムをターゲットアセンブリする NOVOPlasty

 

 次世代シークエンシング(NGS)技術の進化により、様々なアセンブルアルゴリズムが開発されてきたが、オルガネラゲノムのアセンブルに焦点を当てたものはほとんどない。これらのゲノムは、系統研究や食品の同定に利用されており、GenBankに登録されている真核生物ゲノムの中では最も多い。全ゲノムシークエンシング(WGS)データからオルガネラゲノムのアセンブルを行うことが最も正確で手間のかからない方法であるが、このタスクのために特別に設計されたツールがない(論文執筆時点)。我々(著者ら)は、全ゲノムシークエンシング(WGS)データからオルガネラゲノムをアセンブルするシードアンドエクステンドアルゴリズムを開発した。このアルゴリズムは、いくつかの新しいデータ(Gonioctena intermediaAvicennia marina)と公開データ(Arabidopsis thalianaOryza sativa)の全ゲノムIlluminaデータセットでテストされており、アセンブリの精度とカバレッジにおいて既知のアセンブラよりも優れている。ベンチマークでは、NOVOPlastyは30分以内にテストされた全ての環状ゲノムをアセンブルし、最大16GBのメモリを必要とし、99.99%以上の精度を達成した。結論として、NOVOPlastyは、WGSデータから核外ゲノムを高速かつ簡単に抽出できる唯一のデノボアセンブラである。このソフトウェアはオープンソースで、https://github.com/ndierckx/NOVOPlasty からダウンロードできる。

 

  NOVOPlastyはSSAKEやVCAKEのようなストリングオーバーラップアルゴリズムに似たシードエクステンドベースのアセンブラである。配列をハッシュ・テーブルに格納することから始まる。アセンブリはシードによって開始されなければならず、これは反復的に双方向に拡張される。このシード配列は、アセンブリを開始するために使用されるのではなく、NGSデータセットから目的のゲノムの1つの配列リードを取得するために使用される。この戦略は、アセンブリにミスマッチを組み込むことなく、より広い範囲のシード入力を扱うことができる。シードの配列は、1つのシーケンシングリード、保存された遺伝子、または遠い種からの完全なオルガネラゲノムであってもよい。シードの終わりと始まりは、ハッシュテーブルで重複するリードがないかスキャンされ、別々に保存される。推定されるすべての拡張が同定され、その後、それらが正しく配置されているかどうかを確認するために、ペアになっているリードとクロスチェックされる。比較的類似した配列はグループ化され、すべての塩基拡張はオーバーラップしたリード間のコンセンサスによって解決される。複数のコンセンサス拡張が考えられる場合(すなわち、十分なサイズのグループが1つ以上ある場合)、アセンブリは分割され、2つの新しいコンティグが作成される。ほとんどのアセンブラとは異なり、NOVOPlastyはすべてのリードをアセンブルしようとはせず、環状ゲノムが形成されるまで、与えられたシードを拡張する。アセンブルは、長さが予想される範囲内にあり、両端が少なくとも200bp重なっている場合に環状化する。反復領域が検出された場合は、アセンブリが反復領域を抜けるまで環状化は延期される。全ゲノムデータは通常、核外配列のカバレッジが高いため、このアルゴリズムは、1つのリードを完全な環状ゲノムに拡張することができる(論文1の図2)。
 アセンブルへの新しい線形アプローチに加えて、NOVOPlastyはケースベースの調整を組み込むことで、より高品質なアセンブルを実現する。シーケンシングエラーやゲノムエレメントの混入による問題領域は自動的に検出され、パラメータを調整し、適切な戦略を開始することで可能な限り解決される。現在のIlluminaシーケンシング技術(Illumina HiSeqおよびMiSeq)は、例えば、長い一塩基リピート(SNR)ストレッチの後に非常に高いエラー率を持っており、それはその後のシーケンスを信頼性の低いものにしている(ref.16, pubmed)。これらのエラーが発生しやすい領域は、コンティグアセンブリの中断を引き起こす可能性があり、 この線形アセンブリ戦略では特に問題となる。継続性を確保するためのNOVOPlastyの重要な戦略は、これらのSNRストレッチを早期に検出し、コンセンサスを構築する前に、最も誤ったリードを廃棄することにある。SNRストレッチの正確な長さを定義することは、個々のリード間のSNRの長さの強いばらつきと、全体的に低い品質スコアのために、簡単ではない。コンセンサスが解決できない場合、NOVOPlastyはペアエンド情報を使用して問題のある領域を「ジャンプ」し、そこからアセンブリを再開する。SNR領域は両方向からアプローチされるので、リードの非常に誤った部分を回避することができる。ショートリードアセンブリの他の大きな問題の一つは、複雑な反復配列である。カブトムシのミトコンドリアゲノムは多くの場合、長い高度に反復性の高い部分を含んでいるが(ref.17)、葉緑体ゲノムはより短く分散した反復性の高いDNAを含んでいることがある(ref.18)。NOVOPlastyは繰り返し領域に遭遇した場合、繰り返し配列と上流に隣接する領域を決定する。反復領域の前の配列から始まるすべてのリードは、さらなる分析のためにフィルタリングされる。反復領域の長さがリードの長さよりも短い場合には、領域のアセンブリを直接解決することができる。そうでなければ、アルゴリズムは、リファレンスのポイントとして機能し得る反復配列間の小さな変動を検索するために、非常に厳しくパラメータを一時的に調整する。領域が解決できなかった場合、アセンブリは終了し、繰り返し領域に続く配列で新しいコンティグとしてリブートされる。

 

wiki

https://github.com/ndierckx/NOVOPlasty/wiki

 

インストール

mac10.14でテストした。

本体 Github

git clone https://github.com/ndierckx/NOVOPlasty.git
perl NOVOPlasty/NOVOPlasty3.8.3.pl

> perl NOVOPlasty3.8.3.pl

$ perl NOVOPlasty3.8.3.pl 

 

 

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

NOVOPlasty: The Organelle Assembler

Version 3.8.3

Author: Nicolas Dierckxsens, (c) 2015-2019

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

 

Error:Can't open the configuration file, please check the manual!

 

Usage: perl NOVOPlasty3.8.3.pl -c config.txt

 

実行方法

 ランにはconfigファイルを使用する。fastqのパスやリファレンスのパスを記載する。データはraw fastqかgz/bz2圧縮fastqbになっている必要がある。

git clone https://github.com/ndierckx/NOVOPlasty.git
cp NOVOPlasty/config.txt .

cat config.txt

$ cat config.txt 

Project:

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

Project name          = Test

Type                  = mito

Genome Range          = 12000-22000

K-mer                 = 39

Max memory            = 

Extended log          = 0

Save assembled reads  = no

Seed Input            = /path/to/seed_file/Seed.fasta

Extend seed directly  = no

Reference sequence    = /path/to/reference_file/reference.fasta (optional)

Variance detection    = 

Chloroplast sequence  = /path/to/chloroplast_file/chloroplast.fasta (only for "mito_plant" option)

 

Dataset 1:

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

Read Length           = 151

Insert size           = 300

Platform              = illumina

Single/Paired         = PE

Combined reads        = 

Forward reads         = /path/to/reads/reads_1.fastq

Reverse reads         = /path/to/reads/reads_2.fastq

 

Heteroplasmy:

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

MAF                   = 

HP exclude list       = 

PCR-free              = 

 

Optional:

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

Insert size auto      = yes

Insert Range          = 1.9

Insert Range strict   = 1.3

Use Quality Scores    = no

 

 

 

 

Project:

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

Project name         = Choose a name for your project, it will be used for the output files.

Type                 = (chloro/mito/mito_plant) "chloro" for chloroplast assembly, "mito" for mitochondrial assembly and 

                       "mito_plant" for mitochondrial assembly in plants.

Genome Range         = (minimum genome size-maximum genome size) The expected genome size range of the genome.

                       Default value for mito: 12000-20000 / Default value for chloro: 120000-200000

                       If the expected size is know, you can lower the range, this can be useful when there is a repetitive

                       region, what could lead to a premature circularization of the genome.

K-mer                = (integer) This is the length of the overlap between matching reads (Default: 33). 

                       If reads are shorter then 90 bp or you have low coverage data, this value should be decreased down to 23. 

                       For reads longer then 101 bp, this value can be increased, but this is not necessary.

Max memory           = You can choose a max memory usage, suitable to automatically subsample the data or when you have limited                      

                       memory capacity. If you have sufficient memory, leave it blank, else write your available memory in GB

                       (if you have for example a 8 GB RAM laptop, put down 7 or 7.5 (don't add the unit in the config file))

Extended log         = Prints out a very extensive log, could be useful to send me when there is a problem  (0/1).

Save assembled reads = All the reads used for the assembly will be stored in seperate files (yes/no)

Seed Input           = The path to the file that contains the seed sequence.

Extend seed directly = This gives the option to extend the seed directly, in stead of finding matching reads. Only use this when your seed 

                       originates from the same sample and there are no possible mismatches (yes/no)

Reference (optional) = If a reference is available, you can give here the path to the fasta file.

                       The assembly will still be de novo, but references of the same genus can be used as a guide to resolve 

                       duplicated regions in the plant mitochondria or the inverted repeat in the chloroplast. 

                       References from different genus haven't beeen tested yet.

Variance detection   = If you select yes, you should also have a reference sequence (previous line). It will create a vcf file                

                       with all the variances compared to the give reference (yes/no)

Chloroplast sequence = The path to the file that contains the chloroplast sequence (Only for mito_plant mode).

                       You have to assemble the chloroplast before you assemble the mitochondria of plants!

 

Dataset 1:

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

Read Length          = The read length of your reads.

Insert size          = Total insert size of your paired end reads, it doesn't have to be accurate but should be close enough.

Platform             = illumina/ion - The performance on Ion Torrent data is significantly lower

Single/Paired        = For the moment only paired end reads are supported.

Combined reads       = The path to the file that contains the combined reads (forward and reverse in 1 file)

Forward reads        = The path to the file that contains the forward reads (not necessary when there is a merged file)

Reverse reads        = The path to the file that contains the reverse reads (not necessary when there is a merged file)

 

Heteroplasmy:

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

MAF                  = (0.007-0.49) Minor Allele Frequency: If you want to detect heteroplasmy, first assemble the genome without this option. Then give the resulting                         

                       sequence as a reference and as a seed input. And give the minimum minor allele frequency for this option 

                       (0.01 will detect heteroplasmy of >1%)

HP exclude list      = Option not yet available  

PCR-free             = (yes/no) If you have a PCR-free library write yes

 

Optional:

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

Insert size auto     = (yes/no) This will finetune your insert size automatically (Default: yes)

Insert Range         = This variation on the insert size, could lower it when the coverage is very high or raise it when the

                       coverage is too low (Default: 1.9). 

Insert Range strict  = Strict variation to resolve repetitive regions (Default: 1.3).                                

Use Quality Scores   = It will take in account the quality scores, only use this when reads have low quality, like with the    

                       300 bp reads of Illumina (yes/no)

 

 configファイルを指定して実行する。

perl NOVOPlasty3.2.pl -c config.txt

 

引用

NOVOPlasty: de novo assembly of organelle genomes from whole genome data

Nicolas Dierckxsens, Patrick Mardulyn, Guillaume Smits
Nucleic Acids Research, Volume 45, Issue 4, 28 February 2017, Page e18

 

Unraveling heteroplasmy patterns with NOVOPlasty

Nicolas Dierckxsens, Patrick Mardulyn, Guillaume Smits
NAR Genomics and Bioinformatics, Volume 2, Issue 1, March 2020, lqz011

 

 参考


anvi'oを使ってパンゲノム解析を行う

 

anvi'o パンゲノミックワークフローは、3 つの主要なステップで構成されている。

1、anvi-gen-genomes-storage による anvi'o ゲノムストレージ生成

2、anvi-pan-genomeによるパンデータベース生成

3、anvi-display-pan(ゲノムストレージと パンデータベースが必要)を使用して結果を表示

そして、インタラクティブなインターフェースを使用して遺伝子クラスタをコレクションにまとめたりサマリーレポートを作成したりできる。チュートリアルに沿って使い方を見ていく。

 

インストール

公式dockerイメージを使ってubuntu18.04LTSマシンでテストした。

本体 Github

#依存が多いので、condaだと依存チェックに異常な時間がかかる。動かすだけならdockerが楽。

#docker (dockerhub) (link)

#latest (v6)
docker pull meren/anvio:latest

インストールチェック

> anvi-self-test --suite pangenomics

起動

docker run --rm -it -v `pwd`:`pwd` -w `pwd` -p 8080:8080 meren/anvio:latest

anvi-migrate -h

> anvi-migrate -h

usage: anvi-migrate [-h] [--just-do-it] [-t VERSION] DATABASE [DATABASE ...]

 

Migrate an anvi'o database or config file to a newer version.

 

positional arguments:

  DATABASE              Anvi'o database or config file for migration

 

optional arguments:

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

  --just-do-it          Do not bother me with warnings

  -t VERSION, --target-version VERSION

                        Anvi'o will stop upgrading your database when it

                        reaches to this version.

 :: anvi'o v6.2 :: 

anvi-gen-genomes-storage -h

> anvi-gen-genomes-storage -h

usage: anvi-gen-genomes-storage [-h] [-e FILE_PATH] [-i FILE_PATH]

                                [--gene-caller GENE-CALLER] -o GENOMES_STORAGE

 

Create a genome storage from internal or external genomes for a pan genome

analysis.

 

optional arguments:

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

 

EXTERNAL GENOMES:

  External genomes listed as anvi'o contigs databases. As in, you have one

  or more genomes say from NCBI you want to work with, and you created an

  anvi'o contigs database for each one of them.

 

  -e FILE_PATH, --external-genomes FILE_PATH

                        A two-column TAB-delimited flat text file that lists

                        anvi'o contigs databases. The first item in the header

                        line should read 'name', and the second should read

                        'contigs_db_path'. Each line in the file should

                        describe a single entry, where the first column is the

                        name of the genome (or MAG), and the second column is

                        the anvi'o contigs database generated for this genome.

 

INTERNAL GENOMES:

  Genome bins stored in an anvi'o profile databases as collections.

 

  -i FILE_PATH, --internal-genomes FILE_PATH

                        A five-column TAB-delimited flat text file. The header

                        line must contain these columns: 'name', 'bin_id',

                        'collection_id', 'profile_db_path', 'contigs_db_path'.

                        Each line should list a single entry, where 'name' can

                        be any name to describe the anvi'o bin identified as

                        'bin_id' that is stored in a collection.

 

PRO STUFF:

  Things you may not have to change. But you never know (unless you read the

  help).

 

  --gene-caller GENE-CALLER

                        The gene caller to utilize. Anvi'o supports multiple

                        gene callers, and some operations (including this one)

                        requires an explicit mentioning of which one to use.

                        The default is 'prodigal', but it will not be enough

                        if you if you were a rebel and have used `--external-

                        gene-callers` or something.

 

OUTPUT:

  Give it a nice name. Must end with '-GENOMES.db'. This is primarily due to

  the fact that there are other .db files used throughout anvi'o and it

  would be better to distinguish this very special file from them.

 

  -o GENOMES_STORAGE, --output-file GENOMES_STORAGE

                        File path to store results.

 :: anvi'o v6.2 :: 

anvi-pan-genome -h

> anvi-pan-genome -h

 

WARNING

===============================================

If you publish results from this workflow, please do not forget to cite DIAMOND

(doi:10.1038/nmeth.3176), unless you use it with --use-ncbi-blast flag, and MCL

(http://micans.org/mcl/ and doi:10.1007/978-1-61779-361-5_15)

 

usage: anvi-pan-genome [-h] -g GENOMES_STORAGE [-G GENOME_NAMES]

                       [--skip-alignments] [--skip-homogeneity]

                       [--quick-homogeneity] [--align-with ALIGNER]

                       [--exclude-partial-gene-calls] [--use-ncbi-blast]

                       [--minbit MINBIT] [--mcl-inflation INFLATION]

                       [--min-occurrence NUM_OCCURRENCE]

                       [--min-percent-identity PERCENT] [--sensitive]

                       [-n PROJECT_NAME] [--description TEXT_FILE]

                       [-o PAN_DB_DIR] [-W] [-T NUM_THREADS]

                       [--skip-hierarchical-clustering]

                       [--enforce-hierarchical-clustering]

                       [--distance DISTANCE_METRIC] [--linkage LINKAGE_METHOD]

 

A DIAMOND and MCL-based anvi'o workflow for pangenomics. You provide genomes

from anywhere (whether they are external genomes, or anvi'o genome bins in

collections), and it gives you back a pangenome analysis.

 

optional arguments:

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

 

GENOMES:

  The very fancy genomes storage file. This file is generated by the program

  `anvi-genomes-storage`. Please see the online tutorial on pangenomic

  workflow if you don't know how to generate one.

 

  -g GENOMES_STORAGE, --genomes-storage GENOMES_STORAGE

                        Anvi'o genomes storage file

  -G GENOME_NAMES, --genome-names GENOME_NAMES

                        Genome names to 'focus'. You can use this parameter to

                        limit the genomes included in your analysis. You can

                        provide these names as a comma-separated list of

                        names, or you can put them in a file, where you have a

                        single genome name in each line, and provide the file

                        path.

 

PARAMETERS:

  Important stuff Tom never pays attention (but you should).

 

  --skip-alignments     By default, anvi'o attempts to align amino acid

                        sequences in each gene cluster using multiple sequnce

                        alignment via muscle. You can use this flag to skip

                        that step and be upset later.

  --skip-homogeneity    By default, anvi'o attempts to calculate homogeneity

                        values for every gene cluster, given that they are

                        aligned. You can use this flag to have anvi'o skip

                        homogeneity calculations. Anvi'o will ignore this flag

                        if you decide to skip alignments

  --quick-homogeneity   By default, anvi'o will use a homogeneity algorithm

                        that checks for horizontal and vertical geometric

                        homogeneity (along with functional). With this flag,

                        you can tell anvi'o to skip horizontal geometric

                        homogeneity calculations. It will be less accurate but

                        quicker. Anvi'o will ignore this flag if you skip

                        homogeneity calculations or alignments all together.

  --align-with ALIGNER  The multiple sequence alignment program to use when

                        multiple sequence alignment is necessary. To see all

                        available options, use the flag `--list-aligners`.

  --exclude-partial-gene-calls

                        By default, anvi'o includes all partial gene calls

                        from the analysis, which, in some cases, may inflate

                        the number of gene clusters identified and introduce

                        extra heterogeneity within those gene clusters. Using

                        this flag, you can request anvi'o to exclude partial

                        gene calls from the analysis (whether a gene call is

                        partial or not is an information that comes directly

                        from the gene caller used to identify genes during the

                        generation of the contigs database).

  --use-ncbi-blast      This program uses DIAMOND by default, however, if you

                        like, you can use good ol' blastp from NCBI instead.

  --minbit MINBIT       The minimum minbit value. The minbit heuristic

                        provides a mean to set a to eliminate weak matches

                        between two amino acid sequences. We learned it from

                        ITEP (Benedict MN et al, doi:10.1186/1471-2164-15-8),

                        which is a comprehensive analysis workflow for

                        pangenomes, and decided to use it in the anvi'o

                        pangenomic workflow, as well. Briefly, If you have two

                        amino acid sequences, 'A' and 'B', the minbit is

                        defined as 'BITSCORE(A, B) / MIN(BITSCORE(A, A),

                        BITSCORE(B, B))'. So the minbit score between two

                        sequences goes to 1 if they are very similar over the

                        entire length of the 'shorter' amino acid sequence,

                        and goes to 0 if (1) they match over a very short

                        stretch compared even to the length of the shorter

                        amino acid sequence or (2) the match betwen sequence

                        identity is low. The default is 0.5.

  --mcl-inflation INFLATION

                        MCL inflation parameter, that defines the sensitivity

                        of the algorithm during the identification of the gene

                        clusters. More information on this parameter and it's

                        effect on cluster granularity is here:

                        (http://micans.org/mcl/man/mclfaq.html#faq7.2). The

                        default is 2.

  --min-occurrence NUM_OCCURRENCE

                        Do you not want singletons?\ You don't? Well, this

                        parameter will help you get rid of them (along with

                        doubletons, if you want). Anvi'o will remove gene

                        clusters that occur less than the number you set using

                        this parameter from the analysis. The default is 1,

                        which means everything will be kept. If you want to

                        remove singletons, set it to 2, if you want to remove

                        doubletons as well, set it to 3, and so on.

  --min-percent-identity PERCENT

                        Minimum percent identity between the two amino acid

                        sequences for them to have an edge for MCL analysis.

                        This value will be used to filter hits from Diamond

                        search results. Because percent identity is not a

                        predictor of a good match (since it does not

                        communicate many other important factors such as the

                        alignment length between the two sequences and its

                        proportion to the entire length of those involved), we

                        suggest you rely on 'minbit' parameter. But you know

                        what? Maybe you shouldn't listen to anyone, and

                        experiment on your own! The default is 0 percent.

  --sensitive           DIAMOND sensitivity. With this flag you can instruct

                        DIAMOND to be 'sensitive', rather than 'fast' during

                        the search. It is likely the search will take

                        remarkably longer. But, hey, if you are doing it for

                        your final analysis, maybe it should take longer and

                        be more accurate. This flag is only relevant if you

                        are running DIAMOND.

 

OTHERS:

  Sweet parameters of convenience.

 

  -n PROJECT_NAME, --project-name PROJECT_NAME

                        Name of the project. Please choose a short but

                        descriptive name (so anvi'o can use it whenever she

                        needs to name an output file, or add a new table in a

                        database, or name her first born).

  --description TEXT_FILE

                        A plain text file that contains some description about

                        the project. You can use Markdwon syntax. The

                        description text will be rendered and shown in all

                        relevant interfaces, including the anvi'o interactive

                        interface, or anvi'o summary outputs.

  -o PAN_DB_DIR, --output-dir PAN_DB_DIR

                        Directory path for output files

  -W, --overwrite-output-destinations

                        Overwrite if the output files and/or directories

                        exist.

  -T NUM_THREADS, --num-threads NUM_THREADS

                        Maximum number of threads to use for multithreading

                        whenever possible. Very conservatively, the default is

                        1. It is a good idea to not exceed the number of CPUs

                        / cores on your system. Plus, please be careful with

                        this option if you are running your commands on a SGE

                        --if you are clusterizing your runs, and asking for

                        multiple threads to use, you may deplete your

                        resources very fast.

 

ORGANIZING GENE CLUSTERs:

  These are stuff that will change the clustering dendrogram of your gene

  clusters.

 

  --skip-hierarchical-clustering

                        Anvi'o attempts to generate a hierarchical clustering

                        of your gene clusters once it identifies them so you

                        can use `anvi-display-pan` to play with it. But if you

                        want to skip this step, this is your flag.

  --enforce-hierarchical-clustering

                        If you want anvi'o to try to generate a hierarchical

                        clustering of your gene clusters even if the number of

                        gene clusters exceeds its suggested limit for

                        hierarchical clustering, you can use this flag to

                        enforce it. Are you are a rebel of some sorts? Or did

                        computers made you upset? Express your anger towards

                        machine using this flag.

  --distance DISTANCE_METRIC

                        The distance metric for the clustering of gene

                        clusters. If you do not use this flag, the default

                        distance metric will be used for each clustering

                        configuration which is "euclidean".

  --linkage LINKAGE_METHOD

                        The same story with the `--distance`, except, the

                        system default for this one is ward.

 :: anvi'o v6.2 :: 

>anvi-display-pan -h

> anvi-display-pan -h 

usage: anvi-display-pan [-h] -p PAN_DB [-g GENOMES_STORAGE] [-d VIEW_DATA]                                                                                                                                     

                        [-t NEWICK] [-V ADDITIONAL_VIEW]

                        [-A ADDITIONAL_LAYERS] [--view NAME] [--title NAME]

                        [--state-autoload NAME] [--collection-autoload NAME]

                        [--export-svg FILE_PATH] [--skip-init-functions]

                        [--dry-run] [--skip-auto-ordering] [-I IP_ADDR]

                        [-P INT] [--browser-path PATH] [--read-only]

                        [--server-only] [--password-protected]

                        [--user-server-shutdown]

 

Start an anvi'o server to display a pan-genome

 

optional arguments:

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

 

INPUT FILES:

  Input files from the pangenome analysis.

 

  -p PAN_DB, --pan-db PAN_DB

                        Anvi'o pan database

  -g GENOMES_STORAGE, --genomes-storage GENOMES_STORAGE

                        Anvi'o genomes storage file

 

OPTIONAL INPUTS:

  Where the yay factor becomes a reality.

 

  -d VIEW_DATA, --view-data VIEW_DATA

                        A TAB-delimited file for view data

  -t NEWICK, --tree NEWICK

                        NEWICK formatted tree structure

 

ADDITIONAL STUFF:

  Parameters to provide additional layers, views, or layer data.

 

  -V ADDITIONAL_VIEW, --additional-view ADDITIONAL_VIEW

                        A TAB-delimited file for an additional view to be used

                        in the interface. This file should contain all split

                        names, and values for each of them in all samples.

                        Each column in this file must correspond to a sample

                        name. Content of this file will be called 'user_view',

                        which will be available as a new item in the 'views'

                        combo box in the interface

  -A ADDITIONAL_LAYERS, --additional-layers ADDITIONAL_LAYERS

                        A TAB-delimited file for additional layers for splits.

                        The first column of this file must be split names, and

                        the remaining columns should be unique attributes. The

                        file does not need to contain all split names, or

                        values for each split in every column. Anvi'o will try

                        to deal with missing data nicely. Each column in this

                        file will be visualized as a new layer in the tree.

 

VISUALS RELATED:

  Parameters that give access to various adjustements regarding the

  interface.

 

  --view NAME           Start the interface with a pre-selected view. To see a

                        list of available views, use --show-views flag.

  --title NAME          Title for the interface. If you are working with a

                        RUNINFO dict, the title will be determined based on

                        information stored in that file. Regardless, you can

                        override that value using this parameter.

  --state-autoload NAME

                        Automatically load previous saved state and draw tree.

                        To see a list of available states, use --show-states

                        flag.

  --collection-autoload NAME

                        Automatically load a collection and draw tree. To see

                        a list of available collections, use --list-

                        collections flag.

  --export-svg FILE_PATH

                        The SVG output file path.

 

SWEET PARAMS OF CONVENIENCE:

  Parameters and flags that are not quite essential (but nice to have).

 

  --skip-init-functions

                        When declared, function calls for genes will not be

                        initialized (therefore will be missing from all

                        relevant interfaces or output files). The use of this

                        flag may reduce the memory fingerprint and processing

                        time for large datasets.

  --dry-run             Don't do anything real. Test everything, and stop

                        right before wherever the developer said 'well, this

                        is enough testing', and decided to print out results.

  --skip-auto-ordering  When declared, the attempt to include automatically

                        generated orders of items based on additional data is

                        skipped. In case those buggers cause issues with your

                        data, and you still want to see your stuff and deal

                        with the other issue maybe later.

 

SERVER CONFIGURATION:

  For power users.

 

  -I IP_ADDR, --ip-address IP_ADDR

                        IP address for the HTTP server. The default ip address

                        (0.0.0.0) should work just fine for most.

  -P INT, --port-number INT

                        Port number to use for anvi'o services. If nothing is

                        declared, anvi'o will try to find a suitable port

                        number, starting from the default port number, 8080.

  --browser-path PATH   By default, anvi'o will use your default browser to

                        launch the interactive interface. If you would like to

                        use something else than your system default, you can

                        provide a full path for an alternative browser using

                        this parameter, and hope for the best. For instance we

                        are using this parameter to call Google's experimental

                        browser, Canary, which performs better with demanding

                        visualizations.

  --read-only           When the interactive interface is started with this

                        flag, all 'database write' operations will be

                        disabled.

  --server-only         The default behavior is to start the local server, and

                        fire up a browser that connects to the server. If you

                        have other plans, and want to start the server without

                        calling the browser, this is the flag you need.

  --password-protected  If this flag is set, command line tool will ask you to

                        enter a password and interactive interface will be

                        only accessible after entering same password. This

                        option is recommended for shared machines like

                        clusters or shared networks where computers are not

                        isolated.

  --user-server-shutdown

                        Allow users to shutdown an anvi'server via web

                        interface.

 :: anvi'o v6.2 :: 

anvi-split -h

> anvi-split -h

usage: anvi-split [-h] -p PAN_OR_PROFILE_DB [-c CONTIGS_DB]

                  [-g GENOMES_STORAGE] [--skip-variability-tables]

                  [--compress-auxiliary-data] [-C COLLECTION_NAME]

                  [-b BIN_NAME] [-o DIR_PATH] [--list-collections]

                  [--skip-hierarchical-clustering]

                  [--enforce-hierarchical-clustering]

                  [--distance DISTANCE_METRIC] [--linkage LINKAGE_METHOD]

 

Split an anvi'o pan or profile database into smaller, self-contained pieces.

This is usually great when you want to share a subset of an anvi'o project.

You give this guy your databases, and a collection id, and it gives you back

directories of individual projects for each bin that can be treated as self-

contained smaller anvi'o projects. We know you don't read this far into these

help menus, but please remember: you will either need to provide a profile &

contigs database pair, or a pan & genomes storage pair. The rest will be taken

care of. Magic.

 

optional arguments:

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

 

DATABASES:

  You will either provide a PROFILE/CONTIGS or a PAN/GENOMES STORAGE pair

  here.

 

  -p PAN_OR_PROFILE_DB, --pan-or-profile-db PAN_OR_PROFILE_DB

                        Anvi'o pan or profile database (and even genes

                        database in appropriate contexts).

  -c CONTIGS_DB, --contigs-db CONTIGS_DB

                        Anvi'o contigs database generated by 'anvi-gen-

                        contigs'

  -g GENOMES_STORAGE, --genomes-storage GENOMES_STORAGE

                        Anvi'o genomes storage file

 

PROFILE/CONTIGS OPTIONS:

  Some options that are specific to this only.

 

  --skip-variability-tables

                        Processing variability tables in profile database

                        might take a very long time. With this flag you will

                        be asking anvi'o to skip them.

  --compress-auxiliary-data

                        When declared, the auxiliary data file in the

                        resulting output will be compressed. This saves space,

                        but it takes long. Also, if you are planning to

                        compress the entire later using GZIP, it is even

                        useless to do. But you are the boss!

 

COLLECTION:

  You should provide a valid collection name. If you do not provide bin

  names, the program will generate an output for each bin in your collection

  separately.

 

  -C COLLECTION_NAME, --collection-name COLLECTION_NAME

                        Collection name.

  -b BIN_NAME, --bin-id BIN_NAME

                        Bin name you are interested in.

 

OUTPUT:

  Where do we want the resulting split profiles to be stored.

 

  -o DIR_PATH, --output-dir DIR_PATH

                        Directory path for output files

 

EXTRAS:

  Stuff that you rarely need, but you really really need when the time

  comes. Following parameters will aply to each of the resulting anvi'o

  profile that will be split from the mother anvi'o profile.

 

  --list-collections    Show available collections and exit.

  --skip-hierarchical-clustering

                        If you are not planning to use the interactive

                        interface (or if you have other means to add a tree of

                        contigs in the database) you may skip the step where

                        hierarchical clustering of your items are preformed

                        based on default clustering recipes matching to your

                        database type.

  --enforce-hierarchical-clustering

                        If you have more than 25,000 splits in your merged

                        profile, anvi-merge will automatically skip the

                        hierarchical clustering of splits (by setting --skip-

                        hierarchical-clustering flag on). This is due to the

                        fact that computational time required for hierarchical

                        clustering increases exponentially with the number of

                        items being clustered. Based on our experience we

                        decided that 25,000 splits is about the maximum we

                        should try. However, this is not a theoretical limit,

                        and you can overwrite this heuristic by using this

                        flag, which would tell anvi'o to attempt to cluster

                        splits regardless.

  --distance DISTANCE_METRIC

                        The distance metric for the hierarchical clustering.

                        If you do not use this flag, the default distance

                        metric will be used for each clustering configuration

                        which is "euclidean".

  --linkage LINKAGE_METHOD

                        The same story with the `--distance`, except, the

                        system default for this one is ward.

 :: anvi'o v6.2 :: 

 

 

ラン

チュートリアルに従って進める。 


anvi'oのメタゲノムビニングFASTA(internal genome)、ユーザーが用意したFASTA(external genome)を利用できる。

 

1、データのダウンロード。チュートリアルが使っているProchlorococcusのデータになる。

wget https://ndownloader.figshare.com/files/11857577 -O Prochlorococcus_31_genomes.tar.gz
tar -zxvf Prochlorococcus_31_genomes.tar.gz
cd Prochlorococcus_31_genomes

f:id:kazumaxneo:20200404125612p:plain

各ゲノムの.dbファイルがダンロードされる。

 

2、起動と.dbの読み込み

#lauch
docker run --rm -it -v `pwd`:`pwd` -w `pwd` -p 8080:8080 meren/anvio:latest

#h5pyがないと怒られたら導入
pip install h5py

#すでにある.dbを更新
anvi-migrate *.db

更新された.dbファイルがそのままカレントに出力される。

 

3、 全ゲノムの.dbを出力。任意で各ゲノムの追加情報を含むTAB区切りファイル(external-genomes.txt)を指定する。

external-genomes.txtの中身

f:id:kazumaxneo:20200404203008p:plain

テキストを用意したら実行する。

anvi-gen-genomes-storage -e external-genomes.txt -o PROCHLORO-GENOMES.db

f:id:kazumaxneo:20200404130701p:plain

指定したPROCHLORO-GENOMES.dbが出力される。

4、 ゲノムストレージの準備ができたら、anvi-pan-genomeプログラムを使ってパンゲノム解析を実行する。

anvi-pan-genome -g PROCHLORO-GENOMES.db -n PROJECT1 -T 40

ディレクトリ PROJECT1/ができ、中にパンゲノムデータベースPROJECT1-PAN.dbなどが出力される。

 

5、結果を視覚化する。メタゲノムのanvi-interactiveに似たコマンドanvi-display-panを実行する(パンゲノム用の微調整が入っている)。

anvi-display-pan -p PROJECT1/PROJECT1-PAN.db -g PROCHLORO-GENOMES.db 

http://localhost:8080 にアクセスする。

 

Drawで描画

f:id:kazumaxneo:20200404134315p:plain

図は左下のボタンやマウスのホイールから自由に拡大縮小できる。環状ゲノムを表した図ではなく、遺伝子の有無を線で表した図である事に注意する。

 

メニューを隠した。左斜め上の方の、全てのリングが黒くなっている領域がコア遺伝子のクラスタ。左下から右下の方は全てアクセサリ遺伝子。

f:id:kazumaxneo:20200404134409p:plain

遺伝子クラスタリングの結果に基づいてレイヤーの順番(つまり中心から外周までのリングの順番)を並べ換える。Layerタブ=> Order by => gene_cluster frequenciesを選択。

f:id:kazumaxneo:20200404141225p:plain

 

実行すると、クラスタリング結果に基づいてレイヤーの順番が並び替えられ、右上に遺伝子クラスタリング結果のデンドログラムも出現する。

f:id:kazumaxneo:20200404141443p:plain

 

MainタブのItem orderはリング内でのオーダーの指定になる。変更すると、レイヤーの順番は変わらず、1レイヤー内の遺伝子の順番がクラスタリングされる。例えば左斜め上のコア遺伝子クラスタ(真っ黒の部分)が右下に移動したりする。

 

 

 

カスタム設定のプロファイル(JSON)を左下ボタンから設定をexportし、次回起動時にカスタムした設定がデフォルトの状態で視覚化できる。いったんhaltして、以下を実行した。

#自分のprofileをtest1としてexportした場合、--nameをtest1と指定。
anvi-import-state -p PROJECT1/PROJECT1-PAN.db \
--state pan-state.json \
--name default

#再実行。
anvi-display-pan -p PROJECT1/PROJECT1-PAN.db -g PROCHLORO-GENOMES.db

ここではexampleのJSONファイルを読み込ませている。

 

最初から色などがカスタム設定にした状態で描画された。

 

f:id:kazumaxneo:20200404142810p:plain

 

 

例えば左上の共通して保存されている遺伝子群にコア遺伝子と表記をつけたいとする。その場合、まずBinのタブに移動し、

f:id:kazumaxneo:20200404152700p:plain

 

アサインしたい名前をbin_1からcore geneという名前を変更した。色は赤にした。

f:id:kazumaxneo:20200404153050p:plain

 

マウスホイールで少し拡大。中央のデンドログラムのコア遺伝子の枝付近にマウスをホバーすると、下の図のように該当する枝がリアルタイムでハイライト表示されるので。、

f:id:kazumaxneo:20200404153323p:plain

core geneの表記をつけたい枝部分にマウスを合わせる。

 

ホバーされて色が変わった状態で1回左クリックする。その枝の最外周にCore geneという表記がついた。

f:id:kazumaxneo:20200404153650p:plain

(持たないゲノムまでコア遺伝子に含んでいるので間違ってます)

 

左上のメニューもプラスをクリック、名前がaccessory gene(色は青)というタグを作成、残りの枝の最外周にaccessory geneという表記をつけた。

f:id:kazumaxneo:20200404153853p:plain

下にaccessory geneの表記をつけた。

 

大まかな傾向が見られたら、そのうち個別の遺伝子クラスターに興味が出てくる。その場合、関心がある遺伝子クラスターにズームアップし、線の上で右クリックして詳細を調べることができる。

f:id:kazumaxneo:20200404152006p:plain

この図では、一番内側の3レイヤーにしかない線を右クリックしている。右端にウィンドウが出現ている。

右クリック => ウィンドウのinspect gene clusterを選択すると、新しいウィンドウが生成され。そこに下のように配列が表示される。

f:id:kazumaxneo:20200404152221p:plain

線があったゲノム、つまりそのタンパク質が見つかったゲノムだけアミノ酸配列が表示されている。ここでは一番上の3つになる。バッググラウンドの色はリングの色をそのまま反映している。順番も元の通り(上の例だと一番上が円の一番内側)。

 

右クリックからアミノ酸配列をコピーしたりもできる。

f:id:kazumaxneo:20200404150917p:plain

 

 

コア遺伝子、アクセサリ遺伝子などに分けて描画する。いったんhaltして、anvi-splitを実行する。

anvi-display-pan -p PROJECT1/PROJECT1-PAN.db -g PROCHLORO-GENOMES.db

ランが終わったら、 出力ディレクトリのコア遺伝子やアクセサリ遺伝子の.dbを指定して描画する。

引用

Anvi'o: an advanced analysis and visualization platform for 'omics data

Eren AM, Esen ÖC, Quince C, Vineis JH, Morrison HG, Sogin ML, Delmont TO

PeerJ. 2015 Oct 8;3:e1319

 

 関連


作画の参考にした論文

https://www.biorxiv.org/content/10.1101/170639v1.full.pdf

 

パンゲノムグラフから微生物の多様性を調べる PPanGGOLiN

2020 4/10 引用追加、タイトル修正

 

 機能研究、進化研究、疫学研究のために比較ゲノムを使用するには、与えられた種での発現の観点から遺伝子ファミリーを分類する方法が必要である。これらの方法は、通常、分割や最適なクラス数を推論するための多変量統計モデルを持たず、またゲノム構成を考慮していない。本研究では、パンゲノムをモデル化するために、ノードが遺伝子ファミリーを、エッジがゲノム近傍を表すグラフ構造を導入した。この手法はPPanGGOLiNと名付けられ、多変量Bernoulli混合モデルとMarkov Random Fieldを組み合わせた期待値最大化アルゴリズムを用いてノードを分割する。このアプローチでは、グラフのトポロジーとパンゲノム中の遺伝子の有無を考慮して、遺伝子ファミリーを永続的、クラウド、1つまたは複数のシェルパーティションに分類する。本研究では、439種の単離されたゲノムと78種のメタゲノムアセンブリゲノムのパンゲノム分割グラフを解析し、本手法が永続的ゲノムの推定に有効であることを示した。興味深いことに、シェルゲノムはゲノムダイナミクスを理解する上で重要な要素であることが示された。PPanGGOLiNが提案するグラフベースのアプローチは、何千もの株のゲノム全体の多様性をコンパクトな構造で表現するのに有用であり、非常に大規模な比較ゲノムのための有効な基礎を提供する。本ソフトウェアは、https://github.com/labgem/PPanGGOLiN から自由に入手可能である。

 PPanGGOLiNは、ゲノムDNA配列または提供されたゲノムアノテーションのセットから原核生物のパンゲノムを作成し、操作するソフトウェアスイートである。数万のゲノムまでスケールアップできるように設計されている。

PPanGGOLiNは、グラフィカルモデルと統計的手法を用いて、永続ゲノム、シェルゲノム、クラウドゲノムの遺伝子ファミリーを分割し、パンゲノムを構築する。また、その際には、遺伝子ファミリーのグラフを作成するために、タンパク質をコードする遺伝子とそのゲノム近傍の情報を統合し、各ノードが遺伝子ファミリー、各エッジが遺伝的連続性の関係である遺伝子ファミリーのグラフを作成する。この分割法は、グラフ上で隣りあった2つの遺伝子ファミリーが同じ分割に属する可能性が高いことを促進する。これにより、微生物の多様性をナビゲートするために、生物学者が地下鉄地図のようにゲノムをレール上に描く永続ノード、シェルノード、クラウドノードで構成された分割パンゲノムグラフ(PPG)が完成する。

 さらに、PPanGGOLiNに含まれるpanRGP法(Bazin et al. 2020)は、各ゲノムについて、パンゲノムグラフのシェルゲノムとクラウドゲノムからなる遺伝子のクラスターであるRegions of Genome Plasticity(RPP)を予測する。その多くは遺伝子水平伝播(HGT)から発生し、ゲノムアイランド(GI)に対応している。次に、異なるゲノムからのRGPは、保存されている側方の永続的な遺伝子に基づいて、挿入のスポットでグループ化される。

 

wiki

https://github.com/labgem/PPanGGOLiN/wiki

 


 

インストール

ubuntu18.04LTSでテストした。

本体 Github

#bioconda (link)
conda create -n ppanggolin-env -y
conda activate ppanggolin-env
conda install -c bioconda ppanggolin -y

ppanggolin -h

$ ppanggolin -h

usage: ppanggolin [-h] [-v]  ...

 

Depicting microbial species diversity via a Partitioned PanGenome Graph Of Linked Neighbors

 

optional arguments:

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

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

 

subcommands:

  

    Basic:

      workflow      Easy workflow to run a pangenome analysis in one go

      panrgp        Easy workflow to run a pangenome analysis with genomic islands and spots of insertion detection

    

    Expert:

      annotate      Annotate genomes

      cluster       Cluster proteins in protein families

      graph         Create the pangenome graph

      partition     Partition the pangenome graph

      rarefaction   Compute the rarefaction curve of the pangenome

    

    Output:

      draw          Draw figures representing the pangenome through different aspects

      write         Writes 'flat' files representing the pangenome that can be used with other softwares

      info          Prints information about a given pangenome graph file

    

    Regions of genomic Plasticity:

      align         Aligns a genome or a set of proteins to the pangenome gene families representatives and predict informations from it

      rgp           Predicts Regions of Genomic Plasticity in the genomes of your pangenome

      spot          Predicts spots in your pangenome

 

ppanggolin workflow -h

$ ppanggolin workflow -h

usage: ppanggolin workflow [-h] [--fasta FASTA] [--anno ANNO]

                           [--clusters CLUSTERS] [-o OUTPUT]

                           [--basename BASENAME] [--rarefaction]

                           [-K NB_OF_PARTITIONS] [--defrag] [--tmpdir TMPDIR]

                           [--verbose {0,1,2}] [--log LOG] [-c CPU] [-f]

 

Input arguments:

  The possible input arguments :

 

  --fasta FASTA         A tab-separated file listing the organism names, and

                        the fasta filepath of its genomic sequence(s) (the

                        fastas can be compressed). One line per organism. This

                        option can be used alone. (default: None)

  --anno ANNO           A tab-separated file listing the organism names, and

                        the gff filepath of its annotations (the gffs can be

                        compressed). One line per organism. This option can be

                        used alone IF the fasta sequences are in the gff

                        files, otherwise --fasta needs to be used. (default:

                        None)

  --clusters CLUSTERS   a tab-separated file listing the cluster names, the

                        gene IDs, and optionnally whether they are a fragment

                        or not. (default: None)

 

Optional arguments:

  -o OUTPUT, --output OUTPUT

                        Output directory (default:

                        ppanggolin_output_DATE2020-04-10_HOUR14.53.05_PID8487)

  --basename BASENAME   basename for the output file (default: pangenome)

  --rarefaction         Use to compute the rarefaction curves (WARNING: can be

                        time consumming) (default: False)

  -K NB_OF_PARTITIONS, --nb_of_partitions NB_OF_PARTITIONS

                        Number of partitions to use. Must be at least 3. If

                        under 3, it will be detected automatically. (default:

                        -1)

  --defrag              Realign gene families to associated fragments with

                        their non-fragmented gene family. (default: False)

 

Common arguments:

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

  --tmpdir TMPDIR       directory for storing temporary files (default:

                        /var/folders/v0/rdv6fcx11fz1dy2jjvthmx1r0000gn/T)

  --verbose {0,1,2}     Indicate verbose level (0 for warning and errors only,

                        1 for info, 2 for debug) (default: 1)

  --log LOG             log output file (default: stdout)

  -c CPU, --cpu CPU     Number of available cpus (default: 1)

  -f, --force           Force writing in output directory and in pangenome

                        output file. (default: False)

ppanggolin panrgp -h

$ ppanggolin panrgp -h

usage: ppanggolin panrgp [-h] [--fasta FASTA] [--anno ANNO]

                         [--clusters CLUSTERS] [-o OUTPUT]

                         [--basename BASENAME] [--rarefaction]

                         [-K NB_OF_PARTITIONS] [--interest INTEREST]

                         [--defrag] [--tmpdir TMPDIR] [--verbose {0,1,2}]

                         [--log LOG] [-c CPU] [-f]

 

Input arguments:

  The possible input arguments :

 

  --fasta FASTA         A tab-separated file listing the organism names, and

                        the fasta filepath of its genomic sequence(s) (the

                        fastas can be compressed). One line per organism. This

                        option can be used alone. (default: None)

  --anno ANNO           A tab-separated file listing the organism names, and

                        the gff filepath of its annotations (the gffs can be

                        compressed). One line per organism. This option can be

                        used alone IF the fasta sequences are in the gff

                        files, otherwise --fasta needs to be used. (default:

                        None)

  --clusters CLUSTERS   a tab-separated file listing the cluster names, the

                        gene IDs, and optionnally whether they are a fragment

                        or not. (default: None)

 

Optional arguments:

  -o OUTPUT, --output OUTPUT

                        Output directory (default:

                        ppanggolin_output_DATE2020-04-10_HOUR14.53.35_PID8500)

  --basename BASENAME   basename for the output file (default: pangenome)

  --rarefaction         Use to compute the rarefaction curves (WARNING: can be

                        time consumming) (default: False)

  -K NB_OF_PARTITIONS, --nb_of_partitions NB_OF_PARTITIONS

                        Number of partitions to use. Must be at least 3. If

                        under 3, it will be detected automatically. (default:

                        -1)

  --interest INTEREST   Comma separated list of elements to flag when drawing

                        and writing hotspots (default: )

  --defrag              Realign gene families to associated fragments with

                        their non-fragmented gene family. (default: False)

 

Common arguments:

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

  --tmpdir TMPDIR       directory for storing temporary files (default:

                        /var/folders/v0/rdv6fcx11fz1dy2jjvthmx1r0000gn/T)

  --verbose {0,1,2}     Indicate verbose level (0 for warning and errors only,

                        1 for info, 2 for debug) (default: 1)

  --log LOG             log output file (default: stdout)

  -c CPU, --cpu CPU     Number of available cpus (default: 1)

  -f, --force           Force writing in output directory and in pangenome

                        output file. (default: False)

> ppanggolin write -h

$ ppanggolin write -h

usage: ppanggolin write [-h] -p PANGENOME -o OUTPUT [--soft_core SOFT_CORE]

                        [--dup_margin DUP_MARGIN] [--gexf] [--light_gexf]

                        [--csv] [--Rtab] [--projection] [--stats]

                        [--partitions] [--compress] [--json] [--regions]

                        [--spots] [--families_tsv] [--all_genes]

                        [--all_prot_families] [--all_gene_families]

                        [--tmpdir TMPDIR] [--verbose {0,1,2}] [--log LOG]

                        [-c CPU] [-f]

 

Required arguments:

  One of the following arguments is required :

 

  -p PANGENOME, --pangenome PANGENOME

                        The pangenome .h5 file (default: None)

  -o OUTPUT, --output OUTPUT

                        Output directory where the file(s) will be written

                        (default: None)

 

Optional arguments:

  --soft_core SOFT_CORE

                        Soft core threshold to use (default: 0.95)

  --dup_margin DUP_MARGIN

                        minimum ratio of organisms in which the family must

                        have multiple genes for it to be considered

                        'duplicated' (default: 0.05)

  --gexf                write a gexf file with all the annotations and all the

                        genes of each gene family (default: False)

  --light_gexf          write a gexf file with the gene families and basic

                        informations about them (default: False)

  --csv                 csv file format as used by Roary, among others. The

                        alternative gene ID will be the partition, if there is

                        one (default: False)

  --Rtab                tabular file for the gene binary presence absence

                        matrix (default: False)

  --projection          a csv file for each organism providing informations on

                        the projection of the graph on the organism (default:

                        False)

  --stats               tsv files with some statistics for each organism and

                        for each gene family (default: False)

  --partitions          list of families belonging to each partition, with one

                        file per partitions and one family per line (default:

                        False)

  --compress            Compress the files in .gz (default: False)

  --json                Writes the graph in a json file format (default:

                        False)

  --regions             Write the RGP in a tab format, one file per genome

                        (default: False)

  --spots               Write spot summary and a list of all rgp in each spot

                        (default: False)

  --families_tsv        Write a tsv file providing the association between

                        genes and gene families (default: False)

  --all_genes           Write all nucleotic CDS sequences (default: False)

  --all_prot_families   Write Write representative proteic sequences of all

                        the gene families (default: False)

  --all_gene_families   Write representative nucleic sequences of all the gene

                        families (default: False)

 

Common arguments:

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

  --tmpdir TMPDIR       directory for storing temporary files (default:

                        /var/folders/v0/rdv6fcx11fz1dy2jjvthmx1r0000gn/T)

  --verbose {0,1,2}     Indicate verbose level (0 for warning and errors only,

                        1 for info, 2 for debug) (default: 1)

  --log LOG             log output file (default: stdout)

  -c CPU, --cpu CPU     Number of available cpus (default: 1)

  -f, --force           Force writing in output directory and in pangenome

                        output file. (default: False)

他のコマンドは省略

 

 

テストラン

ランにはゲノムのテキストファイルを指定する必要がある。テキストファイルの最初のカラムにはユニークな生物名、2 番目の列には関連する FASTA ファイルのパスが含まれていること。

テキストファイルを指定して実行。

git clone https://github.com/labgem/PPanGGOLiN.git
cd PPanGGOLiN/testingDataset/
ppanggolin workflow --fasta organisms.fasta.list --cpu 12

 

興味のある分類群のパンゲノムを様々な方法で記述したファイルや図が得られる。

出力

f:id:kazumaxneo:20200404112413p:plain

 

PPanGGOLiNはMMseqs2を呼び出して、MMseqs2のgreedy set cover algorithmを使用して、すべての遺伝子配列に対してクラスタリングを実行する(--identity(デフォルト0.8)と --coverage(デフォルト0.8)で調整可能)。

matrix.csv

f:id:kazumaxneo:20200404112254p:plain

  

Ushaped_plot.html

f:id:kazumaxneo:20200404112212p:plain

 

tile_plot.html

タイルプロットは、パンゲノムを構成する生物(x軸)の遺伝子ファミリー(y軸)を表すヒートマップ。タイルは、その遺伝子ファミリーが生物に存在する場合は色付きで、存在しない場合は色なしで表示される。

f:id:kazumaxneo:20200404112125p:plain

拡大

f:id:kazumaxneo:20200404112154p:plain

データ量が多く、ブラウザで開くのが困難な場合、「クラウド」遺伝子ファミリーを除外できる。以下のオプションを使用する。全データが収められたpangenome.h5を対象に以下のコマンドを実行する。

ppanggolin draw -p pangenome.h5 --tile_plot --nocloud

 

 

organisms_statistics.tsv

f:id:kazumaxneo:20200404112320p:plain

 

 

gephiで.gexfファイルと_light.gexfファイルを開く。light.gexfファイルは遺伝子ファミリーをノードとし、遺伝子ファミリー間の関係を記述したエッジを含んでいる。、.gexfファイルも同じものを含むが、各遺伝子や遺伝子ファミリー間の関係についてより詳細な情報を含んでいる。

partitionを選択し、適用をクリック。

f:id:kazumaxneo:20200410154525p:plain

 

レイアウトはForceAtlas2を選択。

f:id:kazumaxneo:20200410154527p:plain

 

視覚化した。ノードには、その遺伝子ファミリーに属する遺伝子の数、最も共通する遺伝子名(アノテーションを提供した場合)、最も一般的な産物名(アノテーションを提供した場合)、それが属するパーティションヌクレオチド単位での平均サイズと中央値、この遺伝子ファミリーを持つ生物の数が格納される。エッジにはパンゲノムに存在する回数が記載される。(wikiより)

f:id:kazumaxneo:20200410154522p:plain

 ノードの色は以下の属性を表す。

f:id:kazumaxneo:20200410155015p:plain

 

 

fastaファイルを"--fasta"で指定する代わりに、PROKKAでアノテーションして得たGFFファイル( gff3、.gbk、.gbff、またはそれらの混合)を指定することもできる、

ppanggolin workflow --anno ORGANISMS_ANNOTATION_LIST --cpu 20

#アノテーションにゲノム配列が含まれていない場合両方のオプションを同時に使用する
ppanggolin workflow --anno ORGANISM_ANNOTATION_LIST --fasta ORGANISM_FASTA_LIST --cpu 20

 

小さなパンゲノム(数百ゲノム)は通常のデスクトップコンピュータで簡単に解析できるが、最大のパンゲノムはかなりの量のRAMを必要とする。例えば、20 656ゲノムでは約1日と120GB RAMを必要とした(Githubより)。

 

追記

ステップバイステップで進めることもできる(解説)。例えばアノテーションのみ行うには以下のように打つ。結果は.h5に格納される。

#annotationのみ
ppanggolin annotate --fasta fasta.list --cpu 12

#.h5のstatusを確認
ppanggolin info -p pangenome.h5 --content --status

#予測された全遺伝子を書き出す
ppanggolin write -p pangenome.h5 -o outdir --all_genes

#予測された全遺伝子を書き出す
ppanggolin write -p pangenome.h5 -o outdir --all_genes

 

引用

PPanGGOLiN: Depicting microbial diversity via a partitioned pangenome graph
Guillaume Gautreau, Adelme Bazin, Mathieu Gachet, Rémi Planel, Laura Burlot, Mathieu Dubois, Amandine Perrin, Claudine Médigue, Alexandra Calteau, Stéphane Cruveiller, Catherine Matias, Christophe Ambroise, Eduardo P. C. Rocha, David Vallenet

PLoS Comput Biol. 2020 Mar; 16(3): e1007732

 

 

PPanGGOLiN: Depicting microbial diversity via a Partitioned Pangenome Graph

Guillaume Gautreau, Adelme Bazin, Mathieu Gachet, Rémi Planel, Laura Burlot, Mathieu Dubois, Amandine Perrin, Claudine Médigue, Alexandra Calteau, Stéphane Cruveiller, Catherine Matias, Christophe Ambroise, Eduardo PC Rocha, David Vallenet

bioRxiv, Posted November 09, 2019