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

RNA-Seqの定量化の二大定番ツールとして、SalmonとKallistoがある。
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")