Edited at

Sailfish/Salmonのブートストラップ法の結果をwasabiでsleuthに渡す

RNA-Seqの定量化の二大定番ツールとして、SalmonKallistoがある。

SalmonはSailfishという以前のツールの拡張である。

これらは、転写産物ごとのカウント値やTPM値などの発現量の単位を、リファレンスへのアライメントを経由せずに高速に計算する( https://www.slideshare.net/antiplastics/rnaskim )。

現状、計算速度、メモリ使用量、ユーザビリティ(使いやすさ、ドキュメントの質、バグへの対応...etc)など、全てにおいて互角の戦いをしている。

また、定量後の値をどのように、その下流の解析に使うのかというのも重要だが、周辺ツールも同程度に充実している。

例えば、tximportというR/Bioconductorパッケージを使えば、ランごとに生成されるSalmon、Kallistoの結果をまとめて、一つのデータフレームに変換してくれる。

転写産物で定量化した結果を、遺伝子レベルに要約して下流に使いたいという要望もあるが、tximportはtx2geneというオプションで遺伝子レベルに要約することもできる。

そもそも、どのように要約するべきか、そもそも要約していいのか、といった議論があるが、この辺りのことは以下の論文で詳細を調べており、この論文で検証したscaledTPMとlengthScaledTPMという二つの要約法が利用可能である。

Michael Love, et. al., Swimming downstream: statistical analysis of differential transcript usage following Salmon quantification, F1000Research, 2018

一方、sleuthというツールは、Kallistoのチームが開発したRパッケージであるが、Kallistoの結果しかサポートしていない。

このツールは、Kallistoの結果をもとに、発現変動遺伝子(DEG)の検出や、ブートストラップ法での定量値の信頼区間の計算を行う。

しかし面白いことに、ここでwasabiというSalmonのチームが開発したツールを使うと、Sailfish/Salmonの結果を、sleuthが受け付けているabundance.h5というHDF5ファイルに変換してくれるため、無理やり(?)Sailfish/Salmonの結果をsleuthに渡すことができてしまう。

DEGの解析は、tximportを使えば、edgeR、DESeq2、Limmaなどの定番の手法に渡せるが、ブートストラップ法の結果の可視化は、現状sleuthしかサポートしていないと思うので、このツールはなかなか便利だと思い、試しに使ってみた。

なお、ここでは、すでに、

analysis/salmon_options/SRR{3670977..3670992}

という16箇所のディレクトリにSailfishの結果(quant.sf、aux_infoなど)が出力されているものとする。

また、現状、ブートストラップ法をしていないSailfish/Salmon/Kallistoの結果は、sleuthは受け付けていないため注意。

https://www.biostars.org/p/274933/

wasabiの使い方は簡単で、prepare_fish_for_sleuth(fishは多分、Sailfish/Salmonをまとめて言っている)という関数に、Sailfish/Salmonの結果が置かれたディレクトリを指定するだけである。

# パッケージダウンロード

install.packages("devtools", repos="http://cran.r-project.org")
devtools::install_github("COMBINE-lab/wasabi")

# パッケージロード
library("wasabi")

# Salmon/Kallisto解析結果ファイル
srr <- paste0("SRR", 3670977:3670992)
file_salmon_options <- file.path("analysis/salmon_options", srr)

# うまくいくとabundance.h5が置かれる => これがsleuthで使える
prepare_fish_for_sleuth(file_salmon_options)

これにより、abundance.h5が生成される。

h5ls analysis/salmon_options/SRR3670977/abundance.h5

# aux Group
# bootstrap Group
# est_counts Dataset {207749}

続いて、このファイルをsleuthに渡す。

ここでは、sample(サンプル名)、condition(実験条件の違い)、path(ファイルの置き場所)というデータフレームを用意して、sleuth_prepという関数に渡す。

# パッケージダウンロード

install.packages("devtools", repos="http://cran.r-project.org")
devtools::install_github("pachterlab/sleuth")

# パッケージロード
library("sleuth")

# Salmon/Kallisto解析結果ファイル
srr <- paste0("SRR", 3670977:3670992)
condition <- c(rep("MAQCA_1", 4), rep("MAQCA_2", 4),
rep("MAQCB_1", 4), rep("MAQCB_2", 4))

# SleuthでBootstrap解析
s2c_salmon_options <- data.frame(sample=srr, condition=condition)
so_salmon_options <- sleuth_prep(s2c_salmon_options,
extra_bootstrap_summary = TRUE, read_bootstrap_tpm = TRUE)

これにより、Sailfish/Salmonの結果がsleuthのRオブジェクトに変換されたので、あとは普通にsleuthの各種関数を適用できる。

例えば、get_bootstrap_summaryという関数では、ブートストラップ法の結果を閲覧でき、plot_bootstrapでは、ブートストラップのばらつきを箱ひげ図で見ることができる。

# ハウスキーピング遺伝子を確認

get_bootstrap_summary(so_salmon_options,
"ENST00000229239.10|ENSG00000111640.15|OTTHUMG00000137379.3|OTTHUMT00000268059.2|GAPDH-201|GAPDH|1285|protein_coding|",
unit="tpm")

plot_bootstrap(so_salmon_options,
"ENST00000229239.10|ENSG00000111640.15|OTTHUMG00000137379.3|OTTHUMT00000268059.2|GAPDH-201|GAPDH|1285|protein_coding|",
color_by="condition", unit="tpm")

スクリーンショット 2019-06-12 15.11.52.png