cole-trapnell-lab · DarwinAwardWinner · Jul 14, 2022 · Jul 14, 2022 · Jul 14, 2022 · Jul 14, 2022
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -66,4 +66,4 @@ biocViews: Sequencing, RNASeq, GeneExpression, DifferentialExpression,
     Infrastructure, DataImport, DataRepresentation, Visualization,
     Clustering, MultipleComparison, QualityControl
 Packaged: 2014-04-07 21:31:54 UTC; sarora
-RoxygenNote: 6.0.1
+RoxygenNote: 7.2.0
diff --git a/NAMESPACE b/NAMESPACE
@@ -88,12 +88,8 @@ importFrom(Biobase,fData)
 importFrom(Biobase,multiassign)
 importFrom(Biobase,pData)
 importFrom(BiocGenerics,"sizeFactors<-")
-importFrom(BiocGenerics,clusterApply)
-importFrom(BiocGenerics,clusterCall)
 importFrom(BiocGenerics,estimateDispersions)
 importFrom(BiocGenerics,estimateSizeFactors)
-importFrom(BiocGenerics,parCapply)
-importFrom(BiocGenerics,parRapply)
 importFrom(BiocGenerics,sizeFactors)
 importFrom(Matrix,readMM)
 importFrom(RANN,nn2)
@@ -146,10 +142,14 @@ importFrom(igraph,subcomponent)
 importFrom(igraph,vertex)
 importFrom(limma,removeBatchEffect)
 importFrom(matrixStats,rowSds)
+importFrom(parallel,clusterApply)
+importFrom(parallel,clusterCall)
 importFrom(parallel,detectCores)
 importFrom(parallel,makeCluster)
 importFrom(parallel,mclapply)
 importFrom(parallel,mcmapply)
+importFrom(parallel,parCapply)
+importFrom(parallel,parRapply)
 importFrom(parallel,splitIndices)
 importFrom(parallel,stopCluster)
 importFrom(plyr,.)
@@ -173,6 +173,7 @@ importFrom(stats,optim)
 importFrom(stats,p.adjust)
 importFrom(stats,prcomp)
 importFrom(stats,qnorm)
+importFrom(stats,quantile)
 importFrom(stats,rnorm)
 importFrom(stats,sd)
 importFrom(stats,setNames)
@@ -183,5 +184,6 @@ importFrom(stringr,str_split)
 importFrom(stringr,str_trim)
 importFrom(tibble,rownames_to_column)
 importFrom(utils,data)
+importFrom(utils,read.delim)
 importFrom(viridis,scale_color_viridis)
 useDynLib(monocle)
diff --git a/R/cds_conversion.R b/R/cds_conversion.R
@@ -116,11 +116,11 @@ exportCDS <- function(monocle_cds, export_to = c('Seurat', 'Scater'), export_all
 #' importCDS(scater_lung_all, import_all = T)
 #' }
 importCDS <- function(otherCDS, import_all = FALSE) {
-  if(class(otherCDS)[1] == 'seurat') {
+  if(is(otherCDS, 'seurat')) {
     requireNamespace("Seurat")
     data <- [email protected]
 
-    if(class(data) == "data.frame") {
+    if(is.data.frame(data)) {
       data <- as(as.matrix(data), "sparseMatrix")
     }
 
@@ -215,7 +215,7 @@ importCDS <- function(otherCDS, import_all = FALSE) {
     }
     monocle_cds@auxClusteringData$seurat <- mist_list
 
-  } else if (class(otherCDS)[1] == 'SCESet') {
+  } else if (is(otherCDS, 'SCESet')) {
     requireNamespace("scater")
 
     message('Converting the exprs data in log scale back to original scale ...')    

diff --git a/R/clustering.R b/R/clustering.R
@@ -72,6 +72,7 @@ clusterGenes<-function(expr_matrix, k, method=function(x){as.dist((1 - cor(Matri
 #' @return an updated CellDataSet object, in which phenoData contains values for Cluster for each cell
 #' @importFrom densityClust densityClust findClusters
 #' @importFrom igraph graph.data.frame cluster_louvain modularity membership
+#' @importFrom stats quantile
 #' @import ggplot2
 #' @importFrom RANN nn2
 #' @references Rodriguez, A., & Laio, A. (2014). Clustering by fast search and find of density peaks. Science, 344(6191), 1492-1496. doi:10.1126/science.1242072
@@ -377,4 +378,4 @@ clusterCells <- function(cds,
 #'   plot_cell_trajectory(cds, color_by = 'Time')
 #'   
 #'   return(list(new_cds = cds, ordering_genes = ordering_genes))
-#' }
+#' }
diff --git a/R/data_io.R b/R/data_io.R
@@ -32,6 +32,7 @@ get_genome_in_matrix_path <- function(matrix_path, genome=NULL) {
 #' @return a new CellDataSet object
 #' @export
 #' @importFrom Matrix readMM
+#' @importFrom utils read.delim
 #' @examples
 #' \dontrun{
 #' # Load from a Cell Ranger output directory
@@ -132,4 +133,4 @@ load_cellranger_data <- function(pipestance_path=NULL, genome=NULL, barcode_filt
                             lowerDetectionLimit=lowerDetectionLimit,
                             expressionFamily=expressionFamily)
   return(gbm)
-}
+}
diff --git a/R/differential_expression.R b/R/differential_expression.R
@@ -25,7 +25,7 @@ diff_test_helper <- function(x,
   x_orig <- x
   disp_guess <- 0
 
-  if (expressionFamily@vfamily %in% c("negbinomial", "negbinomial.size")){
+  if (any(expressionFamily@vfamily %in% c("negbinomial", "negbinomial.size"))){
     if (relative_expr == TRUE)
     {
       x <- x / Size_Factor
@@ -36,23 +36,23 @@ diff_test_helper <- function(x,
       if (is.null(disp_guess) == FALSE && disp_guess > 0 && is.na(disp_guess) == FALSE  ) {
         # FIXME: In theory, we could lose some user-provided parameters here
         # e.g. if users supply zero=NULL or something. 
-        if (expressionFamily@vfamily == "negbinomial")
+        if (any(expressionFamily@vfamily == "negbinomial"))
           expressionFamily <- negbinomial(isize=1/disp_guess)
         else
           expressionFamily <- negbinomial.size(size=1/disp_guess)
       }
     }
-  }else if (expressionFamily@vfamily %in% c("uninormal")){
+  }else if (any(expressionFamily@vfamily %in% c("uninormal"))){
     f_expression <- x
-  }else if (expressionFamily@vfamily %in% c("binomialff")){
+  }else if (any(expressionFamily@vfamily %in% c("binomialff"))){
     f_expression <- x
     #f_expression[f_expression > 1] <- 1
   }else{
     f_expression <- log10(x)
   }
 
   test_res <- tryCatch({
-    if (expressionFamily@vfamily %in% c("binomialff")){
+    if (any(expressionFamily@vfamily %in% c("binomialff"))){
       if (verbose){
         full_model_fit <- VGAM::vglm(as.formula(fullModelFormulaStr), epsilon=1e-1, family=expressionFamily)
         reduced_model_fit <- VGAM::vglm(as.formula(reducedModelFormulaStr), epsilon=1e-1, family=expressionFamily)                         
@@ -78,7 +78,7 @@ diff_test_helper <- function(x,
   error = function(e) { 
     if(verbose)
       print (e);
-      data.frame(status = "FAIL", family=expressionFamily@vfamily, pval=1.0, qval=1.0)
+      data.frame(status = "FAIL", family=expressionFamily@vfamily[1], pval=1.0, qval=1.0)
     #data.frame(status = "FAIL", pval=1.0) 
   }
   )
@@ -139,7 +139,7 @@ differentialGeneTest <- function(cds,
                                  ){
   status <- NA
 
-  if(class(cds)[1] != "CellDataSet") {
+  if(!is(cds, "CellDataSet")) {
     stop("Error cds is not of type 'CellDataSet'")
   }
 
@@ -155,7 +155,7 @@ differentialGeneTest <- function(cds,
   }
 
 
-  if (relative_expr && cds@expressionFamily@vfamily %in% c("negbinomial", "negbinomial.size")){
+  if (relative_expr && cds@any(expressionFamily@vfamily %in% c("negbinomial", "negbinomial.size"))){
     if (is.null(sizeFactors(cds)) || sum(is.na(sizeFactors(cds)))){
       stop("Error: to call this function with relative_expr==TRUE, you must first call estimateSizeFactors() on the CellDataSet.")
     }

diff --git a/R/expr_models.R b/R/expr_models.R
@@ -21,7 +21,7 @@ fit_model_helper <- function(x,
     orig_x <- x
     # FIXME: should we be using this here?
     # x <- x + pseudocount
-    if (expressionFamily@vfamily %in% c("negbinomial", "negbinomial.size")) {
+    if (any(expressionFamily@vfamily %in% c("negbinomial", "negbinomial.size"))) {
         if (relative_expr) {
             x <- x/Size_Factor
         }
@@ -32,14 +32,14 @@ fit_model_helper <- function(x,
             if (is.null(disp_guess) == FALSE && disp_guess >
                 0 && is.na(disp_guess) == FALSE) {
                 size_guess <- 1/disp_guess
-                if (expressionFamily@vfamily == "negbinomial")
+                if (any(expressionFamily@vfamily == "negbinomial"))
                   expressionFamily <- negbinomial(isize=1/disp_guess, ...)
                 else
                   expressionFamily <- negbinomial.size(size=1/disp_guess, ...)
             }
         }
     }
-    else if (expressionFamily@vfamily %in% c("uninormal", "binomialff")) {
+    else if (any(expressionFamily@vfamily %in% c("uninormal", "binomialff"))) {
         f_expression <- x
     }
     else {
@@ -61,12 +61,12 @@ fit_model_helper <- function(x,
         # what the user has specified for expression family
         #print(disp_guess)
         backup_expression_family <- NULL
-        if (expressionFamily@vfamily %in% c("negbinomial", "negbinomial.size")){
+        if (any(expressionFamily@vfamily %in% c("negbinomial", "negbinomial.size"))){
             disp_guess <- calculate_NB_dispersion_hint(disp_func, round(orig_x), expr_selection_func = max)
             backup_expression_family <- negbinomial()
-        }else if (expressionFamily@vfamily %in% c("uninormal")){
+        }else if (any(expressionFamily@vfamily %in% c("uninormal"))){
           backup_expression_family <- NULL
-        }else if (expressionFamily@vfamily %in% c("binomialff")){
+        }else if (any(expressionFamily@vfamily %in% c("binomialff"))){
           backup_expression_family <- NULL
         }else{
           backup_expression_family <- NULL
@@ -150,9 +150,9 @@ fitModel <- function(cds,
 responseMatrix <- function(models, newdata = NULL, response_type="response", cores = 1) {
     res_list <- mclapply(models, function(x) {
       if (is.null(x)) { NA } else {
-          if (x@family@vfamily %in% c("negbinomial", "negbinomial.size")) {
+          if (any(x@family@vfamily %in% c("negbinomial", "negbinomial.size"))) {
               predict(x, newdata = newdata, type = response_type)
-          } else if (x@family@vfamily %in% c("uninormal")) {
+          } else if (any(x@family@vfamily %in% c("uninormal"))) {
               predict(x, newdata = newdata, type = response_type)
           }
           else {
@@ -500,7 +500,7 @@ estimateDispersionsForCellDataSet <- function(cds, modelFormulaStr, relative_exp
   #                              expressionFamily=cds@expressionFamily)
   # }
 
-  if(!(('negbinomial' == cds@expressionFamily@vfamily) || ('negbinomial.size' == cds@expressionFamily@vfamily))){
+  if(!(any('negbinomial' == cds@expressionFamily@vfamily) || any('negbinomial.size' == cds@expressionFamily@vfamily))){
     stop("Error: estimateDispersions only works, and is only needed, when you're using a CellDataSet with a negbinomial or negbinomial.size expression family")
   }
 
@@ -511,7 +511,7 @@ estimateDispersionsForCellDataSet <- function(cds, modelFormulaStr, relative_exp
   options(dplyr.show_progress = T)
 
   # FIXME: this needs refactoring, badly.
-  if (cds@expressionFamily@vfamily %in% c("negbinomial", "negbinomial.size")){
+  if (any(cds@expressionFamily@vfamily %in% c("negbinomial", "negbinomial.size"))){
     if (length(model_terms) > 1 || (length(model_terms) == 1 && model_terms[1] != "1")){
       cds_pdata <- dplyr::group_by_(dplyr::select_(rownames_to_column(pData(cds)), "rowname", .dots=model_terms), .dots=model_terms)
       disp_table <- as.data.frame(cds_pdata %>% do(disp_calc_helper_NB(cds[,.$rowname], cds@expressionFamily, min_cells_detected)))

diff --git a/R/methods-CellDataSet.R b/R/methods-CellDataSet.R
@@ -95,7 +95,7 @@ function(object, modelFormulaStr="~ 1", relative_expr=TRUE, min_cells_detected=1
 #' @importFrom BiocGenerics sizeFactors
 checkSizeFactors <- function(cds)
 {
-  if (cds@expressionFamily@vfamily %in% c("negbinomial", "negbinomial.size"))
+  if (any(cds@expressionFamily@vfamily %in% c("negbinomial", "negbinomial.size")))
   {
     if (is.null(sizeFactors(cds))){
       stop("Error: you must call estimateSizeFactors() before calling this function.")

diff --git a/R/methods-CellTypeHierarchy.R b/R/methods-CellTypeHierarchy.R
@@ -313,11 +313,11 @@ classifyCells <- function(cds, cth, frequency_thresh=NULL, enrichment_thresh=NUL
 #' @export
 calculateMarkerSpecificity <- function(cds, cth, remove_ambig=TRUE, remove_unknown=TRUE){
 
-  if(class(cds)[1] != "CellDataSet") {
+  if(!is(cds, "CellDataSet")) {
     stop("Error cds is not of type 'CellDataSet'")
   }
 
-  if(class(cth)[1] != "CellTypeHierarchy") {
+  if(is(cth, "CellTypeHierarchy")) {
     stop("Error cth is not of type 'CellTypeHierarchy'")
   }
 
@@ -397,14 +397,14 @@ selectTopMarkers <- function(marker_specificities, num_markers = 10){
 #' @importFrom Biobase pData pData<-
 #' @export 
 markerDiffTable <- function (cds, cth, residualModelFormulaStr="~1", balanced=FALSE, reclassify_cells=TRUE, remove_ambig=TRUE, remove_unknown=TRUE, verbose=FALSE, cores=1) {
-  if(class(cds)[1] != "CellDataSet") {
+  if(!is(cds, "CellDataSet")) {
     stop("Error cds is not of type 'CellDataSet'")
   }
-  
-  if(class(cth)[1] != "CellTypeHierarchy") {
+
+  if(is(cth, "CellTypeHierarchy")) {
     stop("Error cth is not of type 'CellTypeHierarchy'")
   }
-  
+
   CellType <- NULL
   if (verbose)
     message("Classifying cells according to markers")

diff --git a/R/order_cells.R b/R/order_cells.R
@@ -1065,7 +1065,7 @@ orderCells <- function(cds,
                        num_paths = NULL,
                        reverse=NULL){
 
-  if(class(cds)[1] != "CellDataSet") {
+  if(!is(cds, "CellDataSet")) {
     stop("Error cds is not of type 'CellDataSet'")
   }
 
@@ -1200,7 +1200,7 @@ normalize_expr_data <- function(cds,
   }
 
   norm_method <- match.arg(norm_method)
-  if (cds@expressionFamily@vfamily %in% c("negbinomial", "negbinomial.size")) {
+  if (any(cds@expressionFamily@vfamily %in% c("negbinomial", "negbinomial.size"))) {
 
     # If we're going to be using log, and the user hasn't given us a pseudocount
     # set it to 1 by default.
@@ -1245,7 +1245,7 @@ normalize_expr_data <- function(cds,
       FM <- Matrix::t(Matrix::t(FM)/sizeFactors(cds))
       FM <- FM + pseudo_expr
     }
-  }else if (cds@expressionFamily@vfamily == "binomialff") {
+  }else if (any(cds@expressionFamily@vfamily == "binomialff")) {
     if (norm_method == "none"){
       #If this is binomial data, transform expression values into TF-IDF scores.
       ncounts <- FM > 0
@@ -1254,7 +1254,7 @@ normalize_expr_data <- function(cds,
     }else{
       stop("Error: the only normalization method supported with binomial data is 'none'")
     }
-  }else if (cds@expressionFamily@vfamily == "Tobit") {
+  }else if (any(cds@expressionFamily@vfamily == "Tobit")) {
     FM <- FM + pseudo_expr
     if (norm_method == "none"){
 
@@ -1263,7 +1263,7 @@ normalize_expr_data <- function(cds,
     }else{
       stop("Error: the only normalization methods supported with Tobit-distributed (e.g. FPKM/TPM) data are 'log' (recommended) or 'none'")
     }
-  }else if (cds@expressionFamily@vfamily == "uninormal") {
+  }else if (any(cds@expressionFamily@vfamily == "uninormal")) {
     if (norm_method == "none"){
       FM <- FM + pseudo_expr
     }else{
@@ -1415,7 +1415,7 @@ reduceDimension <- function(cds,
       #FM <- FM[genes_to_keep,]
       #expression_means <- expression_means[genes_to_keep]
       #expression_vars <- expression_vars[genes_to_keep]
-      # Here✬s how to take the top PCA loading genes, but using
+      # Here's how to take the top PCA loading genes, but using
       # sparseMatrix operations the whole time, using irlba.
 
 
@@ -1618,7 +1618,7 @@ project2MST <- function(cds, Projection_Method){
         projection <- rbind(projection, tmp)
         distance <- c(distance, dist(rbind(Z_i, tmp)))
       }
-      if(class(projection) != 'matrix')
+      if(!is.matrix(projection))
         projection <- as.matrix(projection)
       P[, i] <- projection[which(distance == min(distance))[1], ] #use only the first index to avoid assignment error
     }