Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
Help us understand the problem. What is going on with this article?

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

More than 1 year has passed since last update.

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

antiplastics
書くぜーどんどん書くぜー(書いてない
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away