We should implement the kind of reference-panel based test Darren used in the sci-ATAC-seq paper to find genes that are specifically expressed in each cell type.

I think there's basically never a scenario where we want to be able call partition_cells() without calling cluster_cells() first. I mean conceptually it can work, but it just introduces another bubble in an otherwise linear workflow.

expression_matrix = readRDS(url("http://jpacker-data.s3.amazonaws.com/for-cole/for-cole.ciliated.amphid.neurons.exprs.rds"))
cell_metadata = readRDS(url("http://jpacker-data.s3.amazonaws.com/for-cole/for-cole.ciliated.amphid.neurons.pData.rds"))
gene_annotation = readRDS(url("http://jpacker-data.s3.amazonaws.com/for-cole/for-cole.ciliated.amphid.neurons.fData.rds"))
cds <- new_cell_data_set(expression_matrix,
                         cell_metadata = cell_metadata,
                         gene_metadata = gene_annotation)


cds <- estimate_size_factors(cds)
cds <- preprocess_cds(cds, num_dim = 20)

preprocess_cds raises

Error in assays<-(*tmp*, value = new("SimpleList", listData = list( :
current and replacement dimnames() differ

at

assays(cds)$normalized_exprs <- FM

should check version too

Right now it has a thing called markers_linear, which is confusing and different than other methods that accept "norm_method"

Right now, assumes negative binomial. Replace with generic function to compute CV?

Many of our papers use labels for the root and the branch tips, but these were added in Illustrator. Adding them automatically would make trajectory plots look a lot nicer and it would be easier to call functions that require branch or path names as input (e.g. see #26)

Default policy will be to correct for testing each term across genes. So each term's p values will be adjusted separately.

Only the model fitting functions need to know about this, so we should simplify the class and new_cell_data_set

We should add a volcano plot function that takes tables from fit_models(), or maybe takes output from coefficient_table()

UMAP seems to make much more complex trajectories which need larger values of ncenter than DDRTree. We should re-tune the function that automatically selects the value of this parameter in learn_graph

It would be good to be able to plot cells after reduce_dimensions() but before cluster_cells(), partition_cells(), learn_graph(), order_cells(), etc. Right now if you try to plot after reduce_dimension but before cluster_cells(), you get an error. Maybe this should just be a message, and only then if you are trying to color_by="Cluster"?

If you have a cell_data_set object, calling pre_process_cds or reduce_dimension again after calling learn_graph should clear out old data (clusters, louvain components, etc) that is dependent on reduced dimensional coordinates.

cds <- reduce_dimension(cds)
plot_cells(cds, color_by = "Cluster")

cds <- cluster_cells(cds)
plot_cells(cds, color_by = "Cluster")

cds <- learn_graph(cds)
plot_cells(cds, color_by = "cell.type")

cds <- reduce_dimension(cds)
plot_cells(cds, color_by = "Cluster") # Shows old graph. Bad!

bc ugh

fit_models has checks to validate formula, but they are are too stringent, because they don't catch formula that include functions of columns in colData.

  if (length(model_form[[2]]) == 1) {
    if (!as.character(model_form[[2]]) %in% c(names(colData(cds)), "~", "1", "|", "+", "-", ":", "*", "^", "I")) {
      stop(paste(as.character(model_form[[2]][[i]]), "formula element is missing"))
    }
  } else {
    for(i in 1:length(model_form[[2]])) {
      if (!as.character(model_form[[2]][[i]]) %in% c(names(colData(cds)), "~", "1", "|", "+", "-", ":", "*", "^", "I")) {
        stop(paste(as.character(model_form[[2]][[i]]), "formula element is missing"))
      }
    }
  }

save normalizations in different slots

We're using matlab::eye which we could easily do without bringing in that dependency

It's really hard to see the violins, likely because whenever you have a group that's nearly all zero, ggplot's geom_violin's default policy scales them all to have the same width, which effectively turns the nonzero groups into thin lines.

https://cran.r-project.org/web/packages/uwot/index.html

We should have a short section on using Garnett to classify the L2 worm clusters on the Monocle site

Adding here as note requested by Cole -- on R 3.5.2 (the version of R that is linked against OpenBlas on the cluster), I am not able to run fit_models with anything other than cores = 1, regardless of the dataset. Even just for a small number of genes seems to hang at least in my hands. Here is a toy example that illustrates the problem on my end:

# Make random test matrix
library(Matrix)
library(monocle3)

# Make a fake cds object
genes = 100
cells = 500
mat = matrix(sample(0:1, genes * cells, replace = TRUE), genes, cells)
rownames(mat) = paste0('cell', 1:genes)
colnames(mat) = paste0('gene', 1:cells)

pdata = data.frame('cell'=colnames(mat))
rownames(pdata) = pdata$cell
pdata$category = FALSE
pdata$category[0:floor(cells/2)] = TRUE

fdata = data.frame('gene_short_name'=rownames(mat))
rownames(fdata) = fdata$gene_short_name

cds = monocle3::new_cell_data_set(expression_data = mat, cell_metadata=pdata, gene_metadata=fdata)

# Now try fit models
## Works fine
models = monocle3::fit_models(cds, model_formula_str='~category', cores=1, expression_family='binomial')

## Hangs indefinitely
models = monocle3::fit_models(cds, model_formula_str='~category', cores=2, expression_family='binomial')

Bob reports that plot_pc_variance_explained crashes his machine on an elegans dataset with 220k cells

cds <- load_a549()

set.seed(42)
cds <- preprocess_cds(cds)
cds <- reduce_dimension(cds, umap.fast_sgd=FALSE, cores=1)
cds <- partition_cells(cds)
cds <- learn_graph(cds)

#Negative control: randomize UMAP coordinates:
reducedDims(cds)$UMAP[,1] = reducedDims(cds)$UMAP[sample.int(nr),1]
reducedDims(cds)$UMAP[,2] = reducedDims(cds)$UMAP[sample.int(nr),2]

Calling principal_graph_test with neighbor_graph="principal_graph" returns lots of significant genes. Calling it with neighbor_graph="knn" returns none or a tiny number.

Users should have an easy way to pull out the cells along a lineage of interest. The old monocle3_alpha had a nice function for this: traverseGraphCDS().

In both 2D and 3D. Depends on #10

choose_graph_segments

We need to add/review functions needed to import data from Seurat (v3.0) objects.

By default, cluster_cells should automatically select a resolution parameter that finds the resolution that maximizes modularity while simultaneously minimizing the number of clusters

Several plotting functions (plot_cell_clusters, plot_cell_trajectory, plot_genes_in_pseudotime) need size factor normalized expression matrices. They should use a common subroutine for doing this. Right now they all have different versions of the same basic code, with subtle differences.

@hpliner This no longer works in monocle3 alpha and may need to be changed even more w/r/t slot updating, etc.