macでインフォマティクス

macでインフォマティクス

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

動的に生成されるRスクリプトを用いてバルクRNA-seqの自動探索と可視化を行う Searchlight2

2022 1/5 複数比較の例追記、コマンドの誤字修正

 

 バルクRNA-seqデータが処理されると、すなわちアラインメントされ、発現および差分表が作成されると、生物学的性質の探索、視覚化および解釈が行われる重要なプロセスが残る。可視化・解釈パイプラインを使用しないと、このステップは時間と手間がかかり、多くの場合Rを使用して完了する。市販の可視化・解釈パイプラインは包括的だが、自由に使用できるパイプラインは現在限られている。ここでは、自由に利用できるバルクRNA-seq可視化および解釈パイプラインであるSearchlightを紹介する。Searchlightは、グローバル、パスウェイ、単一遺伝子レベルに焦点を当てた包括的な統計および視覚的分析、3つのワークフローによる生物または実験の複雑さに関係なくほとんどの差分実験デザインとの互換性、レポート、およびユーザーフレンドリーなRスクリプトとShinyアプリによるプロットの下流ユーザー修正のサポートなどを提供する。Searchlightは、現在のベストツール(VIPERとBioJupies)よりも優れた自動化を提供することを示す。標準的なバルクRNA-seq処理パイプラインと並行して、Searchlightを使用すると、論文原稿品質の図まで含めてバルクRNA-seqプロジェクトやtimed re-analysis studyを3時間未満で完了できることを実証した。

 

Githubより

 バルクRNA-seqデータの処理、すなわちアライメント、発現および発現変動遺伝子のテーブルの作成が完了すると、生物学の探究、視覚化、解釈を行う重要なプロセスが残ります。最終的には、報告書、論文、原稿の図表を作成します。この下流工程を完了するための典型的な方法は、手動でコード化したコマンドラインとRベース(または類似の)の分析です。しかし、この方法では、数日から数週間を要することもあり、大変な労力を要します。

 Searchlight2は、バルクRNA-seqの探索、可視化、解釈のパイプラインで、下流域の解析を自動化することを目的としています。標準的なアライメントおよび処理パイプライン(Star2、Hisat2、Kallisto、DEseq2、EdgeRなど)と共に使用すれば、研究者はバルクRNA-seqプロジェクトを数時間のうちに、。Rによる手動解析と区別がつかないほどの水準で、最小限の労力で完了させることが可能です。

複雑なパイプラインを使用したり、理解したりする必要はありません。その強みは

  • 強力で広く利用されている幅広い分析・可視化手法。
  • 発現、差分発現、シグネチャー解析をカバーする独立したワークフローを使用できる。これらは、デザイン、複雑さ、生物学、生物に関係なく、実験との互換性を提供する。また、解析を簡素化する。
  • 包括的なレポートを作成。
  • すべてのプロットにRとR Shinyを使用しているため、バイオインフォマティシャンやウェットラボの科学者がすべてのプロットを視覚的に簡単に変更することができる。
  • R-snippetデータベースにより、ユーザーはすべてのプロットのデフォルトの外観を自分の好みに合わせて変更することができる。
  • Searchlight2は100%自動化されている。
  • Searchlight2は、サンプルシート、発現マトリックス、任意の数の発現差テーブルなど、典型的なRNA-seqダウンストリーム解析の入力を受け付ける。
  • バイオインフォマティシャン、RNA-seqサービスプロバイダー、ベンチサイエンティストが最小限の労力でバルクRNA-seq研究プロジェクトを迅速に進め、さらなる詳細解析や別の解析アプローチにリソースを解放できるように設計されている。
  • Searchlight2はアライメントやリードのカウント、発現量や差分発現量の計算を行わないため、処理パイプラインではないことに留意する。これらのステージは、Searchlight2を使用する前に完了している必要がある。
  • 発現量の行列(TPM、RPKM、Rlogなど)と少なくとも1つの差分発現テーブル(DEseq2、EdgeRなど)があれば、どんな処理パイプラインでもSearchlight2を使用できる。

 

インストール

ubuntu18.04でcondaを使ってR4,0の環境を作ってテストした。

Github

#依存
mamba create -n searchlight2 python=2.7 -y
conda activate searchlight2
mamba install Scipy Numpy -y
mamba install -c bioconda r -y

#Rの依存
install.packages("ggplot2")
install.packages("reshape")
install.packages("amap")
install.packages("amap")
install.packages("gridExtra")
install.packages("gtable")
install.packages("ggally")
install.packages("network")
install.packages("sna")

#shinyを使う場合(一番下で紹介)の依存;Rstudio上でインストールする
install.packages("shiny")
install.packages("shinyFiles")
install.packages("fs")
install.packages("shinycssloaders")
install.packages("graphics")
install.packages("dplyr")
install.packages("GGally")

#Searchlight2本体
cd Searchlight2-master/software/

python Searchlight2.py -h

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

#####       Searchlight 2       #####

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

 

Automated bulk RNA-seq data exploration and visualisation using dynamically generated R-scripts

John J. Cole & Carl S. Goodyear, et al.

University of Glasgow (Scotland) 2021

 

v2.0.0

option -h not recognized

 

 

 

実行方法

Searchlight2は1つのコマンドで実行される。まず、入力ファイルを検証し、それらを1つの「マスター遺伝子テーブル」に結合し、そこから下流の解析が行われる。次に、各ワークフローを繰り返し、中間ファイル、統計解析結果ファイル、プロットごと、ワークフローごとのRスクリプト、プロット、HMTLによるレポート、そして最後にShinyアプリが生成される。

Searchlight2のランには最小構成で4つのファイルが必要。入力のフォーマットには厳密で、ファイル間の対応関係が違うとエラーを起こす。例えば、ファイル間でサンプル名が微妙に間違っているケースは許容されない。

 

必要なファイル1- 発現量のマトリックスTPM、RPKM、Rlogなど)

最初の列は遺伝子ID(Ensembl、Refseqなど)とする。1行目は、最初のセルが"ID "、残りはサンプル名とするヘッダー行。サンプル名は数字で始めることはできず、数字、文字、アンダースコア(_)のみ使用可能。

レポジトリの例

f:id:kazumaxneo:20220104021448p:plain

 

必要なファイル2- Differential expression table

DESeq2やEdgeRなどを使った標準的な発現変動遺伝子のテーブル。遺伝子は行単位で、列は遺伝子ID、log2 fold change、p-value、adjusted p-value(この順)のみに切り詰めたもの。1行目にヘッダ行が必要で、ヘッダは正確に”ID”、”Log2Fold”、”p”、”p.adj”でなければならない。大文字・小文字は区別しない、引用符は無視する。IDの型は、発現マトリックスと同じでなければならない。

レポジトリの例

f:id:kazumaxneo:20220104022039p:plain

 

必要なファイル3- サンプルシート (SS)

標準的なタブ区切りのサンプルシートで、1列目に各サンプル名、2列目に所属するサンプルグループを記載する。ヘッダ行があり、ヘッダは正確に"sample", "sample_group"でなければならない。大文字と小文字の区別はなく、引用符も無視される。サンプルグループの階層が複数ある場合(細胞の種類、治療法、年齢など)、3列名以降に、任意のヘッダー(先頭行)から始めてさらにグループを記載できる。サンプル名とサンプルグループ名は数字で始めることはできず、数字、文字、アンダースコア(_)のみを使用することができる。

f:id:kazumaxneo:20220104022553p:plain

 

必要なファイル4- バックグラウンドファイル(BG)

全遺伝子がリストアップされている典型的なバックグラウンドアノテーションファイル。EnsemblsのBiomartから簡単に作成することができる。遺伝子は行単位で、特定のアノテーションは列単位で記述する。ファイルには以下のアノテーションが必要。Gene ID、Gene Symbol、Chromosome、Start position、Stop position、Biotype(遺伝子タイプ)。ヘッダ行には、”ID”、”Symbol”、”Chromosome”、”Start”、”Stop”、”Biotype”ヘッダを正確に記述する。大文字と小文字の区別はなく、引用符は無視する。遺伝子記号が不明な場合は、IDを使用する。Biotypeがわからない場合は、その列のすべてのセルに "gene "と記入する。

f:id:kazumaxneo:20220104095157p:plain

 

任意で必要なファイル5- 遺伝子セットデータベースファイル

GMT形式の有効な遺伝子セットデータベースファイル(GO、string、KEGGなど)を指する(例;Searchlight2/gene_set_databases/)。1つの遺伝子セットにつき1行で、最初のセルが遺伝子セット名(細胞周期など)、2番目のセルが任意のメモ(ない場合は単にNAと入力)、そして遺伝子セット内の各遺伝子で構成されている。

GOデータベースの例

https://github.com/Searchlight2/Searchlight2/blob/master/gene_set_databases/GO_bp_human.tsv

f:id:kazumaxneo:20220104100157p:plain

 

テストラン

必要な4つのオプションを指定する。過剰発現解析も行うなら--oraオプションで必要なファイル5を指定する。--deと--oraは複数のパラメータを指定する。複数パラメータ書くオプションは、パラメータ間にスペースを入れてはならない

python Searchlight2-master/software/Searchlight2.py --out path=results \
--bg file=Searchlight2/backgrounds/mouse_GRCm38.p6.tsv \
--em file=Searchlight2/sample_datasets/EM.tsv \
--ss file=Searchlight2/sample_datasets/SS.tsv \
--de file=Searchlight2/sample_datasets/WT_vs_KO.tsv,numerator=KO,denominator=WT
--ora file=Searchlight2/gene_set_databases/GO_bp_mouse.tsv,type=GO_bp

(注;すべてのパスは「root」からの完全なものでなければならない)

過剰発現解析が必要ないなら--oraオプションを消す。上流調節因子解析も含める場合は、コマンドに --ura パラメータを追加する。

 

出力例

results/

f:id:kazumaxneo:20220103224136p:plain

results/all_genes/

f:id:kazumaxneo:20220103224216p:plain

results/all_genes/de_workflows/KO_vs_WT/

f:id:kazumaxneo:20220103224252p:plain

 

report.htmlを開く。

 

f:id:kazumaxneo:20220103224104p:plain

f:id:kazumaxneo:20220103234019p:plain

f:id:kazumaxneo:20220103234030p:plain

f:id:kazumaxneo:20220103234045p:plain

f:id:kazumaxneo:20220103234059p:plain

f:id:kazumaxneo:20220103234109p:plain

f:id:kazumaxneo:20220103234116p:plain

f:id:kazumaxneo:20220103234136p:plain

f:id:kazumaxneo:20220103234145p:plain

f:id:kazumaxneo:20220103234158p:plain

f:id:kazumaxneo:20220103234214p:plain

f:id:kazumaxneo:20220103234228p:plain

f:id:kazumaxneo:20220103234240p:plain

f:id:kazumaxneo:20220103234306p:plain

f:id:kazumaxneo:20220103234340p:plain

f:id:kazumaxneo:20220103234403p:plain

f:id:kazumaxneo:20220103234422p:plain

 

 

Rに不慣れなユーザーのためにR shinyの環境も用意されている。利用するには出力のshiny/下にあるserver.rをRstudioで開く。

f:id:kazumaxneo:20220103222344p:plain

server.rが読み込まれたら、Rstudio左上のパネルの右上にあるスタートボタンをクリックする。

f:id:kazumaxneo:20220103222516p:plain

 

ブラウザが立ち上がる。

f:id:kazumaxneo:20220103222701p:plain

 

遺伝子、ワークフロータイプを選択する。

f:id:kazumaxneo:20220103222758p:plain

分析方法を選択するとそのまま視覚化される。

f:id:kazumaxneo:20220103222835p:plain

PCA;各主成分の寄与率の推移

f:id:kazumaxneo:20220103222916p:plain

 

作図パラメータは細かい部分まで調整可能。図のサイズを大きくした。

f:id:kazumaxneo:20220103223405p:plain

 

レポートhtml付きサンプル解析結果がsample_datasets/results.zipに含まれています。興味ある方はアクセスしてみて下さい。

Rの依存関係が正しく導入されていない場合、各図の出力には失敗しますが、Searchlight2の処理は最後まで実行されます。その場合、全体のログを見てもどこが失敗したのかわかりにくいわけですが、Searchlightは出力の各ディレクトリにRの作図コードが含まれているので、Rで実行してどんなエラーが出るか確認してみると原因ははっきりします(例えばcd xxx/xxx/PCA/ してRを立ち上げ、source("pca.R"))。

 

 

2022 1/5 追記

2種類以上のサンプルを持つために複数の比較が必要な場合、Searchlight2はこれらのデータセットを完全にサポートしている。これらは、単に「アップ」または「ダウン」制御を超えたシグネチャを含むことができるため、「シグネチャベース」と呼ばれている。

 

独立した2群間比較を複数ペアで行う。

複数の差分発現ファイル(例:WT vs KO、KO vs KO_rescue)がある場合、その分だけ--deパラメータを追加するだけです。これだけで他の比較も同時に実行できます。

サンプルデータセットだとWT vs KO、KO vs KO_rescueで比較したいので、2つ--deを指定しています。以下のコマンドになります。

python Searchlight2/software/Searchlight2.py --out path=results \
--bg file=Searchlight2/backgrounds/mouse_GRCm38.p6.tsv \
--em file=Searchlight2/sample_datasets/EM.tsv \
--ss file=Searchlight2/sample_datasets/SS.tsv \
--de file=Searchlight2/sample_datasets/WT_vs_KO.tsv,numerator=KO,denominator=WT \
--de file=Searchlight2/sample_datasets/KO_rescue_vs_KO.tsv,numerator=KO_rescue,denominator=KO

出力例

f:id:kazumaxneo:20220105103606p:plain

 

Githubでは、このような独立したいくつかの2群間比較に加えて、3つのグループを一緒に調べる手順についても説明されています。具体的にはmultiple differential expression (MDE)ワークフローを使用することで実現できるようになっており、-mdeパラメータを使用します。-mdeパラメータは、2つ以上の差分比較を同時に行うことを可能にするものです。コマンドが複雑になりますので、レポジトリのIncluding a formal signature analysisのセクションを読んでみて下さい。

引用

Searchlight: automated bulk RNA-seq exploration and visualisation using dynamically generated R scripts
John J Cole, Bekir A Faydaci, David McGuinness, Robin Shaw, Rose A Maciewicz, Neil A Robertson, Carl  Goodyear

BMC Bioinformatics. 2021 Aug 19;22(1):411

 

参考