Merge pull request #117 from pachterlab/devel

bug fix patches

Author: pimentel

R/bootstrap.R

@@ -398,6 +398,7 @@ process_bootstrap <- function(i, samp_name, kal_path,
   if (read_bootstrap_tpm) {
     bs_quant_tpm <- aperm(apply(bs_mat, 1, counts_to_tpm,
                                 eff_len))
+    colnames(bs_quant_tpm) <- colnames(bs_mat)
 
     # gene level code is analogous here to below code
     if (gene_mode) {

R/measurement_error.R

@@ -202,10 +202,10 @@ sleuth_wt <- function(obj, which_beta, which_model = 'full') {
   if ( length(beta_i) == 0 ) {
     stop(paste0("'", which_beta,
         "' doesn't appear in your design. Try one of the following:\n",
-        colnames(d_matrix)))
+        paste(colnames(d_matrix), collapse = ' ')))
   } else if ( length(beta_i) > 1 ) {
     stop(paste0("Sorry. '", which_beta, "' is ambiguous for columns: ",
-        colnames(d_matrix[beta_i])))
+        paste(colnames(d_matrix[beta_i]), collapse = ' ')))
   }
 
   b <- sapply(obj$fits[[ which_model ]]$models,
@@ -237,8 +237,6 @@ sleuth_wt <- function(obj, which_beta, which_model = 'full') {
     qval = p.adjust(pval, method = 'BH')
     )
 
-  res <- dplyr::select(res, -x_group)
-
   obj <- add_test(obj, res, which_beta, 'wt', which_model)
 
   obj

R/sleuth.R

@@ -201,11 +201,19 @@ sleuth_prep <- function(
   msg('reading in kallisto results')
   sample_to_covariates$sample <- as.character(sample_to_covariates$sample)
 
+  if(nrow(sample_to_covariates) == 1 && !is.null(full_model)) {
+    warning("There is only one sample present, but you also provided a model. ",
+            "The model will be set to NULL to prevent downstream errors.\n",
+            "The sample can be viewed using sleuth_live after preparation, ",
+            "but you need more than one sample to run the other aspects of Sleuth.")
+    full_model <- NULL
+  }
+
   kal_dirs <- sample_to_covariates$path
   sample_to_covariates$path <- NULL
 
   msg('dropping unused factor levels')
-  samples_to_covariates <- droplevels(sample_to_covariates)
+  sample_to_covariates <- droplevels(sample_to_covariates)
 
   nsamp <- 0
   # append sample column to data
@@ -280,7 +288,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, , drop = FALSE])
 
     filter_df <- adf(target_id = names(filter_true))
 
@@ -298,7 +306,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, , drop = FALSE])
     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)
@@ -349,6 +357,7 @@ sleuth_prep <- function(
       # 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 <- data.table::data.table(target_mapping)
       target_mapping[target_mapping[[aggregation_column]] == "",
                      aggregation_column] <- NA
       agg_id <- unique(target_mapping[, aggregation_column, with = FALSE])
@@ -446,9 +455,10 @@ sleuth_prep <- function(
     })
 
     # 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.")
+    error_status <- sapply(bs_results, function(x) is(x, "try-error"))
+    if (any(error_status)) {
+      print(attributes(bs_results[error_status])$condition)
+      stop("At least one core from mclapply had an error. See the above error message(s) for more details.")
     }
 
     # mclapply is expected to retun the bootstraps in order; this is a sanity check of that
@@ -471,10 +481,10 @@ sleuth_prep <- function(
     # 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, ]
+      obs_counts <- obs_to_matrix(ret, "scaled_reads_per_base")[which_agg_id, , drop = FALSE]
     } else {
       names(sigma_q_sq) <- which_target_id
-      obs_counts <- obs_to_matrix(ret, "est_counts")[which_target_id, ]
+      obs_counts <- obs_to_matrix(ret, "est_counts")[which_target_id, , drop = FALSE]
     }
 
     sigma_q_sq <- sigma_q_sq[order(names(sigma_q_sq))]
@@ -560,7 +570,7 @@ check_target_mapping <- function(t_id, target_mapping) {
 #' @export
 norm_factors <- function(mat) {
   nz <- apply(mat, 1, function(row) !any(round(row) == 0))
-  mat_nz <- mat[nz, ]
+  mat_nz <- mat[nz, , drop = FALSE]
   p <- ncol(mat)
   geo_means <- exp(apply(mat_nz, 1, function(row) mean(log(row))))
   s <- sweep(mat_nz, 1, geo_means, `/`)
@@ -716,7 +726,7 @@ obs_to_matrix <- function(obj, value_name) {
   rownames(obs_counts) <- obs_counts$target_id
   obs_counts$target_id <- NULL
   obs_counts <- as.matrix(obs_counts)
-  obs_counts <- obs_counts[, obj$sample_to_covariates$sample]
+  obs_counts <- obs_counts[, obj$sample_to_covariates$sample, drop = FALSE]
 
   obs_counts
 }

tests/testthat/helper-setup.R

@@ -10,6 +10,12 @@ target_mapping <- read.table(file.path(data_path, 'target_mappings.txt'), header
 incomplete_mapping <- read.table(file.path(data_path, 'target_mappings_incomplete.txt'), header = TRUE,
   stringsAsFactors = FALSE, sep="\t", quote="")
 
+small_study_map <- data.frame(sample = "small_sample", condition = "test",
+                                path = "small_test_data/kallisto.N",
+                                stringsAsFactors = F)
+small_target_mapping <- read.table('small_test_data/target_mapping.txt', header = TRUE,
+  stringsAsFactors = FALSE, sep="\t", quote="")
+
 study_mapping <- read.table(file.path(data_path, 'study_design.txt'), header = TRUE,
   stringsAsFactors = FALSE)
 study_mapping <- dplyr::select(study_mapping, sample = run, condition)

tests/testthat/small_test_data/kallisto.N/abundance.h5

@@ -0,0 +1,16 @@
+target_id	gene_name
+NM_001168316	UGT3A2
+NM_174914	UGT3A2
+NR_031764	UGT3A2
+NM_004503	HOXC6
+NM_006897	HOXC9
+NM_014212	HOXC11
+NM_014620	HOXC4
+NM_017409	HOXC10
+NM_017410	HOXC13
+NM_018953	HOXC5
+NM_022658	HOXC8
+NM_153633	HOXC4
+NM_153693	HOXC6
+NM_173860	HOXC12
+NR_003084	HOXC5

tests/testthat/small_test_data/abundance.h5 tests/testthat/small_test_data/kallisto/abundance.h5

@@ -55,3 +55,11 @@ test_that("gene level", {
                         target_mapping = incomplete_mapping,
                         aggregation_column = "gene_name"))
 })
+
+test_that(".N target mappings", {
+  expect_warning(result.N <- sleuth_prep(small_study_map,
+                          target_mapping = small_target_mapping))
+  expect_warning(result.N <- sleuth_prep(small_study_map,
+                          target_mapping = small_target_mapping,
+                          aggregation_column = "gene_name"))
+})

tests/testthat/small_test_data/abundance.tsv tests/testthat/small_test_data/kallisto/abundance.tsv

@@ -1,7 +1,7 @@
 context("reading")
 
 test_that("get kallisto path", {
-  dir_name <- "small_test_data"
+  dir_name <- "small_test_data/kallisto"
 
   # the standard case
   file_name <- file.path(dir_name, "abundance.h5")
@@ -28,7 +28,7 @@ test_that("get kallisto path", {
 })
 
 test_that("both read types", {
-  dir_name <- "small_test_data"
+  dir_name <- "small_test_data/kallisto"
 
   h5_file_name <- file.path(dir_name, "abundance.h5")
   kal_h5 <- read_kallisto_h5(h5_file_name, read_bootstrap = FALSE)
@@ -42,7 +42,7 @@ test_that("both read types", {
 })
 
 test_that("generalized read", {
-  dir_name <- "small_test_data"
+  dir_name <- "small_test_data/kallisto"
 
   kal_dir <- read_kallisto(dir_name, read_bootstrap = TRUE)
   h5_file_name <- file.path(dir_name, "abundance.h5")