Kallisto下游推薦的分析軟件是Sleuth哨免，所以我們本次介紹如何使用Sleuth進行相關(guān)的分析。
RNA-seq在完成轉(zhuǎn)錄本的鑒定及定量后呛伴，我們經(jīng)常做的分析之一就是差異表達霜幼。通常我們使用較多的有以下工具：
而Kallisto主要推薦的下游分析軟件是Sleuth這個R包。Sleuth與其他常用的包進行比較其結(jié)果還是很不錯的蜻牢。
下面我們來介紹一下該包的使用。
Sleuth的安裝
source("http://bioconductor.org/biocLite.R")
biocLite("rhdf5")
install.packages("devtools")
devtools::install_github("pachterlab/sleuth")
如果你安裝了conda偏陪，也可以通過命令行的方式安裝conda install --channel bioconda r-sleuth
Sleuth的使用
#設(shè)置kallisto生成的路徑
base_dir <- "/home/617/sleuth_analysis/kallisto_qaunt_output"
#獲取所有的simple_id,具體可以看一些file.path命令的使用
sample_id <- dir(file.path(base_dir))
sample_id
## [1] "WT_1" "WT_2" "WT_3" Mutation_1"
## [5] "Mutation_2" Muatation_3"
#獲取結(jié)果文件的路徑
kal_dirs <- sapply(sample_id, function(id) file.path(base_dir, id))
kal_dirs
## WT_1
## "/home/617/sleuth_analysis/kallisto_qaunt_output/WT_1"
## WT_2
## "/home/617/sleuth_analysis/kallisto_qaunt_output/WT_2"
## WT_3
## "/home/617/sleuth_analysis/kallisto_qaunt_output/WT_3"
## Mutation_1
## "/home/617/sleuth_analysis/kallisto_qaunt_output/Mutation_1"
## Mutation_2
## "/home/617/sleuth_analysis/kallisto_qaunt_output/Mutation_2"
## Mutation_3
## "/home/617/sleuth_analysis/kallisto_qaunt_output/Mutation_3"
#讀取實驗設(shè)計表
s2c <- read.table("./design_matrix.txt", header = TRUE, sep='\t',stringsAsFactors=FALSE)
s2c
## sample condition reps
## 1 WT_1 WT 1
## 2 WT_2 WT 2
## 3 WT_3 WT 3
## 4 Mutation_1 Mutation 1
## 5 Mutation_2 Mutation 2
## 6 Mutation_3 Mutation 3
#與路徑合并
s2c <- dplyr::mutate(s2c, path = kal_dirs)
print(s2c)
## sample condition reps
## 1 WT_1 WT 1
## 2 WT_2 WT 2
## 3 WT_3 WT 3
## 4 Mutation_1 Mutation 1
## 5 Mutation_2 Mutation 2
## 6 Mutation_3 Mutation 3
## path
## 1 /home/617/sleuth_analysis/kallisto_qaunt_output/WT_1
## 2 /home/617/sleuth_analysis/kallisto_qaunt_output/WT_2
## 3 /home/617/sleuth_analysis/kallisto_qaunt_output/WT_3
## 4 /home/617/sleuth_analysis/kallisto_qaunt_output/Mutation_1
## 5 /home/617/sleuth_analysis/kallisto_qaunt_output/IMutation_2
## 6 /home/617/sleuth_analysis/kallisto_qaunt_output/Mutation_3
#獲取轉(zhuǎn)錄本ID信息
#方法一：基于biomart
library(biomaRt)
#change to genename
mart <- biomaRt::useMart(biomart = "ENSEMBL_MART_ENSEMBL",
dataset = "hsapiens_gene_ensembl",
host = 'ensembl.org')
t2g <- biomaRt::getBM(attributes = c("ensembl_transcript_id", "ensembl_gene_id",
"external_gene_name"), mart = mart)
t2g <- dplyr::rename(t2g, target_id = ensembl_transcript_id,
ens_gene = ensembl_gene_id, ext_gene = external_gene_name)
#方法二：下載EnsDb.Hsapiens.v75或者其他版本相應(yīng)的R包
library('EnsDb.Hsapiens.v75')
library('ensembldb')
mmsdb<- EnsDb.Hsapiens.v75
tx.mms<- transcripts(mmsdb, return.type = "DataFrame")
tx.mms<- tx.mms[, c("tx_id", "gene_id")]
gene.mms<- genes(mmsdb, return.type = "DataFrame")
gene.mms<- gene.mms[, c("gene_id", "symbol")]
tx2genes.mms <- merge(tx.mms, gene.mms, by = "gene_id")
t2g <- tx2genes.mms[, 2:3]
#方法三：自己制作相應(yīng)的txt文件第一列為target_id孩饼，第二列為對應(yīng)的gene（可以用ID或者是Symbol）
t2g <- read.table("./t2g.txt", header = TRUE, stringsAsFactors=FALSE)
#讀取Kallisto結(jié)果文件
#如果之前在kallisto quant沒有設(shè)置-b參數(shù)該步會報錯
library(sleuth)
so <- sleuth_prep(s2c, ~ condition, target_mapping = t2g, extra_bootstrap_summary = TRUE)
#如果想要將一個基因的不同轉(zhuǎn)錄本合并，可以按照以下模式導入
so <- sleuth_prep(s2c, ~ condition, target_mapping = t2g, extra_bootstrap_summary = TRUE,gene_mode=T,aggregation_column = 'ext_gene')
#建立模型竹挡，該步通過依據(jù)不同實驗條件下（WT與Mutation）的數(shù)據(jù)構(gòu)建線性模型，以“smooth”每個樣品的原始kallisto豐度估計值立膛。
so <- sleuth_fit(so)
#為了鑒定不同組差異表達的轉(zhuǎn)錄本揪罕，sleuth進行第二次擬合即“reduced”模型，該模型假設(shè)兩種條件下的豐度相等宝泵。然后Sleuth將識別與“full”模型明顯更好擬合的轉(zhuǎn)錄本好啰，即在假設(shè)豐度相等的情況下與reduced模型擬合度低，而加入變量不同的實驗條件后與full模型的擬合度高儿奶，說明該轉(zhuǎn)錄本/基因表達有差異框往。
so <- sleuth_fit(so, ~1, 'reduced')
#Sleuth采取 likelihood ratio test(LRT) 進行鑒定
so <- sleuth_lrt(so, 'reduced', 'full')
#可以使用models()查看不同的模型
models(so)
## [ full ]
## formula: ~condition
## coefficients:
## (Intercept)
## conditionWW
## [ reduced ]
## formula: ~1
## coefficients:
## (Intercept)
# LRT檢驗結(jié)果
sleuth_table <- sleuth_results(so, 'reduced:full', 'lrt', show_all = FALSE)
sleuth_significant <- dplyr::filter(sleuth_table, qval <= 0.05)
#Sleuth還提供了Wald test。該函數(shù)計算每個轉(zhuǎn)錄本上一個特定'beta'系數(shù)的Wald檢驗;這提供了β值闯捎，它們是成對比較的fold change偏差值椰弊。
so <- sleuth_wt(so,'conditionWT','full')
results_table <- sleuth_results(so, 'conditionWT', test_type = 'wt')
sleuth_significant <- dplyr::filter(results_table, qval <= 0.05)
#如何導出基因表達量TPM
sleuth_geneexp<-sleuth_to_matrix(so,'obs_norm', 'tpm')
以上就是Sleuth 的應(yīng)用。如果你依然還是想要使用經(jīng)典的DESeq2瓤鼻，你可以這樣導入文件：
library(tximport)
library(DESeq2)
txi<-tximport(files=kal_dirs,type='kallisto',tx2gene=t2g)
dds1<-DESeqDataSetFromTximport(txi,s2c,~condition)
然后就可使用DESeq2進行下游的分析秉版。
