kallisto/sleuth使いたい人生だった
[更新版を書いています]
自分のメモ用と，あとは要望があったので
最近インフォ界隈で発現量解析のソフトで盛り上がりを見せているソフトの[kallisto]とそのパイプライン下流の発現変動の有意差検定を行う[sleuth]．早さと正確さが売りになっているみたいです．
patcherラボのブログ記事を見るとそれぞれのアルゴリズムがわかりますkallisto sleuth
僕なりの理解ですと，kallistoではトランスリプトームをk-mer (デフォだと31)にわけてde Bruijn graphを作って，それに対してfastqを当ててる感じです．有名なtophat/bowtieなどは全bpをアライメントしていくので非常に遅くなってしまうのですが，kallistoはこれらのアルゴリズムと比べて爆発的に早くなっています．そのため，なんとbootstrapを使えるようになったのです．
最近うちのラボでもこのkallistoが流行っているんですけど，このパイプライン下流にedgeR/DESeqを使うとあんまりよくないという記事を見たので，せっかくだからpatcherラボが作ったkallisto用の解析ソフトのsleuthも使ってみたいという感じで個々数日いじくりまわしていました．このsleuthもcuffdiffと比べて爆発的に早く，あくまでも個人的な感覚ですけど約10倍近く早いです．アルゴリズム的には，RNA-Seqのfastqをkallistoでトランスクリプトームにマップして得られるexpected_countsから(bootstrapの結果を使って)統計的にtechnical noiseとbiological noise を除いた各transcriptの真のcount値を算出して，(正規分布に乗っ取っているはずなので)そこから有意差を検定しています．
そんな二つのソフトのインストールから使い方をとりあえず載せておきます．
Kallisto
まだ生まれたてのソフトなので (最新がver0.42.3)かなり更新が早いです．
Downloadより最新版のbundle，もしくはsource codeをDLしてください．なおlinux機だと数verが無いので，野良buildしなければ行けません．ですが，僕みたいにrequirement (gcc ver > 4.8; zlib; cmake; HDF5; 僕はgccが4.4でした)が入らない場合は古いバージョンで我慢してください．
Bundle (linux; version 0.42.1)
% wget
https://github.com/pachterlab/kallisto/releases/download/v0.42.1/kallisto_linux-v0.42.1.tar.gz
% tar zxfv https://github.com/pachterlab/kallisto/releases/download/v0.42.1/kallisto_linux-v0.42.1.tar.gz
% cd kallisto_linux-v0.42.1
## ソフトは統一して~/opt/binに入れています
% cp kallisto ~/opt/bin
## Transcriptome index作成
% kallisto index -i <index_name> <transcript_fasta>
# Paired end
% kallisto quant -i <index_name> -o <output_dir> -b <bootstrap_num> R1.fq R2.fq
# Single read
% kallisto quant -i <index_name> -o <output_dir> -b <bootstrap_num> --single -l <fragment_length>
single endですと--singleと-l を入れないと行けません．fragment lengthはシークエンスライブラリーを作成した人に聞いてください．以下使う可能性が高いオプションです．なおこの先のsleuthは-bを使います (制作者は大体100あれば大丈夫という風に言っています)．
-t thread数の指定 (v0.42.2以降)
--pseudobam alignmentのSAMでの出力
--bias bias補正を行う
-b bootstrap数
これを前サンプルに対して行い(僕はshell scriptで全部bgで一気に流しています)，一つのdirectoryにまとめておきます．このまとめておくのがsleuthに必要です．ここでは~/k_resultsとします．
sleuth
kallistoより遅れて出されているんですけど同様に更新が早いのでupdateに気をつけてください．現状v0.27.3です．
とりあえずインストールからです．
source("http://bioconductor.org/biocLite.R")
biocLite("rhdf5")
install.packages("devtools")
devtools::install_github("pachterlab/sleuth")
僕はこのdevtoolsをインストールするときにgit2rが入らなくて結構苦労しました，結局libiconvを再びインストールしたら通りました．
ここから解析に入りますが，sample-to-conditionのファイルが必要です (一例です)．
sample condition
A1 Control
A2 Control
A3 Control
B1 Sample
B2 Sample
B3 Sample
そして解析に入ります．
## ライブラリの読み込み
library(sleuth)
## 結果ファイルのdirectory
base_dir <- "~/k_results"
# サンプル名の取得
sample_id <- dir(file.path(base_dir))
# サンプル名とそれぞれのkallisto結果directoryの連結
kal_dirs <- sapply(sample_id, function(id) file.path(base_dir, id))
# sample-to-condition table
s2c<-read.delim("sample_table.txt", header=T)
# kallistoの出力をcondition依存でsleuth object (so)として読み込む
so <- sleuth_prep(kal_dirs, s2c, ~ condition)
# <標準出力>
reading in kallisto results
............
normalizing est_counts
xxxxxx targets passed the filter
normalizing tpm
merging in metadata
normalizing bootstrap samples
# </標準出力>
# soをモデルにfitさせる
so <- sleuth_fit(so)
#<標準出力>
summarizing bootstraps
fitting measurement error models
shrinkage estimation
computing variance of betas
# </標準出力>
# fitされたモデルの確認．
models(so)
# <標準出力>
> models(so)
[ full ]
formula: ~condition
coefficients:
(Intercept)
conditionSample
no tests found.
# </標準出力>
# test
so <- sleuth_test(so, which_beta = 'conditionSample')
models(so)
# <標準出力>
> models(so)
[ full ]
formula: ~condition
coefficients:
(Intercept)
conditionSample
tests:
conditionSample
# </標準出力>
# 結果の出力
results_table <- sleuth_results(so, 'conditionSample')
ここで色々とわかりにくいのがmodels(so)とsleuth_test，sleuth_resultsです．
models(so)はsleuth_fitでfitされたモデルの一覧をくれます．そのときの結果 (上記のsample_to_conditionを一例に)はこのようになります．
> models(so)
[ full ]
formula: ~condition
coefficients:
(Intercept)
conditionSample
no tests found
これはControlとSampleの二群間比較なのでこうなりますが，複数群間でやろうとするとこうなります．
> models(so)
[ full ]
formula: ~condition
coefficients:
(Intercept)
conditionSample1
conditionSample2
conditionSample3
・
・
・
no tests found
# </標準出力>ですが，一つ注意しないといけないのは，これらはすべて"Control"との比較でのモデルです．現在のversionではこうですが，今後多群間の前組み合わせの比較が実装されるかも知れません．どのconditionXXXがどのサンプル間を比較しているのかを確認するにはdesign_matirx()を使えばわかります．
またsleuth_testでは，which_betaを指定する事でどの(design_matrixに基づいた)組み合わせで検定を行うかを指定します．*****一回のsleuth_testでは一つのwhich_betaしか指定できませんが，which_betaを変えながら複数回sleuth_testを行う事は可能です*****．そして気づくのに時間かかりましたが，このsleuth_testしなければsleuth_resultsで書き出せません．
ここまでいけば，q-value (BH補正されたp-value)がでるのでそこから機能解析やエンリッチメント解析などに進めます．また，sleuthにはデフォルトで様々なplot関数があるので，これを使ってもよし，shinyというさむしんぐ (sleuth_live)を使ってGUIで解析する事ができます．
この記事で一番言いたかったのは，
Kallistoは早い．
Sleuthは早い．
どっちもラップトップでも早い．
やべぇ．
です．
