statOmics / tradeSeq

TRAjectory-based Differential Expression analysis for SEQuencing data
Other
237 stars 27 forks source link

Unable to determine root cell #38

Closed Seandelao closed 4 years ago

Seandelao commented 4 years ago

Hello,

Thank you for the wonderful package, I am very excited to utilize it for my data sets!

I have created a cds object from an integrated Seurat object, but an having issues with calculating pseudotime by following the code provided by the tradeSeq tutorial. I can calculate it fine following Monocle3's method in their tutorial, but then I don't know how to modify the tradeSeq tutorial to extract the pseudotime and cellweight values.

When I run: exprs_human<- GetAssayData(Alpha_Beta, slot="counts", assay = "RNA") counts <- as.matrix(exprs_human)

meta_data <- Alpha_Beta@meta.data keep <- ("CellfindR") meta_data <- meta_data[keep] names(meta_data)[names(meta_data)=="CellfindR"] <- "cellType"

df <- data.frame(cells = colnames(counts), cellType = meta_data) rownames(df) <- df$cells

cds <- new_cell_data_set(counts, cell_metadata = df, gene_metadata = data.frame(gene_short_name = rownames(counts), row.names = rownames(counts)))

cds <- preprocess_cds(cds, method = "PCA") cds <- reduce_dimension(cds, preprocess_method = "PCA", reduction_method = "UMAP")

plot_cells(cds, label_groups_by_cluster = FALSE, cell_size = 1, color_cells_by = "cellType")

umaps = Embeddings(Alpha_Beta, assay="integrated",reduction="umap") reducedDims(cds)$UMAP <- umaps plot_cells(cds,color_cells_by="cellType")

cds <- cluster_cells(cds, reduction_method = "UMAP")

cds <- learn_graph(cds) plot_cells(cds, label_groups_by_cluster = FALSE, cell_size = 1, color_cells_by = "cellType")

celltype <- as.character(meta_data$cellType) cell_ids <- which(celltype == "3.1.0") closest_vertex <- cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex closest_vertex <- closest_vertex[colnames(cds), 1] closest_vertex <- closest_vertex[cell_ids] closestvertex <- paste0("Y", closest_vertex) root <- names(which(igraph::degree(principal_graph(cds)[["UMAP"]]) == 1)) root <- root[root %in% closest_vertex] cds <- order_cells(cds, root_pr_nodes = root)

I get the error: Error in data_matrix[, ((((i - 1) * block_size) + 1):(ncol(data_matrix)))] : subscript out of bounds

Everything seems to be working well until I run the line: root <- root[root %in% closest_vertex]

in which it seems to return an empty character string.

Any idea what might be happening? Or, alternatively, how I might be able to extract the pseudotime and cell weights from the cds object after selecting the root cell by the current method suggested by Monocle3's tutorial? I don't have the data 'root' if I use their method...

Sorry if that was confusing, happy to clarify! Thanks again.

-Sean

HectorRDB commented 4 years ago

Hi, Thanks for running tradeSeq, I'm glad you like it.

For now, the monocle3 / tradeSeq part is very experimental since monocle3 itself is still in its beta version. So the code in the monocle / tradeSeq vignette is resorting to two hacks:

So now, to your issue. I would guess that the problem is that there is no root point (i.e a terminal node in the graph) in that cluster. That should appear when running those lines

closest_vertex <- paste0("Y_", closest_vertex)
root <- names(which(igraph::degree(principal_graph(cds)[["UMAP"]]) == 1))
root <- root[root %in% closest_vertex]

There is not intersect between the roots and the closest_vertex.

This code is inspired from https://cole-trapnell-lab.github.io/monocle3/docs/trajectories/, the get_earliest_principal_node function.

I would recommend either

Let me know if this is clear.

Seandelao commented 4 years ago

Hi @HectorRDB

That makes perfect sense, thank you! I will try both options. -Sean

HectorRDB commented 4 years ago

Great let me know how that goes. Also @Seandelao if you end up opening an issue on monocle3, could you please refer to this one as well so future tradeSeq users can know where to look. Thanks