diff --git a/DESCRIPTION b/DESCRIPTION
index 76c4ef0..1304bf7 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,9 @@
Package: sleuth
Title: Tools for investigating RNA-Seq
-Version: 0.28.0
+Version: 0.28.1
Authors@R: c(person("Harold", "Pimentel", , "haroldpimentel@gmail.com", role = c("aut", "cre")))
-Description: Investigate transcript abundance from "kallisto" and differential expression analysis from RNA-Seq data.
+Description: Investigate transcript abundance from "kallisto" and differential
+ expression analysis from RNA-Seq data.
License: GPL-3
LazyData: true
URL: https://github.com/pachterlab/sleuth
@@ -17,7 +18,9 @@ Imports:
tidyr,
reshape2,
rhdf5,
+ parallel,
lazyeval,
+ matrixStats,
shiny
Suggests:
MASS,
@@ -25,3 +28,4 @@ Suggests:
testthat,
knitr
VignetteBuilder: knitr
+RoxygenNote: 5.0.1
diff --git a/NAMESPACE b/NAMESPACE
index 9e1e74b..b2781e9 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,5 +1,6 @@
# Generated by roxygen2 (4.1.1): do not edit by hand
+S3method("$<-",sleuth)
S3method(bias_table,kallisto)
S3method(bias_table,sleuth)
S3method(get_bootstraps,kallisto)
@@ -14,6 +15,7 @@ S3method(print,sleuth)
S3method(print,sleuth_model)
S3method(summary,sleuth)
S3method(tests,sleuth)
+export("transform_fxn<-")
export(basic_filter)
export(bias_table)
export(bs_sigma_summary)
@@ -25,6 +27,7 @@ export(get_bootstrap_summary)
export(get_bootstraps)
export(get_quantile)
export(kallisto_table)
+export(log_transform)
export(melt_bootstrap_sleuth)
export(models)
export(norm_factors)
@@ -63,6 +66,9 @@ export(sliding_window_grouping)
export(tests)
export(tpm_to_alpha)
export(transcripts_from_gene)
+export(transform_status)
+export(transform_status.sleuth)
+export(transform_status.sleuth_model)
import(dplyr)
importFrom(data.table,fread)
importFrom(lazyeval,interp)
diff --git a/R/bootstrap.R b/R/bootstrap.R
index 793362e..7bf8b92 100644
--- a/R/bootstrap.R
+++ b/R/bootstrap.R
@@ -145,7 +145,7 @@ aggregate_bootstrap <- function(kal, mapping, split_by = "gene_id",
if ( any(!complete.cases(mapping)) ) {
warning("Found some NAs in mapping. Removing them.")
- mapping <- mapping[complete.cases(mapping),]
+ mapping <- mapping[complete.cases(mapping), ]
}
m_bs <- melt_bootstrap(kal, column)
@@ -316,8 +316,8 @@ sample_bootstrap <- function(obj, n_samples = 100L) {
# matrix sample
for (s in 1:n_samples) {
for (idx in 1:nrow(which_samp)) {
- b <- which_samp[idx,s]
- sample_mat[[s]][,idx] <- obj$kal[[idx]]$bootstrap[[b]]$est_counts
+ b <- which_samp[idx, s]
+ sample_mat[[s]][, idx] <- obj$kal[[idx]]$bootstrap[[b]]$est_counts
}
}
@@ -359,9 +359,143 @@ dcast_bootstrap.kallisto <- function(obj, units, nsamples = NULL) {
mat <- matrix(NA_real_, nrow = n_features, ncol = length(which_bs))
for (j in seq_along(which_bs)) {
- mat[ ,j] <- obj[[ "bootstrap" ]][[which_bs[j]]][[ units ]]
+ mat[, j] <- obj[[ "bootstrap" ]][[which_bs[j]]][[ units ]]
}
rownames(mat) <- obj[["bootstrap"]][[1]][["target_id"]]
mat
}
+
+# Function to process bootstraps for parallelization
+process_bootstrap <- function(i, samp_name, kal_path,
+ num_transcripts, est_count_sf,
+ read_bootstrap_tpm, gene_mode,
+ extra_bootstrap_summary,
+ target_id, mappings, which_ids,
+ aggregation_column, transform_fxn)
+{
+ dot(i)
+ bs_quants <- list()
+
+ num_bootstrap <- as.integer(rhdf5::h5read(kal_path$path,
+ "aux/num_bootstrap"))
+ if (num_bootstrap == 0) {
+ stop(paste0("File ", kal_path, " has no bootstraps.",
+ "Please generate bootstraps using \"kallisto quant -b\"."))
+ }
+
+ # TODO: only perform operations on filtered transcripts
+ eff_len <- rhdf5::h5read(kal_path$path, "aux/eff_lengths")
+ bs_mat <- read_bootstrap_mat(fname = kal_path$path,
+ num_bootstraps = num_bootstrap,
+ num_transcripts = num_transcripts,
+ est_count_sf = est_count_sf)
+
+ if (read_bootstrap_tpm) {
+ bs_quant_tpm <- aperm(apply(bs_mat, 1, counts_to_tpm,
+ eff_len))
+
+ # gene level code is analogous here to below code
+ if (gene_mode) {
+ colnames(bs_quant_tpm) <- target_id
+ # Make bootstrap_num an explicit column; each is treated as a "sample"
+ bs_tpm_df <- data.frame(bootstrap_num = c(1:num_bootstrap),
+ bs_quant_tpm, check.names = F)
+ rm(bs_quant_tpm)
+ # Make long tidy table; this step is much faster
+ # using data.table melt rather than tidyr gather
+ tidy_tpm <- data.table::melt(bs_tpm_df, id.vars = "bootstrap_num",
+ variable.name = "target_id",
+ value.name = "tpm")
+ tidy_tpm <- data.table::as.data.table(tidy_tpm)
+ rm(bs_tpm_df)
+ tidy_tpm$target_id <- as.character(tidy_tpm$target_id)
+ tidy_tpm <- merge(tidy_tpm, mappings, by = "target_id",
+ all.x = T)
+ # Data.table dcast uses non-standard evaluation
+ # So quote the full casting formula to make sure
+ # "aggregation_column" is interpreted as a variable
+ # see: http://stackoverflow.com/a/31295592
+ quant_tpm_formula <- paste("bootstrap_num ~",
+ aggregation_column)
+ bs_quant_tpm <- data.table::dcast(tidy_tpm,
+ quant_tpm_formula, value.var = "tpm",
+ fun.aggregate = sum)
+ bs_quant_tpm <- as.matrix(bs_quant_tpm[, -1])
+ rm(tidy_tpm) # these tables are very large
+ }
+ bs_quant_tpm <- aperm(apply(bs_quant_tpm, 2,
+ quantile))
+ colnames(bs_quant_tpm) <- c("min", "lower", "mid",
+ "upper", "max")
+ ret$bs_quants[[samp_name]]$tpm <- bs_quant_tpm
+ }
+
+ if (gene_mode) {
+ # I can combine target_id and eff_len
+ # I assume the order is the same, since it's read from the same kallisto
+ # file and each kallisto file has the same order
+ eff_len_df <- data.frame(target_id, eff_len,
+ stringsAsFactors = F)
+ # make bootstrap number an explicit column to facilitate melting
+ bs_df <- data.frame(bootstrap_num = c(1:num_bootstrap),
+ bs_mat, check.names = F)
+ rm(bs_mat)
+ # data.table melt function is much faster than tidyr's gather function
+ # output is a long table with each bootstrap's value for each target_id
+ tidy_bs <- data.table::melt(bs_df, id.vars = "bootstrap_num",
+ variable.name = "target_id",
+ value.name = "est_counts")
+ rm(bs_df)
+ # not sure why, but the melt function always returns a factor,
+ # even when setting variable.factor = F, so I coerce target_id
+ tidy_bs$target_id <- as.character(tidy_bs$target_id)
+ # combine the long tidy table with eff_len and aggregation mappings
+ # note that bootstrap number is treated as "sample" here
+ # for backwards compatibility
+ tidy_bs <- dplyr::select(tidy_bs, target_id,
+ est_counts, sample = bootstrap_num)
+ tidy_bs <- merge(data.table::as.data.table(tidy_bs),
+ data.table::as.data.table(eff_len_df), by = "target_id",
+ all.x = T)
+ tidy_bs <- merge(tidy_bs, mappings, by = "target_id",
+ all.x = T)
+ # create the median effective length scaling factor for each gene
+ scale_factor <- tidy_bs[, scale_factor := median(eff_len),
+ by = eval(parse(text=aggregation_column))]
+ # use the old reads_per_base_transform method to get gene scaled counts
+ scaled_bs <- reads_per_base_transform(tidy_bs,
+ scale_factor$scale_factor,
+ aggregation_column,
+ mappings)
+ # this step undoes the tidying to get back a matrix format
+ # target_ids here are now the aggregation column ids
+ bs_mat <- data.table::dcast(scaled_bs, sample ~ target_id,
+ value.var = "scaled_reads_per_base")
+ # this now has the same format as the transcript matrix
+ # but it uses gene ids
+ bs_mat <- as.matrix(bs_mat[, -1])
+ rm(tidy_bs, scaled_bs)
+ }
+
+ if (extra_bootstrap_summary) {
+ bs_quant_est_counts <- aperm(apply(bs_mat, 2,
+ quantile))
+ colnames(bs_quant_est_counts) <- c("min", "lower",
+ "mid", "upper", "max")
+ ret$bs_quants[[samp_name]]$est_counts <- bs_quant_est_counts
+ }
+
+ bs_mat <- transform_fxn(bs_mat)
+ # If bs_mat was made at gene-level, already has column names
+ # If at transcript-level, need to add target_ids
+ if(!gene_mode) {
+ colnames(bs_mat) <- target_id
+ }
+ # all_sample_bootstrap[, i] bootstrap point estimate of the inferential
+ # variability in sample i
+ # NOTE: we are only keeping the ones that pass the filter
+ bootstrap_result <- matrixStats::colVars(bs_mat[, which_ids])
+
+ list(index = i, bs_quants = bs_quants, bootstrap_result = bootstrap_result)
+}
diff --git a/R/gene_analysis.R b/R/gene_analysis.R
new file mode 100644
index 0000000..df36a01
--- /dev/null
+++ b/R/gene_analysis.R
@@ -0,0 +1,14 @@
+propagate_transcript_filter <- function(filter_df, target_mapping,
+ grouping_column) {
+
+ filtered_target_mapping <- dplyr::inner_join(as.data.table(filter_df), # nolint
+ as.data.table(target_mapping), by = 'target_id') # nolint
+
+ filtered_target_mapping <- dplyr::select_(filtered_target_mapping,
+ grouping_column)
+
+ data.table::setnames(filtered_target_mapping, grouping_column, 'target_id')
+ filtered_target_mapping <- dplyr::distinct(filtered_target_mapping)
+
+ filtered_target_mapping
+}
diff --git a/R/likelihood.R b/R/likelihood.R
index d648346..845c62c 100644
--- a/R/likelihood.R
+++ b/R/likelihood.R
@@ -66,6 +66,16 @@ sleuth_lrt <- function(obj, null_model, alt_model) {
model_exists(obj, null_model)
model_exists(obj, alt_model)
+ if(!obj$fits[[alt_model]]$transform_synced) {
+ stop("Model '", alt_model, "' was not computed using the sleuth object's",
+ " current transform function. Please rerun sleuth_fit for this model.")
+ }
+
+ if(!obj$fits[[null_model]]$transform_synced) {
+ stop("Model '", null_model, "' was not computed using the sleuth object's",
+ " current transform function. Please rerun sleuth_fit for this model.")
+ }
+
if ( !likelihood_exists(obj, null_model) ) {
obj <- compute_likelihood(obj, null_model)
}
diff --git a/R/matrix.R b/R/matrix.R
index 36e9bb8..c53b9be 100644
--- a/R/matrix.R
+++ b/R/matrix.R
@@ -25,6 +25,11 @@
#' ("obs_norm" or "obs_raw")
#' @param which_units character vector of length one. Which units to use ("tpm"
#' or "est_counts")
+#' @return a \code{list} with an attribute 'data', which contains a matrix of target_ids
+#' and transcript expression in \code{which_units}
+#' @examples
+#' sleuth_matrix <- sleuth_to_matrix(sleuth_obj, 'obs_norm', 'tpm')
+#' head(sleuth_matrix$data) # look at first 5 transcripts, sorted by name
#' @export
sleuth_to_matrix <- function(obj, which_df, which_units) {
if ( !(which_df %in% c("obs_norm", "obs_raw")) ) {
diff --git a/R/measurement_error.R b/R/measurement_error.R
index af070f7..df3a0b8 100644
--- a/R/measurement_error.R
+++ b/R/measurement_error.R
@@ -30,16 +30,22 @@
#' \code{\link{models}} helpful.
#'
#' @param obj a \code{sleuth} object
-#' @param formula a formula specifying the design to fit OR a design matrix
-#' @param fit_name the name to store the fit in the sleuth
-#' object
+#' @param formula an R formula specifying the design to fit OR a design matrix.
+#' If you are interested in only fitting the model that was specified in \code{sleuth_prep}
+#' you do not need to specify it again (will be fit as the 'full' model).
+#' @param fit_name the name to store the fit in the sleuth object (at so$fits$fit_name).
+#' If \code{NULL}, the model will be named 'full'.
#' @param ... additional arguments passed to \code{sliding_window_grouping} and
#' \code{shrink_df}
-#' @return a sleuth object with updated attributes
+#' @return a sleuth object with updated attributes.
+#' @examples # If you specified the formula in sleuth_prep, you can simply run to run the full model
+#' so <- sleuth_fit(so)
+#' # The intercept only model can be fit like this
+#' so <- sleuth_fit(so, ~1, 'reduced')
#' @seealso \code{\link{models}} for seeing which models have been fit,
#' \code{\link{sleuth_prep}} for creating a sleuth object,
-#' \code{\link{sleuth_wt}} to test whether a coefficient is zero or
-#' \code{\link{sleuth_lrt}} to compare two different models
+#' \code{\link{sleuth_wt}} to test whether a coefficient is zero,
+#' \code{\link{sleuth_lrt}} to test nested models.
#' @export
sleuth_fit <- function(obj, formula = NULL, fit_name = NULL, ...) {
stopifnot( is(obj, 'sleuth') )
@@ -104,13 +110,12 @@ sleuth_fit <- function(obj, formula = NULL, fit_name = NULL, ...) {
l_smooth <- dplyr::mutate(l_smooth,
smooth_sigma_sq_pmax = pmax(smooth_sigma_sq, sigma_sq))
-
msg('computing variance of betas')
beta_covars <- lapply(1:nrow(l_smooth),
function(i) {
- row <- l_smooth[i,]
+ row <- l_smooth[i, ]
with(row,
- covar_beta(smooth_sigma_sq_pmax + sigma_q_sq, X, A)
+ covar_beta(smooth_sigma_sq_pmax + sigma_q_sq, X, A)
)
})
names(beta_covars) <- l_smooth$target_id
@@ -124,7 +129,8 @@ sleuth_fit <- function(obj, formula = NULL, fit_name = NULL, ...) {
summary = l_smooth,
beta_covars = beta_covars,
formula = formula,
- design_matrix = X)
+ design_matrix = X,
+ transform_synced = TRUE)
class(obj$fits[[fit_name]]) <- 'sleuth_model'
@@ -152,14 +158,24 @@ model_exists <- function(obj, which_model, fail = TRUE) {
#' every transcript.
#'
#' @param obj a \code{sleuth} object
-#' @param which_beta a character string of length one denoting which beta to
-#' test
+#' @param which_beta a character string of denoting which grouping to test.
+#' For example, if you have a model fit to 'treatment,' with values of neg_ctl, pos_ctl,
+#' and drug, you would need to run \code{sleuth_wt} once each for pos_ctl and drug
#' @param which_model a character string of length one denoting which model to
#' use
#' @return an updated sleuth object
+#' @examples # Assume we have a sleuth object with a model fit to both genotype and drug,
+#' models(so)
+#' # formula: ~genotype + drug
+#' # coefficients:
+#' # (Intercept)
+#' # genotypeKO
+#' # drugDMSO
+#' so <- sleuth_wt(so, 'genotypeKO')
+#' so <- sleuth_wt(so, 'drugDMSO')
#' @seealso \code{\link{models}} to view which models have been fit and which
#' coefficients can be tested, \code{\link{sleuth_results}} to get back
-#' a data.frame of the results
+#' a \code{data.frame} of the results
#' @export
sleuth_wt <- function(obj, which_beta, which_model = 'full') {
stopifnot( is(obj, 'sleuth') )
@@ -169,6 +185,11 @@ sleuth_wt <- function(obj, which_beta, which_model = 'full') {
substitute(obj), ") for a list of fitted models")
}
+ if(!obj$fits[[which_model]]$transform_synced) {
+ stop("Model '", which_model, "' was not computed using the sleuth object's",
+ " current transform function. Please rerun sleuth_fit for this model.")
+ }
+
d_matrix <- obj$fits[[which_model]]$design_matrix
# get the beta index
@@ -250,9 +271,10 @@ me_model_by_row <- function(obj, design, bs_summary) {
stopifnot( all.equal(names(bs_summary$sigma_q_sq), rownames(bs_summary$obs_counts)) )
stopifnot( length(bs_summary$sigma_q_sq) == nrow(bs_summary$obs_counts))
+# obs_counts <- obj$transform_fxn(bs_summary$obs_counts)
models <- lapply(1:nrow(bs_summary$obs_counts),
function(i) {
- me_model(design, bs_summary$obs_counts[i,], bs_summary$sigma_q_sq[i])
+ me_model(design, bs_summary$obs_counts[i, ], bs_summary$sigma_q_sq[i])
})
names(models) <- rownames(bs_summary$obs_counts)
@@ -288,7 +310,7 @@ me_heteroscedastic_by_row <- function(obj, design, samp_bs_summary, obs_counts)
models <- lapply(1:nrow(bs_summary$obs_counts),
function(i) {
- res <- me_white_model(design, obs_counts[i,], sigma_q_sq[i,], A)
+ res <- me_white_model(design, obs_counts[i, ], sigma_q_sq[i, ], A)
res$df$target_id <- rownames(obs_counts)[i]
res
})
@@ -338,8 +360,14 @@ me_white_var <- function(df, sigma_col, sigma_q_col, X, tXX_inv) {
res
}
+
+
#' @export
-bs_sigma_summary <- function(obj, transform = identity) {
+bs_sigma_summary <- function(obj, transform = identity, norm_by_length = FALSE) {
+ # if (norm_by_length) {
+ # scaling_factor <- get_scaling_factors(obj$obs_raw)
+ # reads_per_base_transform()
+ # }
obs_counts <- obs_to_matrix(obj, "est_counts")
obs_counts <- transform( obs_counts )
@@ -357,6 +385,117 @@ bs_sigma_summary <- function(obj, transform = identity) {
list(obs_counts = obs_counts, sigma_q_sq = bs_sigma)
}
+# transform reads into reads per base
+#
+#
+reads_per_base_transform <- function(reads_table, scale_factor_input,
+ collapse_column = NULL,
+ mapping = NULL,
+ norm_by_length = TRUE) {
+
+ reads_table <- data.table::as.data.table(reads_table)
+
+ if (is(scale_factor_input, 'data.frame')) {
+ # message('USING NORMALIZATION BY EFFECTIVE LENGTH')
+ # browser()
+ scale_factor_input <- data.table::as.data.table(dplyr::select(scale_factor_input, target_id,
+ sample, scale_factor))
+ reads_table <- merge(reads_table, scale_factor_input,
+ by = c('sample', 'target_id'), all.x=T)
+ } else {
+ reads_table[, scale_factor := scale_factor_input]
+ }
+ # browser()
+ reads_table[, reads_per_base := est_counts / eff_len]
+ reads_table[, scaled_reads_per_base := scale_factor * reads_per_base]
+
+ if (!is.null(collapse_column)) {
+ mapping <- data.table::as.data.table(mapping)
+ # old stuff
+ if (!(collapse_column %in% colnames(reads_table))) {
+ reads_table <- merge(reads_table, mapping, by = 'target_id', all.x=T)
+ }
+ # browser()
+ # reads_table <- dplyr::left_join(reads_table, mapping, by = 'target_id')
+
+ rows_to_remove <- !is.na(reads_table[[collapse_column]])
+ reads_table <- reads_table[rows_to_remove]
+ if ('sample' %in% colnames(reads_table)) {
+ reads_table <- reads_table[, j = list(scaled_reads_per_base = sum(scaled_reads_per_base)),
+ by = list(sample, eval(parse(text=collapse_column)))]
+ } else {
+ reads_table <- reads_table[, j = list(scaled_reads_per_base = sum(scaled_reads_per_base)),
+ by = eval(parse(text=collapse_column))]
+ }
+
+ data.table::setnames(reads_table, 'parse', 'target_id')
+ }
+
+ as_df(reads_table)
+}
+
+gene_summary <- function(obj, which_column, transform = identity,
+ norm_by_length = TRUE, num_cores=2) {
+ # stopifnot(is(obj, 'sleuth'))
+ msg(paste0('aggregating by column: ', which_column))
+ obj_mod <- obj
+ if (norm_by_length) {
+ tmp <- obj$obs_raw
+ # tmp <- as.data.table(tmp)
+ tmp <- dplyr::left_join(
+ data.table::as.data.table(tmp),
+ data.table::as.data.table(obj$target_mapping),
+ by = 'target_id')
+ tmp <- dplyr::group_by_(tmp, 'sample', which_column)
+ scale_factor <- dplyr::mutate(tmp, scale_factor = median(eff_len))
+ } else {
+ scale_factor <- median(obj_mod$obs_norm_filt$eff_len)
+ }
+ # scale_factor <- median(obj_mod$obs_norm_filt$eff_len)
+ obj_mod$obs_norm_filt <- reads_per_base_transform(obj_mod$obs_norm_filt,
+ scale_factor, which_column, obj$target_mapping, norm_by_length)
+ obj_mod$obs_norm <- reads_per_base_transform(obj_mod$obs_norm,
+ scale_factor, which_column, obj$target_mapping, norm_by_length)
+
+ obs_counts <- obs_to_matrix(obj_mod, "scaled_reads_per_base")
+ obs_counts <- transform(obs_counts)
+
+ # NEW CODE: Switched mclapply to be on second lapply
+ # This means that it's applied on the order of one bootstrap
+ # This creates a much smaller memory footprint
+ msg("starting mclapply process now")
+ obj_mod$kal <- lapply(seq_along(obj_mod$kal),
+ function(i) {
+ k <- obj_mod$kal[[i]]
+ current_sample <- obj_mod$sample_to_covariates$sample[i]
+ msg(paste('aggregating across sample: ', current_sample))
+ k$bootstrap <- parallel::mclapply(seq_along(k$bootstrap), function(j) {
+ b <- k$bootstrap[[j]]
+ b <- dplyr::mutate(b, sample = current_sample)
+ reads_per_base_transform(b, scale_factor, which_column,
+ obj$target_mapping, norm_by_length)
+ }, mc.cores = num_cores)
+
+ k
+ })
+
+ bs_summary <- sleuth_summarize_bootstrap_col(obj_mod, "scaled_reads_per_base",
+ transform)
+
+ bs_summary <- dplyr::group_by(bs_summary, target_id)
+ # FIXME: the column name 'bs_var_est_counts' is incorrect. should actually rename it above
+ bs_summary <- dplyr::summarise(bs_summary,
+ sigma_q_sq = mean(bs_var_scaled_reads_per_base))
+
+ bs_summary <- as_df(bs_summary)
+
+ bs_sigma <- bs_summary$sigma_q_sq
+ names(bs_sigma) <- bs_summary$target_id
+ bs_sigma <- bs_sigma[rownames(obs_counts)]
+
+ list(obs_counts = obs_counts, sigma_q_sq = bs_sigma)
+}
+
me_model <- function(X, y, sigma_q_sq) {
n <- nrow(X)
degrees_free <- n - ncol(X)
diff --git a/R/misc.R b/R/misc.R
index f840e78..37b802b 100644
--- a/R/misc.R
+++ b/R/misc.R
@@ -25,7 +25,7 @@ get_quantile <- function(data, which_col, lwr, upr, ignore_zeroes = TRUE) {
data$iqr <- valid
if (ignore_zeroes) {
- valid <- data[,which_col] > 0
+ valid <- data[, which_col] > 0
if (sum(valid) == 0) {
# in this case, everything in this bin is zero, so we just return the
# entire bin
@@ -36,7 +36,7 @@ get_quantile <- function(data, which_col, lwr, upr, ignore_zeroes = TRUE) {
data$iqr[!valid] <- FALSE
}
- data$iqr[valid] <- data[valid,which_col] %>%
+ data$iqr[valid] <- data[valid, which_col] %>%
ecdf(.)(.) %>%
ifelse(lwr <= . & . <= upr)
@@ -49,7 +49,7 @@ sliding_window_grouping <- function(data, x_col, y_col,
data <- as.data.frame(data)
- data <- mutate(data,x_ecdf = ecdf(data[,x_col])(data[,x_col]))
+ data <- mutate(data, x_ecdf = ecdf(data[, x_col])(data[, x_col]))
data <- mutate(data, x_group = cut(x_ecdf, n_bins))
data <- group_by(data, x_group)
@@ -63,7 +63,7 @@ sliding_window_grouping <- function(data, x_col, y_col,
shrink_df <- function(data, shrink_formula, filter_var) {
data <- as.data.frame(data)
s_formula <- substitute(shrink_formula)
- fit <- eval(loess(s_formula, data[data[,filter_var],]))
+ fit <- eval(loess(s_formula, data[data[, filter_var], ]))
data.frame(data, shrink = predict(fit, data))
}
@@ -119,17 +119,17 @@ apply_all_pairs <- function(mat, fun) {
all_pairs <- utils::combn(ids, 2)
for (i in 1:ncol(all_pairs)) {
- j <- all_pairs[1,i]
- k <- all_pairs[2,i]
+ j <- all_pairs[1, i]
+ k <- all_pairs[2, i]
- cur <- fun(mat[,j], mat[,k])
+ cur <- fun(mat[, j], mat[, k])
- res[j,k] <- cur
- res[k,j] <- cur
+ res[j, k] <- cur
+ res[k, j] <- cur
}
for (i in 1:length(ids)) {
- res[i,i] <- fun(mat[,i], mat[,i])
+ res[i, i] <- fun(mat[, i], mat[, i])
}
diff --git a/R/model.R b/R/model.R
index 7822afe..0c9edc7 100644
--- a/R/model.R
+++ b/R/model.R
@@ -33,9 +33,17 @@ print.sleuth_model <- function(obj) {
#' View which models have been fit
#'
-#' View which models have been fit
+#' @description View which models have been fit. sleuth fits data using R formulas
#'
-#' @param obj
+#' @param obj a sleuth object, containing kallisto results, usually made by sleuth_prep
+#' @return an R formula showing what has been fit
+#' @examples # imagine you have a set of samples from input and IP, and input has been set to intercept
+#' models(so)
+#' # [full]
+#' # formula: ~condition
+#' # coefficients:
+#' # (Intercept)
+#' # conditionIP
#' @export
models <- function(obj, ...) {
UseMethod('models')
@@ -46,7 +54,7 @@ models.sleuth <- function(obj, verbose = TRUE) {
# TODO: output a new in between models for readability
if (verbose) {
for (x in names(obj$fits)) {
- cat('[ ', x,' ]\n')
+ cat('[ ', x, ' ]\n')
models(obj$fits[[x]])
}
}
@@ -60,6 +68,40 @@ models.sleuth_model <- function(obj) {
print(obj)
}
+#' Check Transform Sync Status of Sleuth Fits
+#'
+#' This method prints out the sync status for all fits of \code{sleuth} object
+#' If the sleuth object's transform function was changed after sleuth_fit was used,
+#' the user will need to redo sleuth_fit for any fits already done.
+#'
+#' @param obj a \code{sleuth} object.
+#' @return a print out of each fit with the transform sync status.
+#' @export
+transform_status <- function(obj) {
+ useMethod('transform_status')
+}
+
+#' @export
+transform_status.sleuth <- function(obj, verbose=TRUE) {
+ if (is.null(obj$fits))
+ stop("sleuth obj has no fits.")
+
+ if (verbose) {
+ for (x in names(obj$fits)) {
+ cat('[ ', x, ' ]\n')
+ models(obj$fits[[x]]$transform_synced)
+ }
+ }
+
+
+ invisible(obj$fits)
+}
+
+#' @export
+transform_status.sleuth_model <- function(obj) {
+ print(obj$transform_synced)
+}
+
#' Extract design matrix
#'
#' Getter method for extracting a design matrix from a sleuth object
@@ -232,16 +274,31 @@ tests.sleuth <- function(obj, lrt = TRUE, wt = TRUE) {
#' This function extracts Wald test results from a sleuth object.
#'
#' @param obj a \code{sleuth} object
-#' @param test a character denoting the test to extract
+#' @param test a character string denoting the test to extract. Possible tests can be found by using \code{models(obj)}.
#' @param which_model a character string denoting the model. If extracting a wald test, use the model name. If extracting a likelihood ratio test, use 'lrt'.
#' @param rename_cols if \code{TRUE} will rename some columns to be shorter and
#' consistent with vignette
#' @param show_all if \code{TRUE} will show all transcripts (not only the ones
#' passing filters). The transcripts that do not pass filters will have
#' \code{NA} values in most columns.
-#' @return a \code{data.frame}
-#' @seealso \code{\link{sleuth_wt}}/\code{\link{sleuth_lrt}} to compute tests, \code{\link{models}} to
+#' @return a \code{data.frame} with the following columns:
+#' @return target_id: transcript name, e.g. "ENSXX#####" (dependent on the transcriptome used in kallisto)
+#' @return pval: p-value of the chosen model
+#' @return qval: false discovery rate adjusted p-value, using Benjamini-Hochberg (see \code{\link{p.adjust}})
+#' @return b: 'beta' value (effect size). Technically a biased estimator of the fold change
+#' @return se_b: standard error of the beta
+#' @return mean_obs: mean of natural log counts of observations
+#' @return var_obs: variance of observation
+#' @return tech_var: technical variance of observation from the bootstraps
+#' @return sigma_sq: raw estimator of the variance once the technical variance has been removed
+#' @return smooth_sigma_sq: smooth regression fit for the shrinkage estimation
+#' @return final_simga_sq: max(sigma_sq, smooth_sigma_sq); used for covariance estimation of beta
+#' @seealso \code{\link{sleuth_wt}} and \code{\link{sleuth_lrt}} to compute tests, \code{\link{models}} to
#' view which models, \code{\link{tests}} to view which tests were performed (and can be extracted)
+#' @examples
+#' models(sleuth_obj) # for this example, assume the formula is ~condition,
+#' and a coefficient is IP
+#' results_table <- sleuth_results(sleuth_obj, 'conditionIP')
#' @export
sleuth_results <- function(obj, test, test_type = 'wt',
which_model = 'full', rename_cols = TRUE, show_all = TRUE) {
@@ -292,7 +349,7 @@ sleuth_results <- function(obj, test, test_type = 'wt',
)
}
- if (show_all) {
+ if (show_all && !obj$gene_mode) {
tids <- adf(target_id = obj$kal[[1]]$abundance$target_id)
res <- dplyr::left_join(
data.table::as.data.table(tids),
@@ -301,12 +358,28 @@ sleuth_results <- function(obj, test, test_type = 'wt',
)
}
- if ( !is.null(obj$target_mapping) ) {
+ if ( !is.null(obj$target_mapping) && !obj$gene_mode) {
res <- dplyr::left_join(
data.table::as.data.table(res),
data.table::as.data.table(obj$target_mapping),
by = 'target_id')
}
+
+ if (!is.null(obj$target_mapping) && obj$gene_mode) {
+ # after removing the target_id column
+ # there are several redundant columns for each gene
+ # this line gets the unique line for each gene
+ target_mapping <- unique(dplyr::select(
+ obj$target_mapping,
+ -target_id))
+ # this line uses dplyr's "left_join" syntax for "by"
+ # to match "target_id" from the "res" table,
+ # and the gene_column from the target_mapping table.
+ res <- dplyr::left_join(data.table::as.data.table(res),
+ data.table::as.data.table(target_mapping),
+ by = c("target_id" = obj$gene_column))
+ }
+
res <- as_df(res)
dplyr::arrange(res, qval)
diff --git a/R/plots.R b/R/plots.R
index e60c431..6229f59 100644
--- a/R/plots.R
+++ b/R/plots.R
@@ -50,6 +50,11 @@ plot_mean_var <- function(obj,
substitute(obj), ')')
}
+ if(!obj$fits[[which_model]]$transform_synced) {
+ stop("Model '", which_model, "' was not computed using the sleuth object's",
+ " current transform function. Please rerun sleuth_fit for this model.")
+ }
+
df <- obj$fits[[which_model]]$summary
p <- ggplot(df, aes(mean_obs, sqrt(sqrt(sigma_sq_pmax))))
@@ -73,7 +78,7 @@ plot_mean_var <- function(obj,
#' @param pc_x integer denoting the principle component to use for the x-axis
#' @param pc_y integer denoting the principle component to use for the y-axis
#' @param use_filtered if TRUE, use filtered data. otherwise, use all data
-#' @param units either 'est_counts' or 'tpm'
+#' @param units either 'est_counts' ('scaled_reads_per_base' for gene_mode) or 'tpm'
#' @param text_labels if TRUE, use text labels instead of points
#' @param color_by a variable to color by. if NA, then will leave all as 'black'
#' @return a gpplot object
@@ -90,6 +95,16 @@ plot_pca <- function(obj,
...) {
stopifnot( is(obj, 'sleuth') )
+ if (obj$gene_mode & units == 'est_counts') {
+ warning(paste("your sleuth object is in gene mode,",
+ "but you selected 'est_counts'. Selecting 'scaled_reads_per_base'..."))
+ units <- 'scaled_reads_per_base'
+ } else if (!obj$gene_mode & units == 'scaled_reads_per_base') {
+ warning(paste("your sleuth object is not in gene mode,",
+ "but you selected 'scaled_reads_per_base'. Selecting 'est_counts'..."))
+ units <- 'scaled_reads_per_base'
+ }
+
mat <- NULL
if (use_filtered) {
mat <- spread_abundance_by(obj$obs_norm_filt, units,
@@ -99,9 +114,8 @@ plot_pca <- function(obj,
obj$sample_to_covariates$sample)
}
- pca_res <- prcomp(mat)
-
- pcs <- as_df(pca_res$rotation[, c(pc_x, pc_y)])
+ pca_res <- prcomp(t(mat))
+ pcs <- as_df(pca_res$x[, c(pc_x, pc_y)])
pcs$sample <- rownames(pcs)
rownames(pcs) <- NULL
@@ -136,7 +150,7 @@ plot_pca <- function(obj,
#' @param use_filtered if TRUE, use filtered data. otherwise, use all data
#' @param sample user input on which sample and which PC's contribute the most
#' @param PC principal component to view sample's contribution to that PC
-#' @param units either 'est_counts' or 'tpm'
+#' @param units either 'est_counts' ('scaled_reads_per_base' for gene_mode) or 'tpm'
#' @param pc_count # of PC's
#' @param scale scale or not
#' @param pca_loading_abs default true, to see all PC's magnitude (recommended)
@@ -155,6 +169,16 @@ plot_loadings <- function(obj,
stopifnot( is(obj, 'sleuth') )
#filtering?? doesn't work right now
+ if (obj$gene_mode & units == 'est_counts') {
+ warning(paste("your sleuth object is in gene mode,",
+ "but you selected 'est_counts'. Selecting 'scaled_reads_per_base'..."))
+ units <- 'scaled_reads_per_base'
+ } else if (!obj$gene_mode & units == 'scaled_reads_per_base') {
+ warning(paste("your sleuth object is not in gene mode,",
+ "but you selected 'scaled_reads_per_base'. Selecting 'est_counts'..."))
+ units <- 'scaled_reads_per_base'
+ }
+
# mat <- NULL
# if (use_filtered) {
# mat <- spread_abundance_by(obj$obs_norm_filt, units)
@@ -175,7 +199,7 @@ plot_loadings <- function(obj,
#given a sample
if (!is.null(sample)) {
toggle <- TRUE
- loadings <- pca_calc$x[sample,]
+ loadings <- pca_calc$x[sample, ]
if (pca_loading_abs) {
loadings <- abs(loadings)
loadings <- sort(loadings, decreasing = TRUE)
@@ -186,7 +210,7 @@ plot_loadings <- function(obj,
#given a PC, which samples contribute the most?
if (!toggle) {
- loadings <- pca_calc$x[,pc_input]
+ loadings <- pca_calc$x[, pc_input]
if (pca_loading_abs) {
loadings <- abs(loadings)
loadings <- sort(loadings, decreasing = TRUE)
@@ -235,7 +259,7 @@ plot_loadings <- function(obj,
#'
#' @param obj a \code{sleuth} object
#' @param use_filtered if TRUE, use filtered data. otherwise, use all data
-#' @param units either 'est_counts' or 'tpm'
+#' @param units either 'est_counts' ('scaled_reads_per_base' for gene_mode) or 'tpm'
#' @param pca_number user input on how many PC to display, otherwise default is 5
#' @param scale determines scaling
#' @param PC_relative gives the option to compare subsequent principal components and their contributions
@@ -249,6 +273,17 @@ plot_pc_variance <- function(obj,
PC_relative = NULL,
...) {
+ stopifnot( is(obj, 'sleuth') )
+ if (obj$gene_mode & units == 'est_counts') {
+ warning(paste("your sleuth object is in gene mode,",
+ "but you selected 'est_counts'. Selecting 'scaled_reads_per_base'..."))
+ units <- 'scaled_reads_per_base'
+ } else if (!obj$gene_mode & units == 'scaled_reads_per_base') {
+ warning(paste("your sleuth object is not in gene mode,",
+ "but you selected 'scaled_reads_per_base'. Selecting 'est_counts'..."))
+ units <- 'scaled_reads_per_base'
+ }
+
# mat <- NULL
# if (use_filtered) {
# mat <- spread_abundance_by(obj$obs_norm_filt, units)
@@ -257,7 +292,7 @@ plot_pc_variance <- function(obj,
# }
mat <- spread_abundance_by(obj$obs_norm_filt, units)
- pca_calc <- prcomp(mat, scale = scale) #PCA calculations
+ pca_calc <- prcomp(t(mat), scale. = scale) #PCA calculations
#computation
eigenvalues <- (pca_calc$sdev) ^ 2
@@ -273,20 +308,21 @@ plot_pc_variance <- function(obj,
}
pc_df <- data.frame(PC_count = 1:colsize, var = var_explained) #order here matters
- if(!is.null(PC_relative)) {
+ if (!is.null(PC_relative)) {
pc_df <- data.frame(PC_count = 1:length(eigenvalues), var = var_explained2)
- pc_df <- pc_df[PC_relative:nrow(pc_df),]
+ pc_df <- pc_df[PC_relative:nrow(pc_df), ]
if (!is.null(pca_number) && (PC_relative + pca_number <= length(eigenvalues))) {
- pc_df <- pc_df[1:pca_number,]
+ pc_df <- pc_df[1:pca_number, ]
} else if (PC_relative + 5 >= length(eigenvalues)) {
- pc_df <- pc_df[1:nrow(pc_df),]
+ pc_df <- pc_df[1:nrow(pc_df), ]
}
}
p <- ggplot(pc_df, aes(x = PC_count, y = var)) + geom_bar(stat = "identity")
p <- p + scale_x_continuous(breaks = 1:length(eigenvalues))
p <- p + ylab("% of variance") + xlab("principal components")
+ p <- p + ylim(0, 100)
p
@@ -299,7 +335,7 @@ plot_pc_variance <- function(obj,
#'
#' @param obj a \code{sleuth} object
#' @param use_filtered if TRUE, use filtered data. otherwise use all data
-#' @param units either 'est_counts' or 'tpm'
+#' @param units either 'est_counts' ('scaled_reads_per_base' for gene_mode) or 'tpm'
#' @param trans a string pointing to a function to use for the transformation.
#' @param grouping a string from the columns of \code{sample_to_covariates} in
#' the sleuth object for which to group and color by
@@ -315,6 +351,17 @@ plot_group_density <- function(obj,
offset = 1
) {
+ stopifnot( is(obj, 'sleuth') )
+ if (obj$gene_mode & units == 'est_counts') {
+ warning(paste("your sleuth object is in gene mode,",
+ "but you selected 'est_counts'. Selecting 'scaled_reads_per_base'..."))
+ units <- 'scaled_reads_per_base'
+ } else if (!obj$gene_mode & units == 'scaled_reads_per_base') {
+ warning(paste("your sleuth object is not in gene mode,",
+ "but you selected 'scaled_reads_per_base'. Selecting 'est_counts'..."))
+ units <- 'scaled_reads_per_base'
+ }
+
res <- kallisto_table(obj, use_filtered = use_filtered, include_covariates = TRUE)
# res <- NULL
# if (use_filtered) {
@@ -354,7 +401,7 @@ plot_group_density <- function(obj,
#' @param which_sample a character string matching a sample in
#' \code{obj$sample_to_covariates}
#' @param use_filtered if TRUE, use filtered data. Otherwise use all data
-#' @param units either \code{'est_counts'} or \code{'tpm'}
+#' @param units either \code{'est_counts'} (\code{'scaled_reads_per_base'} for gene_mode) or \code{'tpm'}
#' @param trans a string pointing to a function to use for the transformation.
#' @param offset the offset so that transformations such as log don't compute
#' -Inf. If NULL, then will not add an offset
@@ -368,6 +415,18 @@ plot_sample_density <- function(obj,
offset = 1
) {
res <- NULL
+
+ stopifnot( is(obj, 'sleuth') )
+ if (obj$gene_mode & units == 'est_counts') {
+ warning(paste("your sleuth object is in gene mode,",
+ "but you selected 'est_counts'. Selecting 'scaled_reads_per_base'..."))
+ units <- 'scaled_reads_per_base'
+ } else if (!obj$gene_mode & units == 'scaled_reads_per_base') {
+ warning(paste("your sleuth object is not in gene mode,",
+ "but you selected 'scaled_reads_per_base'. Selecting 'est_counts'..."))
+ units <- 'scaled_reads_per_base'
+ }
+
if (use_filtered) {
res <- obj$obs_norm_filt
} else {
@@ -398,7 +457,7 @@ plot_sample_density <- function(obj,
#' @param sample_x the string corresponding to the sample name in \code{obj$sample_to_covariates}
#' @param sample_y same as \code{sample_x} but for the y-axis
#' @param use_filtered if TRUE, use filtered data. otherwise, use all data
-#' @param units either 'est_counts' or 'tpm'
+#' @param units either 'est_counts' ('scaled_reads_per_base' for gene_mode) or 'tpm'
#' @param offset a linear offset to help deal with zeroes if transforming the abundances
#' @param point_alpha the alpha on the points
#' @param xy_line if TRUE, plot the xy_line
@@ -424,6 +483,17 @@ plot_scatter <- function(obj,
ylim = NULL
) {
+ stopifnot( is(obj, 'sleuth') )
+ if (obj$gene_mode & units == 'est_counts') {
+ warning(paste("your sleuth object is in gene mode,",
+ "but you selected 'est_counts'. Selecting 'scaled_reads_per_base'..."))
+ units <- 'scaled_reads_per_base'
+ } else if (!obj$gene_mode & units == 'scaled_reads_per_base') {
+ warning(paste("your sleuth object is not in gene mode,",
+ "but you selected 'scaled_reads_per_base'. Selecting 'est_counts'..."))
+ units <- 'scaled_reads_per_base'
+ }
+
abund <- NULL
if (use_filtered) {
abund <- spread_abundance_by(obj$obs_norm_filt, units,
@@ -438,19 +508,22 @@ plot_scatter <- function(obj,
abund <- dplyr::mutate(abund, target_id = rownames(abund))
if (!is.null(trans)) {
- sample_x <- paste0( trans, '( ', sample_x)
- sample_y <- paste0( trans, '( ', sample_y)
+ sample_x <- paste0( trans, '( `', sample_x)
+ sample_y <- paste0( trans, '( `', sample_y)
}
if ( (!is.null(offset) && !is.na(offset)) && offset != 0 ) {
off <- deparse(eval(offset))
- sample_x <- paste0(sample_x, ' + ', off)
- sample_y <- paste0(sample_y, ' + ', off)
+ sample_x <- paste0(sample_x, '` + ', off)
+ sample_y <- paste0(sample_y, '` + ', off)
}
- if (!is.null(trans)) {
+ if (!is.null(trans) & !is.null(offset)) {
sample_x <- paste0(sample_x, ' )')
sample_y <- paste0(sample_y, ' )')
+ } else {
+ sample_x <- paste0(sample_x, '` )')
+ sample_y <- paste0(sample_y, '` )')
}
p <- ggplot(abund, aes_string(sample_x, sample_y))
@@ -502,6 +575,11 @@ plot_vars <- function(obj,
) {
stopifnot( is(obj, 'sleuth') )
+ if(!obj$fits[[which_model]]$transform_synced) {
+ stop("Model '", which_model, "' was not computed using the sleuth object's",
+ " current transform function. Please rerun sleuth_fit for this model.")
+ }
+
cur_summary <- NULL
if (is.null(test)) {
@@ -610,7 +688,7 @@ plot_ma <- function(obj, test, test_type = 'wt', which_model = 'full',
#'
#' @param obj a sleuth object that contains a bootstrap summary (see \code{\link{get_bootstrap_summary}})
#' @param target_id a character vector of length 1 indicating the target_id (transcript or gene name depending on aggregation mode)
-#' @param units a character vector of either 'est_counts' or 'tpm'
+#' @param units a character vector of either 'est_counts' ('scaled_reads_per_base' for gene_mode) or 'tpm'
#' @param color_by a column in the sample to covariates to color by
#' @param x_axis_angle the angle of the x-axis labels
#' @return a ggplot2 object
@@ -622,6 +700,15 @@ plot_bootstrap <- function(obj,
x_axis_angle = 50
) {
stopifnot( is(obj, 'sleuth') )
+ if (obj$gene_mode & units == 'est_counts') {
+ warning(paste("your sleuth object is in gene mode,",
+ "but you selected 'est_counts'. Selecting 'scaled_reads_per_base'..."))
+ units <- 'scaled_reads_per_base'
+ } else if (!obj$gene_mode & units == 'scaled_reads_per_base') {
+ warning(paste("your sleuth object is not in gene mode,",
+ "but you selected 'scaled_reads_per_base'. Selecting 'est_counts'..."))
+ units <- 'scaled_reads_per_base'
+ }
df <- get_bootstrap_summary(obj, target_id, units)
@@ -870,7 +957,7 @@ plot_qq <- function(obj, test, test_type = 'wt', which_model = 'full',
#'
#' @param transcripts a vector of strings containing a list of transcripts to be plotted in a heatmap
#' @param obj a \code{sleuth} object
-#' @param units a string specifying which units to use, either tpm or est_counts
+#' @param units a string specifying which units to use, either tpm or est_counts (scaled_reads_per_base for gene_mode)
#' @param trans a string specifying a function to transform the data by
#' @return a \code{ggplot} object
#' @export
@@ -879,15 +966,26 @@ plot_transcript_heatmap <- function(obj,
units = 'tpm',
trans = 'log',
offset = 1) {
+
+ if (obj$gene_mode & units == 'est_counts') {
+ warning(paste("your sleuth object is in gene mode,",
+ "but you selected 'est_counts'. Selecting 'scaled_reads_per_base'..."))
+ units <- 'scaled_reads_per_base'
+ } else if (!obj$gene_mode & units == 'scaled_reads_per_base') {
+ warning(paste("your sleuth object is not in gene mode,",
+ "but you selected 'scaled_reads_per_base'. Selecting 'est_counts'..."))
+ units <- 'scaled_reads_per_base'
+ }
+
if(!all(transcripts %in% obj$obs_norm$target_id)) {
stop("Couldn't find the following transcripts: ",
paste(transcripts[!(transcripts %in% obj$obs_norm$target_id)], collapse = ", "),
"\n\tIt is highly likely that some of them were filtered out.")
}
- tabd_df <- obj$obs_norm[obj$obs_norm$target_id %in% transcripts,]
+ tabd_df <- obj$obs_norm[obj$obs_norm$target_id %in% transcripts, ]
- if(units == 'tpm') {
+ if (units == 'tpm') {
tabd_df <- dplyr::select(tabd_df, target_id, sample, tpm)
tabd_df <- reshape2::dcast(tabd_df, target_id ~sample, value.var = 'tpm')
} else if (units == 'est_counts') {
@@ -912,16 +1010,16 @@ plot_transcript_heatmap <- function(obj,
}
-#' Heatmap of expression
-#'
-#' Plot all of the points in an expression matrix
-#'
-#' @param exMat the expression matrix
-#' @param clustRows if TRUE, cluster the rows by hierarchical clustering.
-#' @param clustCols if TRUE, cluster the columns by hierarchical clustering.
-#' @param rowNames if TRUE, print the row names on the plot
-#' @param colNames if TRUE, print the column names on the plot
-#' @return a ggplot object
+# Heatmap of expression
+#
+# Plot all of the points in an expression matrix
+#
+# @param exMat the expression matrix
+# @param clustRows if TRUE, cluster the rows by hierarchical clustering.
+# @param clustCols if TRUE, cluster the columns by hierarchical clustering.
+# @param rowNames if TRUE, print the row names on the plot
+# @param colNames if TRUE, print the column names on the plot
+# @return a ggplot object
ggPlotExpression <- function(exMat, clustRows = TRUE, clustCols = TRUE,
rowNames = TRUE, colNames = TRUE) {
if (is(exMat, 'matrix')) {
@@ -953,25 +1051,25 @@ ggPlotExpression <- function(exMat, clustRows = TRUE, clustCols = TRUE,
legend.direction = 'vertical',
legend.position = 'top',
legend.background = element_rect(fill = "gray95", colour = "black", size = 0.5, linetype = 1),
- axis.title=element_blank())
+ axis.title = element_blank())
if (rowNames)
- p <- p + theme(axis.text.x=element_text(angle = 90, size=14))
+ p <- p + theme(axis.text.x = element_text(angle = 90, size = 14))
else
- p <- p + theme(axis.text.x=element_text(size=0))
+ p <- p + theme(axis.text.x = element_text(size = 0))
if (colNames)
- p <- p + theme(axis.text.y=element_text(size=14))
+ p <- p + theme(axis.text.y = element_text(size = 14))
else
- p <- p + theme(axis.text.y=element_text(size=0))
+ p <- p + theme(axis.text.y = element_text(size = 0))
p
#list(plot = p, rowOrder = rowOrder, colOrder = colOrder)
}
-#' Order by dendrogram
-#'
-#' @param mat a matrix where the rows are observations and the columns are different dimensions on the matrix
-#' @return a vector of label orderings
+# Order by dendrogram
+#
+# @param mat a matrix where the rows are observations and the columns are different dimensions on the matrix
+# @return a vector of label orderings
orderByDendrogram <- function(mat) {
hc <- hclust(dist(mat))
dc <- as.dendrogram(hc)
diff --git a/R/read_write.R b/R/read_write.R
index 4557bd2..1cc927e 100644
--- a/R/read_write.R
+++ b/R/read_write.R
@@ -115,7 +115,8 @@ read_kallisto_h5 <- function(fname, read_bootstrap = TRUE, max_bootstrap = NULL)
msg("Only reading ", max_bootstrap, " bootstrap samples")
num_bootstrap <- max_bootstrap
}
- bs_samples <- lapply(0:(num_bootstrap[1] - 1), function(i) {
+ bs_samples <- lapply(0:(num_bootstrap[1] - 1),
+ function(i) {
.read_bootstrap_hdf5(fname, i, abund)
})
} else {
@@ -165,10 +166,11 @@ read_bootstrap_mat <- function(fname,
num_bootstraps,
num_transcripts,
est_count_sf) {
- bs_mat <- matrix(nrow=num_bootstraps, ncol=num_transcripts)
- for (i in 1:nrow(bs_mat)) {
- bs_mat[i, ] <- rhdf5::h5read(fname, paste0("bootstrap/bs", i - 1)) / est_count_sf
+ bs_mat <- matrix(ncol = num_bootstraps, nrow = num_transcripts)
+ for (i in 1:ncol(bs_mat)) {
+ bs_mat[, i] <- rhdf5::h5read(fname, paste0("bootstrap/bs", i - 1)) / est_count_sf
}
+ bs_mat <- t(bs_mat)
target_id <- as.character(rhdf5::h5read(fname, "aux/ids"))
colnames(bs_mat) <- target_id
@@ -285,7 +287,7 @@ gtf_gene_names <- function(gtf_attr) {
for (i in 1:length(all_attr)) {
j <- 1
- while ((nchar(gene_id[i]) < 1 || nchar(trans_id[i]) < 1) &&
+ while ( (nchar(gene_id[i]) < 1 || nchar(trans_id[i]) < 1) &&
j <= length(all_attr[[i]]) ) {
if (all_attr[[i]][j] == "gene_id") {
gene_id[i] <- all_attr[[i]][j + 1] %>%
diff --git a/R/shiny.R b/R/shiny.R
index df0fc22..f81b8f6 100644
--- a/R/shiny.R
+++ b/R/shiny.R
@@ -150,7 +150,7 @@ sleuth_live <- function(obj, settings = sleuth_live_settings(),
textInput('hm_transcripts', label = 'enter target ids: ', value = '')
),
column(3,
- selectInput('hm_units', label = 'units:', choices = c('est_counts','tpm'), selected = 'tpm')
+ selectInput('hm_units', label = 'units:', choices = c('est_counts', 'tpm'), selected = 'tpm')
),
column(3,
textInput('hm_trans', label = 'tranform: ', value = 'log')
@@ -161,7 +161,7 @@ sleuth_live <- function(obj, settings = sleuth_live_settings(),
actionButton('hm_go', 'view')
)
),
- tags$style(type='text/css', "#hm_go {margin-top: 25px}"),
+ tags$style(type = 'text/css', "#hm_go {margin-top: 25px}"),
fluidRow(plotOutput('hm_plot')),
fluidRow(uiOutput("download_hm_plt_button"))
),
@@ -1241,15 +1241,15 @@ sleuth_live <- function(obj, settings = sleuth_live_settings(),
result
})
output$table_type <- renderUI({
- if(!is.null(obj$target_mapping)) {
+ if (!is.null(obj$target_mapping)) {
selectInput('pop_genes', label = 'table type: ',
- choices = list('transcript table' = 1, 'gene table' = 2),
+ choices = list('target table' = 1, 'gene table' = 2),
selected = 1)
}
})
output$group_by <- renderUI({
- if(!is.null(input$pop_genes) && input$pop_genes == 2) {
+ if (!is.null(input$pop_genes) && input$pop_genes == 2) {
selectInput('mappingGroup',
label = 'group by: ',
choices = names(obj$target_mapping)[2:length(names(obj$target_mapping))])
@@ -1264,7 +1264,7 @@ sleuth_live <- function(obj, settings = sleuth_live_settings(),
wb <- possible_tests[1]
}
- if(!is.null(input$mappingGroup) && (input$pop_genes == 2)) {
+ if (!is.null(input$mappingGroup) && (input$pop_genes == 2)) {
mg <- input$mappingGroup
test_table <- sleuth_gene_table(obj, wb, input$settings_test_type,
input$which_model_de, mg)
@@ -1297,7 +1297,7 @@ sleuth_live <- function(obj, settings = sleuth_live_settings(),
possible_tests <- list_tests(obj, input$settings_test_type)
which_test <- possible_tests[1]
}
- if(!is.null(input$mapping_group_lrt) && (input$pop_genes_lrt == 2)) {
+ if (!is.null(input$mapping_group_lrt) && (input$pop_genes_lrt == 2)) {
mg <- input$mapping_group_lrt
sleuth_gene_table(obj, which_test, input$settings_test_type,
which_group = mg)
@@ -1307,14 +1307,14 @@ sleuth_live <- function(obj, settings = sleuth_live_settings(),
})
output$table_type_lrt <- renderUI({
- if(!is.null(obj$target_mapping)) {
+ if (!is.null(obj$target_mapping)) {
selectInput('pop_genes_lrt', label = 'table type: ',
choices = list('transcript table' = 1, 'gene table' = 2),
selected = 1)
}
})
output$group_by_lrt <- renderUI({
- if(!is.null(input$pop_genes_lrt) && input$pop_genes_lrt == 2) {
+ if (!is.null(input$pop_genes_lrt) && input$pop_genes_lrt == 2) {
selectInput('mapping_group_lrt',
label = 'group by: ',
choices = names(obj$target_mapping)[2:length(names(obj$target_mapping))])
@@ -1351,13 +1351,13 @@ sleuth_live <- function(obj, settings = sleuth_live_settings(),
### Gene Viewer
# the name of the gene supplied
gv_var_text <- eventReactive(input$gv_go, {
- if(!is.null(obj$target_mapping)) {
+ if (!is.null(obj$target_mapping)) {
input$gv_var_input
}
})
output$gv_gene_column <- renderUI({
- if(!is.null(obj$target_mapping)) {
+ if (!is.null(obj$target_mapping)) {
selectInput('gv_gene_colname',
label = 'genes from: ',
choices = names(obj$target_mapping)[2:length(names(obj$target_mapping))])
@@ -1394,7 +1394,7 @@ sleuth_live <- function(obj, settings = sleuth_live_settings(),
output$gv_var_plts <- renderUI({
gv_plot_list <- lapply(1:input$gv_maxplots,
function(i) {
- gv_plotname <- paste("plot", i, sep="")
+ gv_plotname <- paste("plot", i, sep = "")
plotOutput(gv_plotname)
})
@@ -1404,9 +1404,9 @@ sleuth_live <- function(obj, settings = sleuth_live_settings(),
for (i in 1:15) {
local({
my_i <- i
- gv_plotname <- paste("plot", my_i, sep="")
+ gv_plotname <- paste("plot", my_i, sep = "")
output[[gv_plotname]] <- renderPlot({
- if(!is.null(obj$target_mapping) && !is.na(gv_var_list()[my_i])) {
+ if (!is.null(obj$target_mapping) && !is.na(gv_var_list()[my_i])) {
plot_bootstrap(obj,
gv_var_list()[my_i],
units = input$gv_var_units,
@@ -1417,7 +1417,7 @@ sleuth_live <- function(obj, settings = sleuth_live_settings(),
}
output$no_genes_message <- renderUI({
- if(is.null(obj$target_mapping)) {
+ if (is.null(obj$target_mapping)) {
HTML(' You need to add genes to your sleuth object to use the gene viewer.
',
' To add genes to your sleuth object, see the ',
'sleuth getting started guide.')
@@ -1432,7 +1432,7 @@ sleuth_live <- function(obj, settings = sleuth_live_settings(),
default_hm_plot_height <- 400
hm_plot_height <- function() {
- if(length(hm_transcripts()) > 5) {
+ if (length(hm_transcripts()) > 5) {
length(hm_transcripts()) * 60
} else {
default_hm_plot_height
@@ -1572,7 +1572,7 @@ enclosed_brush <- function(df, brush) {
xbool <- brush$xmin <= df[[xvar]] & df[[xvar]] <= brush$xmax
ybool <- brush$ymin <= df[[yvar]] & df[[yvar]] <= brush$ymax
- df[xbool & ybool,]
+ df[xbool & ybool, ]
}
#' settings for sleuth_live
diff --git a/R/sleuth.R b/R/sleuth.R
index dd1b0f6..4b19e9f 100644
--- a/R/sleuth.R
+++ b/R/sleuth.R
@@ -29,12 +29,24 @@ basic_filter <- function(row, min_reads = 5, min_prop = 0.47) {
mean(row >= min_reads) >= min_prop
}
+#' Natural log and offset transformation
+#'
+#' The default transformation function for converting the normalized counts.
+#'
+#' @param x numeric that must be >=0. represents an individual observed count (one transcript in one sample).
+#' @param offset numeric offset to prevent taking the log of 0.
+#' @return log(x + offset)
+#' @export
+log_transform <- function(x, offset=0.5) {
+ log(x + offset)
+}
+
# currently defunct
filter_df_all_groups <- function(df, fun, group_df, ...) {
grps <- setdiff(colnames(group_df), 'sample')
res <- sapply(unique(grps),
function(g) {
- apply(filter_df_by_groups(df, fun, group_df[,c('sample', g)], ...), 1, any)
+ apply(filter_df_by_groups(df, fun, group_df[, c('sample', g)], ...), 1, any)
})
vals <- apply(res, 1, any)
@@ -51,7 +63,7 @@ filter_df_by_groups <- function(df, fun, group_df, ...) {
function(g) {
valid_samps <- grps %in% g
valid_samps <- group_df$sample[valid_samps]
- apply(df[,valid_samps], 1, fun, ...)
+ apply(df[, valid_samps], 1, fun, ...)
})
}
@@ -62,15 +74,15 @@ filter_df_by_groups <- function(df, fun, group_df, ...) {
#' accounting for covariates, sequencing depth, technical and biological
#' variance.
#'
-#' @param sample_to_covariates is a \code{data.frame} which contains a mapping
+#' @param sample_to_covariates a \code{data.frame} which contains a mapping
#' from \code{sample} (a column) to some set of experimental conditions or
#' covariates. The column \code{path} is also required, which is a character
#' vector where each element points to the corresponding kallisto output directory. The column
#' \code{sample} should be in the same order as the corresponding entry in
#' \code{path}.
-#' @param full_model is a \code{formula} which explains the full model (design)
+#' @param full_model an R \code{formula} which explains the full model (design)
#' of the experiment OR a design matrix. It must be consistent with the data.frame supplied in
-#' \code{sample_to_covariates}.
+#' \code{sample_to_covariates}. You can fit multiple covariates by joining them with '+' (see example)
#' @param filter_fun the function to use when filtering.
#' @param target_mapping a \code{data.frame} that has at least one column
#' 'target_id' and others that denote the mapping for each target. if it is not
@@ -83,15 +95,26 @@ filter_df_by_groups <- function(df, fun, group_df, ...) {
#' @param ... additional arguments passed to the filter function
#' @param norm_fun_counts a function to perform between sample normalization on the estimated counts.
#' @param norm_fun_tpm a function to perform between sample normalization on the TPM
+#' @param aggregation_column a string of the column name in \code{\link{target_mapping}} to aggregate targets
#' @param read_bootstrap_tpm read and compute summary statistics on bootstraps on the TPM.
#' NOTE: Unnecessary for typical analyses
#' @param extra_bootstrap_summary if \code{TRUE}, compute extra summary
#' statistics needed for some plots (e.g. \code{\link{plot_bootstrap}}).
#' NOTE: Unnecessary for typical analyses
+#' @param transformation_function the transformation that should be applied
+#' to the normalized counts. Default is \code{'log(x+0.5)'} (i.e. natural log with 0.5 offset)
+#' NOTE: be sure you know what you're doing before you change this.
+#' @param num_cores an integer of the number of computer cores mclapply should use
+#' to speed up sleuth preparation
#' @return a \code{sleuth} object containing all kallisto samples, metadata,
#' and summary statistics
-#' @seealso \code{\link{sleuth_fit}} to fit a model, \code{\link{sleuth_wt}} to
-#' test whether a coeffient in the model is zero or \code{\link{sleuth_lrt}} to compare different models
+#' @examples # Assume we have run kallisto on a set of samples, and have two treatments,
+#' genotype and drug.
+#' colnames(s2c)
+#' # [1] "sample" "genotype" "drug" "path"
+#' so <- sleuth_prep(s2c, ~genotype + drug)
+#' @seealso \code{\link{sleuth_fit}} to fit a model, \code{\link{sleuth_wt}} or
+#' \code{\link{sleuth_lrt}} to perform hypothesis testing
#' @export
sleuth_prep <- function(
sample_to_covariates,
@@ -101,8 +124,12 @@ sleuth_prep <- function(
max_bootstrap = NULL,
norm_fun_counts = norm_factors,
norm_fun_tpm = norm_factors,
+ aggregation_column = NULL,
+ gene_mode = FALSE,
read_bootstrap_tpm = FALSE,
extra_bootstrap_summary = FALSE,
+ transformation_function = log_transform,
+ num_cores = 2L,
...) {
##############################
@@ -115,7 +142,7 @@ sleuth_prep <- function(
}
if (!is(full_model, "formula") && !is(full_model, "matrix")) {
- stop(paste0("'",substitute(full_model), "' (full_model) must be a formula or a matrix"))
+ stop(paste0("'", substitute(full_model), "' (full_model) must be a formula or a matrix"))
}
if (!("sample" %in% colnames(sample_to_covariates))) {
@@ -128,7 +155,7 @@ sleuth_prep <- function(
"' (sample_to_covariates) must contain a column named 'path'")
}
- if ( !is.null(target_mapping) && !is(target_mapping, 'data.frame')) {
+ if (!is.null(target_mapping) && !is(target_mapping, 'data.frame')) {
stop(paste0("'", substitute(target_mapping),
"' (target_mapping) must be a data.frame or NULL"))
} else if (is(target_mapping, 'data.frame')){
@@ -138,27 +165,37 @@ sleuth_prep <- function(
}
}
- if ( !is.null(max_bootstrap) && max_bootstrap <= 0 ) {
+ if (!is.null(max_bootstrap) && max_bootstrap <= 0 ) {
stop("max_bootstrap must be > 0")
}
- if ( any(is.na(sample_to_covariates)) ) {
+ if (any(is.na(sample_to_covariates))) {
warning("Your 'sample_to_covariance' data.frame contains NA values. This will likely cause issues later.")
}
- if ( is(full_model, "matrix") &&
- nrow(full_model) != nrow(sample_to_covariates) ) {
+ if (is(full_model, "matrix") &&
+ nrow(full_model) != nrow(sample_to_covariates)) {
stop("The design matrix number of rows are not equal to the number of rows in the sample_to_covariates argument.")
}
- if ( !is(norm_fun_counts, 'function') ) {
+ if (!is(norm_fun_counts, 'function')) {
stop("norm_fun_counts must be a function")
}
- if ( !is(norm_fun_tpm, 'function') ) {
+ if (!is(norm_fun_tpm, 'function')) {
stop("norm_fun_tpm must be a function")
}
+ if (!is.null(aggregation_column) && is.null(target_mapping)) {
+ stop(paste("You provided a 'aggregation_column' to aggregate by,",
+ "but not a 'target_mapping'. Please provided a 'target_mapping'."))
+ }
+
+ if (is.null(num_cores) || is.na(suppressWarnings(as.integer(num_cores))) ||
+ num_cores < 1 || num_cores > parallel::detectCores()) {
+ stop("num_cores must be an integer between 1 and the number of cores on your machine")
+ }
+
# TODO: ensure transcripts are in same order -- if not, report warning that
# kallisto index might be incorrect
@@ -170,6 +207,10 @@ sleuth_prep <- function(
kal_dirs <- sample_to_covariates$path
sample_to_covariates$path <- NULL
+
+ msg('dropping unused factor levels')
+ samples_to_covariates <- droplevels(sample_to_covariates)
+
nsamp <- 0
# append sample column to data
kal_list <- lapply(seq_along(kal_dirs),
@@ -177,7 +218,8 @@ sleuth_prep <- function(
nsamp <- dot(nsamp)
path <- kal_dirs[i]
suppressMessages({
- kal <- read_kallisto(path, read_bootstrap = FALSE, max_bootstrap = max_bootstrap)
+ kal <- read_kallisto(path, read_bootstrap = FALSE,
+ max_bootstrap = max_bootstrap)
})
kal$abundance <- dplyr::mutate(kal$abundance,
sample = sample_to_covariates$sample[i])
@@ -192,10 +234,10 @@ sleuth_prep <- function(
obs_raw <- dplyr::bind_rows(lapply(kal_list, function(k) k$abundance))
design_matrix <- NULL
- if ( is(full_model, 'formula') ) {
+ if (is(full_model, 'formula')) {
design_matrix <- model.matrix(full_model, sample_to_covariates)
} else {
- if ( is.null(colnames(full_model)) ) {
+ if (is.null(colnames(full_model))) {
stop("If matrix is supplied, column names must also be supplied.")
}
design_matrix <- full_model
@@ -203,6 +245,17 @@ sleuth_prep <- function(
rownames(design_matrix) <- sample_to_covariates$sample
obs_raw <- dplyr::arrange(obs_raw, target_id, sample)
+
+ ###
+ # try to deal with weird ensemble names
+ ###
+ if (!is.null(target_mapping)) {
+ tmp_names <- data.frame(target_id = kal_list[[1]]$abundance$target_id,
+ stringsAsFactors = FALSE)
+ target_mapping <- check_target_mapping(tmp_names, target_mapping)
+ rm(tmp_names)
+ }
+
ret <- list(
kal = kal_list,
kal_versions = kal_versions,
@@ -211,12 +264,16 @@ sleuth_prep <- function(
bootstrap_summary = NA,
full_formula = full_model,
design_matrix = design_matrix,
- target_mapping = target_mapping
+ target_mapping = target_mapping,
+ gene_mode = !is.null(aggregation_column),
+ gene_column = aggregation_column,
+ transform_fxn = transformation_function
)
# TODO: eventually factor this out
normalize <- TRUE
- if ( normalize ) {
+ if (normalize ) {
+
msg("normalizing est_counts")
est_counts_spread <- spread_abundance_by(obs_raw, "est_counts",
sample_to_covariates$sample)
@@ -224,7 +281,7 @@ sleuth_prep <- function(
filter_true <- filter_bool[filter_bool]
msg(paste0(sum(filter_bool), ' targets passed the filter'))
- est_counts_sf <- norm_fun_counts(est_counts_spread[filter_bool,])
+ est_counts_sf <- norm_fun_counts(est_counts_spread[filter_bool, ])
filter_df <- adf(target_id = names(filter_true))
@@ -242,7 +299,7 @@ sleuth_prep <- function(
msg("normalizing tpm")
tpm_spread <- spread_abundance_by(obs_raw, "tpm",
sample_to_covariates$sample)
- tpm_sf <- norm_fun_tpm(tpm_spread[filter_bool,])
+ tpm_sf <- norm_fun_tpm(tpm_spread[filter_bool, ])
tpm_norm <- as_df(t(t(tpm_spread) / tpm_sf))
tpm_norm$target_id <- rownames(tpm_norm)
tpm_norm <- tidyr::gather(tpm_norm, sample, tpm, -target_id)
@@ -253,12 +310,12 @@ sleuth_prep <- function(
obs_norm <- dplyr::arrange(obs_norm, target_id, sample)
tpm_norm <- dplyr::arrange(tpm_norm, target_id, sample)
- stopifnot( all.equal(obs_raw$target_id, obs_norm$target_id) &&
- all.equal(obs_raw$sample, obs_norm$sample) )
+ stopifnot(all.equal(obs_raw$target_id, obs_norm$target_id) &&
+ all.equal(obs_raw$sample, obs_norm$sample))
suppressWarnings({
- if ( !all.equal(dplyr::select(obs_norm, target_id, sample),
- dplyr::select(tpm_norm, target_id, sample)) ) {
+ if (!all.equal(dplyr::select(obs_norm, target_id, sample),
+ dplyr::select(tpm_norm, target_id, sample))) {
stop('Invalid column rows. In principle, can simply join. Please report error.')
}
@@ -270,13 +327,6 @@ sleuth_prep <- function(
# add in eff_len and len
obs_norm <- dplyr::bind_cols(obs_norm, dplyr::select(obs_raw, eff_len, len))
- msg("normalizing bootstrap samples")
-
- ret$kal <- lapply(seq_along(ret$kal), function(i) {
- normalize_bootstrap(ret$kal[[i]],
- tpm_size_factor = tpm_sf[i],
- est_counts_size_factor = est_counts_sf[i])
- })
obs_norm <- as_df(obs_norm)
ret$obs_norm <- obs_norm
@@ -286,6 +336,7 @@ sleuth_prep <- function(
ret$obs_norm_filt <- dplyr::semi_join(obs_norm, filter_df, by = 'target_id')
ret$tpm_sf <- tpm_sf
+ #### This code through the for loop is a candidate for moving to another function
path <- kal_dirs[1]
kal_path <- get_kallisto_path(path)
target_id <- as.character(rhdf5::h5read(kal_path$path, "aux/ids"))
@@ -293,62 +344,138 @@ sleuth_prep <- function(
ret$bs_quants <- list()
which_target_id <- ret$filter_df$target_id
- all_sample_bootstrap <- matrix(NA_real_,
- nrow = length(which_target_id),
- ncol = length(ret$kal))
- transformation_function <- function(x) log(x + 0.5)
+
+ if (ret$gene_mode) {
+ msg(paste0("aggregating by column: ", aggregation_column))
+ # Get list of IDs to aggregate on (usually genes)
+ # Also get the filtered list and update the "filter_df" and "filter_bool"
+ # variables for the sleuth object
+ target_mapping[target_mapping[[aggregation_column]] == "",
+ aggregation_column] <- NA
+ agg_id <- unique(target_mapping[, aggregation_column, with = F])
+ agg_id <- agg_id[[1]]
+ agg_id <- agg_id[!is.na(agg_id)]
+ mappings <- dplyr::select_(target_mapping, "target_id", aggregation_column)
+ mappings <- data.table::as.data.table(mappings)
+ which_tms <- which(mappings$target_id %in% which_target_id)
+ which_agg_id <- unique(mappings[which_tms, aggregation_column, with = F])
+ which_agg_id <- which_agg_id[[1]]
+ which_agg_id <- which_agg_id[!is.na(which_agg_id)]
+ filter_df <- adf(target_id = which_agg_id)
+ filter_bool <- agg_id %in% which_agg_id
+
+ msg(paste0(length(which_agg_id), " genes passed the filter"))
+
+ # Taken from gene_summary; scale normalized observed counts to "reads/base"
+ norm_by_length <- TRUE
+ tmp <- data.table::as.data.table(ret$obs_raw)
+ tmp <- merge(tmp, mappings,
+ by = "target_id", all.x=T)
+ scale_factor <- tmp[, scale_factor := median(eff_len),
+ by=list(sample,eval(parse(text=aggregation_column)))]
+ obs_norm_gene <- reads_per_base_transform(ret$obs_norm,
+ scale_factor, aggregation_column, mappings, norm_by_length)
+ # New code: get gene-level TPM (simple sum of normalized transcript TPM)
+ tmp <- data.table::as.data.table(tpm_norm)
+ tmp <- merge(tmp, mappings,
+ by = "target_id", all.x = T)
+ if (any(is.na(tmp[[aggregation_column]]))) {
+ rows_to_remove <- is.na(tmp[[aggregation_column]])
+ num_missing <- length(unique(tmp[rows_to_remove, target_id]))
+ warning(num_missing, " target_ids are missing annotations for the aggregation_column: ",
+ aggregation_column, ".\nThese target_ids will be dropped from the gene-level analysis.",
+ "\nIf you did not expect this, check your 'target_mapping' table for missing values.")
+ tmp <- tmp[!rows_to_remove]
+ }
+ tpm_norm_gene <- tmp[, j = list(tpm = sum(tpm)),
+ by = list(sample, eval(parse(text = aggregation_column)))]
+ data.table::setnames(tpm_norm_gene, 'parse', 'target_id')
+ tpm_norm_gene <- as_df(tpm_norm_gene)
+
+ # Same steps as above to add TPM column to "obs_norm" table
+ obs_norm_gene <- dplyr::arrange(obs_norm_gene, target_id, sample)
+ tpm_norm_gene <- dplyr::arrange(tpm_norm_gene, target_id, sample)
+
+ stopifnot(all.equal(dplyr::select(obs_norm_gene, target_id, sample),
+ dplyr::select(tpm_norm_gene, target_id, sample)))
+ suppressWarnings({
+ if ( !all.equal(dplyr::select(obs_norm_gene, target_id, sample),
+ dplyr::select(tpm_norm_gene, target_id, sample), check.attributes=F) ) {
+ stop('Invalid column rows. In principle, can simply join. Please report error.')
+ }
+
+ # obs_norm <- dplyr::left_join(obs_norm, data.table::as.data.table(tpm_norm),
+ # by = c('target_id', 'sample'))
+ obs_norm_gene <- dplyr::bind_cols(obs_norm_gene, dplyr::select(tpm_norm_gene, tpm))
+ })
+
+ # These are the updated gene-level variables
+ ret$filter_df <- adf(target_id = which_agg_id)
+ ret$filter_bool <- agg_id %in% which_agg_id
+ ret$obs_norm <- obs_norm_gene
+ ret$obs_norm_filt <- dplyr::semi_join(obs_norm_gene, filter_df, by = 'target_id')
+
+ rm(obs_norm, tpm_norm, obs_norm_gene, tpm_norm_gene)
+
+ # This is the gene-level version of the matrix
+ all_sample_bootstrap <- matrix(NA_real_,
+ nrow = length(which_agg_id),
+ ncol = length(ret$kal))
+ which_ids <- which_agg_id
+ } else {
+ all_sample_bootstrap <- matrix(NA_real_,
+ nrow = length(which_target_id),
+ ncol = length(ret$kal))
+ which_ids <- which_target_id
+ }
msg('summarizing bootstraps')
- for(i in 1:length(kal_dirs)) {
- dot(i)
+ bs_results <- parallel::mclapply(seq_along(kal_dirs), function(i) {
samp_name <- sample_to_covariates$sample[i]
kal_path <- get_kallisto_path(kal_dirs[i])
+ process_bootstrap(i, samp_name, kal_path,
+ num_transcripts, est_counts_sf[[i]],
+ read_bootstrap_tpm, ret$gene_mode,
+ extra_bootstrap_summary,
+ target_id, mappings, which_ids, ret$gene_column,
+ ret$transform_fxn)
+ }, mc.cores = num_cores)
+
+ # if mclapply results in an error (a warning is shown), then print error and stop
+ if (is(bs_results[[1]], "try-error")) {
+ print(attributes(bs_results[[1]])$condition)
+ stop("mclapply had an error. See the above error message for more details.")
+ }
- num_bootstrap <- as.integer(
- rhdf5::h5read(kal_path$path, "aux/num_bootstrap")
- )
- if (num_bootstrap == 0) {
- stop(paste0('File ', kal_path,
- ' has no bootstraps. Please generate bootstraps using "kallisto quant -b".')
- )
- }
+ # mclapply is expected to retun the bootstraps in order; this is a sanity check of that
+ indices <- sapply(bs_results, function(result) result$index)
+ stopifnot(identical(indices, order(indices)))
- # TODO: only perform operations on filtered transcripts
- eff_len <- rhdf5::h5read(kal_path$path, "aux/eff_lengths")
- bs_mat <- read_bootstrap_mat(fname = kal_path$path,
- num_bootstraps = num_bootstrap,
- num_transcripts = num_transcripts,
- est_count_sf = est_counts_sf[[i]])
-
- if (extra_bootstrap_summary) {
- bs_quant_est_counts <- aperm(apply(bs_mat, 2, quantile))
- colnames(bs_quant_est_counts) <- c("min", "lower", "mid", "upper",
- "max")
- ret$bs_quants[[samp_name]] <- list(est_counts = bs_quant_est_counts)
- }
+ if(read_bootstrap_tpm | extra_bootstrap_summary) {
+ ret$bs_quants <- lapply(bs_results, function(result) result$bs_quants)
+ names(ret$bs_quants) <- sample_to_covariates$sample
+ }
- if (read_bootstrap_tpm) {
- bs_quant_tpm <- aperm(apply(bs_mat, 1, counts_to_tpm, eff_len))
- bs_quant_tpm <- aperm(apply(bs_quant_tpm, 2, quantile))
- colnames(bs_quant_tpm) <- c("min", "lower", "mid", "upper", "max")
- ret$bs_quants[[samp_name]]$tpm <- bs_quant_tpm
- }
+ all_sample_bootstrap <- sapply(bs_results, function(result) result$bootstrap_result)
+ rownames(all_sample_bootstrap) <- which_ids
- bs_mat <- transformation_function(bs_mat)
- colnames(bs_mat) <- target_id
- # all_sample_bootstrap[, i] bootstrap point estimate of the inferential
- # variability in sample i
- # NOTE: we are only keeping the ones that pass the filter
- all_sample_bootstrap[, i] <- apply(bs_mat[, which_target_id], 2, var)
- } # end summarize bootstraps
+ # end summarize bootstraps
msg('')
sigma_q_sq <- rowMeans(all_sample_bootstrap)
- names(sigma_q_sq) <- which_target_id
- sigma_q_sq <- sigma_q_sq[order(names(sigma_q_sq))]
- obs_counts <- obs_to_matrix(ret, "est_counts")[which_target_id, ]
- obs_counts <- transformation_function(obs_counts)
+ # This is the rest of the gene_summary code
+ if (ret$gene_mode) {
+ names(sigma_q_sq) <- which_agg_id
+ obs_counts <- obs_to_matrix(ret, "scaled_reads_per_base")[which_agg_id, ]
+ } else {
+ names(sigma_q_sq) <- which_target_id
+ obs_counts <- obs_to_matrix(ret, "est_counts")[which_target_id, ]
+ }
+
+ sigma_q_sq <- sigma_q_sq[order(names(sigma_q_sq))]
+ obs_counts <- ret$transform_fxn(obs_counts)
+ obs_counts <- obs_counts[order(rownames(obs_counts)),]
ret$bs_summary <- list(obs_counts = obs_counts, sigma_q_sq = sigma_q_sq)
}
@@ -363,19 +490,62 @@ check_kal_pack <- function(kal_list) {
versions <- sapply(kal_list, function(x) attr(x, 'kallisto_version'))
u_versions <- unique(versions)
- if ( length(u_versions) > 1) {
- warning('More than one version of kallisto was used: ', u_versions)
+ if (length(u_versions) > 1) {
+ warning('More than one version of kallisto was used: ',
+ paste(u_versions, collapse = "\n"))
}
ntargs <- sapply(kal_list, function(x) attr(x, 'num_targets'))
u_ntargs <- unique(ntargs)
- if ( length(u_ntargs) > 1 ) {
+ if (length(u_ntargs) > 1 ) {
stop('Inconsistent number of transcripts. Please make sure you used the same index everywhere.')
}
list(versions = u_versions)
}
+# this function is mostly to deal with annoying ENSEMBL transcript names that
+# have a training .N to keep track of version number
+#
+# @return the target_mapping if an intersection is found. a target_mapping that
+# matches \code{t_id} if no matching is found
+check_target_mapping <- function(t_id, target_mapping) {
+ t_id <- data.table::as.data.table(t_id)
+ target_mapping <- data.table::as.data.table(target_mapping)
+
+ tmp_join <- dplyr::inner_join(t_id, target_mapping, by = 'target_id')
+
+ if (!nrow(tmp_join)) {
+ t_id <- dplyr::mutate(t_id, target_id_old = sub('\\.[0-9]+', '', target_id))
+ target_mapping_upd <- dplyr::rename(target_mapping, target_id_old = target_id)
+
+ tmp_join <- dplyr::inner_join(t_id, target_mapping_upd,
+ c('target_id_old'))
+
+ if (!nrow(tmp_join)) {
+ stop("couldn't solve nonzero intersection")
+ }
+
+ target_mapping <- dplyr::left_join(target_mapping,
+ unique(dplyr::select(tmp_join, target_id_new = target_id,
+ target_id = target_id_old),
+ by = NULL),
+ by = 'target_id')
+ # if we couldn't recover, use the old one
+ target_mapping <- dplyr::mutate(target_mapping,
+ target_id_new = ifelse(is.na(target_id_new), target_id, target_id_new))
+ target_mapping <- dplyr::select(target_mapping, -target_id)
+ target_mapping <- dplyr::rename(target_mapping, target_id = target_id_new)
+
+ warning(paste0('intersection between target_id from kallisto runs and ',
+ 'the target_mapping is empty. attempted to fix problem by removing .N ',
+ 'from target_id, then merging back into target_mapping.',
+ ' please check obj$target_mapping to ensure this new mapping is correct.'))
+ }
+
+ target_mapping
+}
+
#' Normalization factors
#'
#' More or less the DESeq size factors
@@ -386,9 +556,9 @@ check_kal_pack <- function(kal_list) {
#' @export
norm_factors <- function(mat) {
nz <- apply(mat, 1, function(row) !any(round(row) == 0))
- mat_nz <- mat[nz,]
+ mat_nz <- mat[nz, ]
p <- ncol(mat)
- geo_means <- exp(apply(mat_nz, 1, function(row) mean(log(row)) ))
+ geo_means <- exp(apply(mat_nz, 1, function(row) mean(log(row))))
s <- sweep(mat_nz, 1, geo_means, `/`)
sf <- apply(s, 2, median)
@@ -467,7 +637,7 @@ spread_abundance_by <- function(abund, var, which_order) {
result <- as.matrix(var_spread)
- result[, which_order]
+ result[, which_order, drop = FALSE]
}
#' @export
@@ -602,11 +772,23 @@ summary.sleuth <- function(obj, covariates = TRUE) {
#' @param test_type either 'wt' for Wald test or 'lrt' for likelihood ratio test
#' @param which_model a character string denoting which model to use
#' @param which_group a character string denoting which gene group to use
-#' @return a \code{data.frame} containing gene names, transcript names, and significance
+#' @return a \code{data.frame} with the following columns
+#'
+#' @return gene name; if ext_gene name specified, it will be legible gene name.
+#' If ens_gene name, it will be an Ensemble gene (assuming user followed the vignette)
+#' @return most_sig_trancript: Most significant transcript for the given gene
+#' @return pval: p-value for the test chosen
+#' @return qval: False discovery rate normalized p-value (Benjamini-Hochberg, see: \code{\link{p.adjust}})
+#' @return num_transcripts: Total number of transcripts for the given gene
+#' @return list_of_transcripts: All transcripts associated with this gene
+#' @examples sleuth_genes <- sleuth_gene_table(so, 'conditionIP', test_type ='wt',
+#' which_group = 'ext_gene')
+#' head(sleuth_genes) # show info for first 5 genes
+#' sleuth_genes[1:5, 6] # show transcripts for first 5 genes
#' @export
sleuth_gene_table <- function(obj, test, test_type = 'lrt', which_model = 'full', which_group = 'ens_gene') {
- if(is.null(obj$target_mapping)) {
+ if (is.null(obj$target_mapping)) {
stop("This sleuth object doesn't have added gene names.")
}
popped_gene_table <- sleuth_results(obj, test, test_type, which_model)
@@ -629,9 +811,9 @@ sleuth_gene_table <- function(obj, test, test_type = 'lrt', which_model = 'full'
filter_empty <- !filter_empty
popped_gene_table <- popped_gene_table[filter_empty, ]
- # popped_gene_table = popped_gene_table[!popped_gene_table[,1] == "",]
- # popped_gene_table = popped_gene_table[!is.na(popped_gene_table[,1]),] #gene_id
- # popped_gene_table = popped_gene_table[!is.na(popped_gene_table$qval),]
+ # popped_gene_table = popped_gene_table[!popped_gene_table[,1] == "", ]
+ # popped_gene_table = popped_gene_table[!is.na(popped_gene_table[,1]), ] #gene_id
+ # popped_gene_table = popped_gene_table[!is.na(popped_gene_table$qval), ]
popped_gene_table
}
@@ -653,11 +835,62 @@ sleuth_gene_table <- function(obj, test, test_type = 'lrt', which_model = 'full'
#' @export
transcripts_from_gene <- function(obj, test, test_type,
which_model, gene_colname, gene_name) {
+
+ # FIXME: this is a work around
+ obj$gene_mode <- FALSE
+
table <- sleuth_results(obj, test, test_type, which_model)
table <- dplyr::select_(table, ~target_id, gene_colname, ~qval)
table <- dplyr::arrange_(table, gene_colname, ~qval)
- if(!(gene_name %in% table[,2])) {
+ if (!(gene_name %in% table[, 2])) {
stop("Couldn't find gene ", gene_name)
}
- table$target_id[table[,2] == gene_name]
+ table$target_id[table[, 2] == gene_name]
+}
+
+#' Change sleuth transform function
+#'
+#' Replace the transformation function of a sleuth object
+#'
+#' NOTE: if you change the transformation function after having done a fit,
+#' the fit(s) will have to be redone using the new transformation.
+#' @examples transform_fxn(x) <- function(x) log2(x+0.5)
+#' @export
+`transform_fxn<-` <- function(obj, fxn) {
+ stopifnot(is.function(fxn))
+ obj$transform_fxn <- fxn
+ if(!is.null(obj$fits)) {
+ warning(paste("Your sleuth object has fits based on the old transform function.",
+ "Please rerun sleuth_prep and sleuth_fit."))
+ obj$fits <- lapply(obj$fits, function(x) {
+ x$transform_synced <- FALSE
+ x
+ })
+ }
+ obj
+}
+
+
+#' Extend internal '$<-' for sleuth object
+#'
+#' This extension is mainly to address case where
+#' transform_fxn is changed by user.
+#' This function informs user that the fits need to be redone
+#' and updates those fits.
+#' Otherwise it acts normally.
+#' @examples obj$transform_fxn <- function(x) log2(x+0.5)
+#' @export
+`$<-.sleuth` <- function(obj, name, value) {
+ obj[[name]] <- value
+ if(name=="transform_fxn") {
+ if(!is.null(obj$fits)) {
+ warning(paste("Your sleuth object has fits based on the old transform function.",
+ "Please rerun sleuth_prep and sleuth_fit."))
+ obj$fits <- lapply(obj$fits, function(x) {
+ x$transform_synced <- FALSE
+ x
+ })
+ }
+ }
+ obj
}
diff --git a/inst/doc/intro.html b/inst/doc/intro.html
index 9af0481..a4f6873 100644
--- a/inst/doc/intro.html
+++ b/inst/doc/intro.html
@@ -10,7 +10,7 @@
-
+
Introduction to sleuth
@@ -72,7 +72,7 @@
@@ -169,61 +169,104 @@ Example
## 4 ~/Downloads/cuffdiff2_data_kallisto_results/results/SRR493369/kallisto
## 5 ~/Downloads/cuffdiff2_data_kallisto_results/results/SRR493370/kallisto
## 6 ~/Downloads/cuffdiff2_data_kallisto_results/results/SRR493371/kallisto
-Now the “sleuth object” can be constructed. This requires three commands that (1) load the kallisto processed data into the object (2) estimate parameters for the sleuth response error measurement model and (3) perform differential analysis (testing). On a laptop the three steps should take about 2 minutes altogether.
+Now the “sleuth object” can be constructed. This requires four commands that (1) load the kallisto processed data into the object (2) estimate parameters for the sleuth response error measurement (full) model (3) estimate parameters for the sleuth reduced model, and (4) perform differential analysis (testing). On a laptop the four steps should take about a few minutes altogether.
First type
so <- sleuth_prep(s2c, ~ condition)
## reading in kallisto results
-## ......
-## normalizing est_counts
-## 50844 targets passed the filter
-## normalizing tpm
-## merging in metadata
-## normalizing bootstrap samples
-## summarizing bootstraps
-then
+## reading in kallisto results
+## ......
+##
+## normalizing est_counts
+## 50844 targets passed the filter
+## normalizing tpm
+## merging in metadata
+## normalizing bootstrap samples
+## summarizing bootstraps
+then fit the full model
-## shrinkage estimation
-## computing variance of betas
-and finally
-so <- sleuth_wt(so, 'conditionscramble')
In general, one can see the possible tests that could be performed using the which_beta parameter in sleuth_wt and examining the coefficients:
+## fitting measurement error models
+## shrinkage estimation
+## computing variance of betas
+Next, we fit the reduced model. In this case, the reduced model is the intercept-only model:
+so <- sleuth_fit(so, ~1, 'reduced')
## fitting measurement error models
+## shrinkage estimation
+## computing variance of betas
+and finally perform the test:
+so <- sleuth_lrt(so, 'reduced', 'full')
In general, we can test models that are nested using the likelihood ratio test. Viewing models which have been fit can be done using the models() function.
## [ full ]
## formula: ~condition
## coefficients:
## (Intercept)
-## conditionscramble
+## conditionscramble
+## [ reduced ]
+## formula: ~1
+## coefficients:
+## (Intercept)
+
+
Including gene names into transcript-level analysis
At this point the sleuth object constructed from the kallisto runs has information about the data, the experimental design, the kallisto estimates, the model fit, and the testing. In other words it contains the entire analysis of the data. There is, however, one piece of information that can be useful to add in, but that is optional. In reading the kallisto output sleuth has no information about genes, but this can be added allowing for searching and analysis by gene instead of transcript.
Since the example was constructed with the ENSEMBL human transcriptome, we will add gene names from ENSEMBL using biomaRt (there are other ways to do this as well):
First, install biomaRt with
source("http://bioconductor.org/biocLite.R")
biocLite("biomaRt")
Then collect gene names with
-
mart <- biomaRt::useMart(biomart = "ensembl", dataset = "hsapiens_gene_ensembl")
-
## Creating a generic function for 'nchar' from package 'base' in package 'S4Vectors'
+
mart <- biomaRt::useMart(biomart = "ENSEMBL_MART_ENSEMBL",
+ dataset = "hsapiens_gene_ensembl",
+ host = 'ensembl.org')
and add them into the sleuth table with
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)
so <- sleuth_prep(s2c, ~ condition, target_mapping = t2g)
-
## reading in kallisto results
-## ......
-## normalizing est_counts
-## 50844 targets passed the filter
-## normalizing tpm
-## merging in metadata
-## normalizing bootstrap samples
-## summarizing bootstraps
+
## reading in kallisto results
+
## ......
+
##
+
## normalizing est_counts
+
## 50844 targets passed the filter
+
## normalizing tpm
+
## merging in metadata
+
## normalizing bootstrap samples
+
## summarizing bootstraps
-
## shrinkage estimation
-## computing variance of betas
-
so <- sleuth_wt(so, which_beta = 'conditionscramble')
+
## fitting measurement error models
+
## shrinkage estimation
+
## computing variance of betas
+
so <- sleuth_fit(so, ~1, 'reduced')
+
## fitting measurement error models
+
## shrinkage estimation
+
## computing variance of betas
+
so <- sleuth_lrt(so, 'reduced', 'full')
This addition of metadata to transcript IDs is very general, and can be used to add in other information.
The best way to view the results is to generate the Shiny webpage that allows for exploratory data analysis:
To generate a table of results for analysis within R type
-
results_table <- sleuth_results(so, 'conditionscramble')
+
results_table <- sleuth_results(so, 'reduced:full', test_type = 'lrt')
+
+
Gene level analysis
+
Assuming biomaRt has been installed as in the previous step, sleuth can also be run in the ‘gene aggregation’ mode. In addition to requiring a target_mapping, a string which references a column to aggregate by, (aggregation_column). In our example, we could use ens_gene or ext_gene. It is preferable to use ens_gene as ext_gene tends to be a bit redundant.
+
The modified sleuth prep command looks as follows:
+
so <- sleuth_prep(s2c, ~condition, target_mapping = t2g,
+ aggregation_column = 'ens_gene')
+
## reading in kallisto results
+
## ......
+
##
+
## normalizing est_counts
+
## 50844 targets passed the filter
+
## normalizing tpm
+
## merging in metadata
+
## normalizing bootstrap samples
+
## summarizing bootstraps
+
## aggregating by column: ens_gene
+
## 14412 genes passed the filter
+
The remainder of the pipeline is unchanged. When running sleuth_live or sleuth_results, the gene column you selected will be listed as target_id, rather than the transcript name.
+
sleuth_prep might take a bit longer here because aggregation has to be done amongst all of the bootstraps as well. We are currently working to speed this up and expect to release it along with several memory improvements in the next version (0.28.2). One way to speed this up is to use more processes (assuming you have more cores). Internally, sleuth uses the function parallel::mclapply. You can set the number of cores as follows:
+
# set the number of available cores to 4
+options(mc.cores = 4L)
+