Sleuth code

Previously, we have derived transcript abundances from Kallisto for the samples (1, 3, 6, 13, 15, 18).
The conditions to compare are:

• Medium (1,13) vs 0.2MOI (6,18)

• Medium (1,13) vs polyl25 (3,15)

## Sleuth

# Specify the base directory of the results
base_dir <- "/pica/v3/g2016025/G2_proj/data"

#  Retrieve the list of samples
sample_id <- dir(file.path(base_dir,"kallisto.results"))

# Construct appropriate paths
 kal_dirs <- sapply(sample_id, function(id) file.path(base_dir, "kallisto.results", id))


# Create the experimental design matrix (here for the first comparison)
 s2c <- read.table(file.path(base_dir, "kallisto.results/kallisto.table.txt"), header = TRUE, stringsAsFactors=FALSE)

# Add the paths as an extra column
s2c <- dplyr::mutate(s2c, path = kal_dirs)

# Load the kallisto processed data into the object
so <- sleuth_prep(s2c, ~ condition)

# Estimate parameters for the general linear model (full)
so <- sleuth_fit(so)

# Estimate parameters for the sleuth reduced model (intersection only)
so <- sleuth_fit(so, ~1, 'reduced')

# Estimate the differential expression via the Wald test (in order to extract beta values, comparable
# to fold change)
so <- sleuth_wt(so,"MOI")

# Obtain results
results_table <- sleuth_results(so, 'MOI', test_type = 'wt')

## To acquire gene-centric results

# Construct mart object
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)

# Generate gene-centric sleuth model
so <- sleuth_prep(s2c, ~condition, target_mapping = t2g,aggregation_column = 'ens_gene')

#...Subsequent steps are the same as above