kallisto-sleuth workflow in dealing with sequencing data

Needless to say, kallisto-sleuth workflow opens up a new era in dealing with sequencing data; it dramatically scales down the computational complexity to the point which sequencing can be handled by a PC. I find this article is a good startpoint to understand the mechanism of kallisto

As the downstream of kallisto, sleuth has been optimised to deal with kallisto output. An attactive concept from sleuth is “variance decoupling”, to “decouple biological variance from inferential variance.” Inferential variance, according to online method, consists of “errors in mesurement and subsequent inference”, which can be obtained through bootstrapping from kallisto; inferential variance is eliminated, or decoupled away before shrinkage (and later added back when estimating the fixed-effect parameter). This articles walks through R code that realises this decoupling process.

Dataset used and method to “open up packages”

cd /your_dirs_to_Ecoli_kallisto_files

kallisto index -i transcripts.idx Ecoli_transcripts.fasta

kallisto quant -b 100 -i transcripts.idx -o LB1 --single -l 500 -s 50 LB1.fastq.gz
kallisto quant -b 100 -i transcripts.idx -o LB2 --single -l 500 -s 50 LB2.fastq.gz
kallisto quant -b 100 -i transcripts.idx -o LB3 --single -l 500 -s 50 LB3.fastq.gz
kallisto quant -b 100 -i transcripts.idx -o MG1 --single -l 500 -s 50 MG1.fastq.gz
kallisto quant -b 100 -i transcripts.idx -o MG2 --single -l 500 -s 50 MG2.fastq.gz
kallisto quant -b 100 -i transcripts.idx -o MG3 --single -l 500 -s 50 MG3.fastq.gz

## setwd(paste0(getwd(), "/content/package")) ## needed to run this code chunk alone instead of Rmd compilation. 
if (packageVersion("sleuth") != "0.29.0"){
    devtools::install_github("pachterlab/sleuth@v0.29.0")
}

library(tidyverse)
library(sleuth)
library(kableExtra)

source("../../static/sleuth/sleuth_internal_structure.R")

## read the experiment design matrix
t <- read.table("../../static/sleuth/tbl.txt", header = T)
t$path <- as.character(t$path)

obj <<- test <- PPP_sleuth_prep(t, extra_bootstrap_summary = TRUE, num_cores = 1)

## Warning: `select_()` is deprecated as of dplyr 0.7.0.
## Please use `select()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

test <- FFF_sleuth_fit(test, ~condition, 'full')
test <- RRR_sleuth_fit(test, ~1, 'reduced')
test <- LLL_sleuth_lrt(test, 'reduced', 'full')

The walkthrough

options(knitr.table.format = "html", kableExtra.auto_format = TRUE, kableExtra.font_size = 8)
fs <- 8

sleuth_prep output is a list of 19; the observed count and bootstrapping information is stored in bs_summary.

str(obj$bs_summary)

## List of 2
##  $ obs_counts: num [1:1481, 1:6] 4.62 3.13 3.53 2.46 5.6 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:1481] "b0002" "b0003" "b0004" "b0006" ...
##   .. ..$ : chr [1:6] "SRR493366" "SRR493367" "SRR493368" "SRR493369" ...
##  $ sigma_q_sq: Named num [1:1481] 0.02105 0.07362 0.03354 0.049 0.00279 ...
##   ..- attr(*, "names")= chr [1:1481] "b0002" "b0003" "b0004" "b0006" ...

Within sleuth_fit, the function sleuth:::me_model wrapped in sleuth:::me_model_by_row elimiates inferential variance from total variance. Inferrential variance (variable sigma_q_sq) is obtained from rowMean of the bs_summary.

## modified from me_model and me_model_by_row
me_model_by_row <- lapply(1:nrow(bs_summary$obs_counts), function(i) {
  function(X, bs_summary$obs_counts[i, ], bs_summary$sigma_q_sq[i]) {
    n <- nrow(X)
    degrees_free <- n - ncol(X)
    ols_fit <- lm.fit(X, y)
    rss <- sum(ols_fit$residuals^2) # residual sum sq
    
    ## (raw) biological_var = total_var - inferential_var
    sigma_sq <- rss/(degrees_free) - sigma_q_sq 
    
    ## collectin result 
    mean_obs <- mean(y)
    var_obs <- var(y)
    list(ols_fit = ols_fit, b1 = ols_fit$coefficients[2], rss = rss, 
         sigma_sq = sigma_sq, sigma_q_sq = sigma_q_sq, mean_obs = mean_obs, 
         var_obs = var_obs)
    }
})

Subsequently, sigma_sq_pmax is obtained from sigma_sq of each gene by sigma_sq_pmax = pmax(sigma_sq,0). This is to prevent negative variance of count.

head(FFF_mes_df) %>%
  kable(.) %>%
  kable_styling(font_size = fs) 

sleuth:::sliding_window_grouping (uses ecdf) and sleuth:::shrink_df perform loess prediction using biological variance.

head(FFF_swg) %>%
  kable(., "html") %>%
  kable_styling(font_size = fs)

head(FFF_shrink1) %>%
  kable(., "html") %>%
  kable_styling(font_size = fs)

A power of 4 is performed on shrink variable to produce smooth_sigma_sq.

head(FFF_shrink_select) %>%
  kable(., "html") %>%
  kable_styling(font_size = fs)

The maximun of variance was obtained between raw and shrunken estimate.

head(FFF_shrink_mutate) %>%
  kable(., "html") %>%
  kable_styling(font_size = fs)

The fixed-effect beta is estimated with sleuth:::covar_beta using total vairance, which is the shrunken biological variance smooth_sigma_sq_pmax plus inference variance sigma_q_sq.

In the later part, the likelihood of each model is calculated by

all_likelihood <- sapply(nrow, function(i) {
    cur_mu <- fitted.values
    obs <- residuals + cur_mu
    cur_summary <- obj$fits[[which_model]]$summary
    cur_var <- smooth_sigma_sq_pmax + sigma_q_sq
    sum(dnorm(obs, mean = cur_mu, sd = sqrt(cur_var),log = TRUE)) ## normal distribution!!
})

… with test stats as \[2 \times (l_{full}- l_{null})\] P-value is from pchisq and q-value from p.adjust(pval, method = "BH").