very simple sleuth/DESeq2 comparison at transcript level

sleuth_DESeq2.R

suppressPackageStartupMessages({
  library('Biobase')
  library('DESeq2')
  library('sleuth')
})

base_dir <- "~/Downloads/cuffdiff2_data_kallisto_results"

sample_id <- dir(file.path(base_dir,"results"))

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

s2c <- read.table(file.path(base_dir, "hiseq_info.txt"), header = TRUE, stringsAsFactors=FALSE)
s2c <- dplyr::select(s2c, sample = run_accession, condition)
s2c <- dplyr::mutate(s2c, path = kal_dirs)

mart <- biomaRt::useMart(biomart = "ENSEMBL_MART_ENSEMBL",
    dataset = "hsapiens_gene_ensembl",
    host = "may2015.archive.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)

suppressMessages(so <- sleuth_prep(s2c, ~ condition, target_mapping = t2g))

so <- sleuth_fit(so)

so <- sleuth_wt(so, 'conditionscramble')

sr <- sleuth_results(so, 'conditionscramble')
sr_sig <- dplyr::filter(sr, qval <= 0.01)

###
# DESeq2
###


obs_raw <- sleuth:::spread_abundance_by(so$obs_raw, "est_counts")
obs_raw <- round(obs_raw)
mode(obs_raw) <- 'integer'
s2c <- as.data.frame(s2c)
rownames(s2c) <- s2c$sample

cds <- ExpressionSet(obs_raw, AnnotatedDataFrame(s2c))
dds <- DESeqDataSetFromMatrix(exprs(cds), DataFrame(pData(cds)), ~ condition)
dds <- dds[ rowSums(counts(dds)) > 1, ]
dds <- DESeq(dds,quiet=TRUE)
ds_res <- results(dds)
ds_res$target_id <- rownames(ds_res)
ds_res <- as.data.frame(ds_res)
ds_sig <- dplyr::filter(ds_res, padj <= 0.01)

###
# comparison
###

nrow(sr_sig)
nrow(ds_sig)

length(intersect(sr_sig$target_id, ds_sig$target_id))