Open cohenp05 opened 5 years ago
Hello @cohenp05,
Could you provide a minimum reproducible example for which your problem occurs? For now, I haven't managed to replicate it using the following code.
Robrecht
library(dyno)
library(tidyverse)
dataset <- dyntoy::generate_dataset()
traj <- infer_trajectory(dataset, ti_slingshot())
metric_ids <- c("correlation", "him", "featureimp_wcor", "F1_branches")
dyneval::calculate_metrics(dataset = dataset, model = traj, metrics = metric_ids, expression_source = dataset$expression)
# A tibble: 1 x 12
time_waypointedgeodesic correlation time_correlation him time_him time_featureimp featureimp_cor featureimp_wcor time_mapping_branches recovery_branches relevance_branches F1_branches
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 0.0878 0.949 0.00342 1 0.0263 1.68 0.797 0.687 0.104 1 1 1
It's a bit weird that calculate_metrics doesn't error on https://github.com/dynverse/dyneval/blob/master/R/calculate_metrics.R#L49 because the dataset in your case does not contain cell waypoints (yet).
Apart from what @rcannood suggested, could you perhaps also try to add waypoints by:
monocyte_dyn <- dynwrap::add_cell_waypoints(monocyte_dyn)
and see whether it errors in that case?
Hi @cohenp05
Apart from what @rcannood said, could you show us the names()
of the dataset and model object? It looks like that neither of these objects contain a trajectory (i.e. milestone network, milestone percentages, ...)
Hello,
Thanks for your help. I've made an example below, but the code is still very lengthy and certainly overly verbose. Sorry about that. Anyway here is what I've written, and I still get the same error.
In order to run this I'm using a dataset that I had previously analyzed in the R package Seurat and also in Monocle. Those objects are called dataset_seurat and dataset_monocle.
library(dyno)
library(dyneval)
library(dyntoy)
#Check if docker is correctly installed
dynwrap::test_docker_installation(detailed = TRUE)
library(tidyverse)
library(Seurat)
library(Matrix)
library(monocle)
library(uwot)
library(rlist)
#Pull counts data from Seurat, log2 normalize, and load into Dynwrap
dataset_counts <- GetAssayData(object = dataset_seurat, slot = "counts", assay = "RNA")
dataset_norm <- log2(dataset_counts+1)
dataset_norm <- as(dataset_norm, "dgCMatrix")
dataset_norm <- t(dataset_norm)
dataset_counts <- t(dataset_counts)
dataset_dyn <- wrap_expression(expression = dataset_norm, counts = dataset_counts)
#Pull barcodes and pseudotime assignments from Monocle object of corresponding data
barcode <- rownames(pData(dataset_monocle))
pseudotime <- pData(dataset_monocle)$Pseudotime
cellsinpseudotime <- data.frame(cells = barcode, pseudotime = pseudotime)
#Use only barcodes that are shared between dynwrap object and monocle object
cellsindyn <- dataset_dyn$cell_ids
cellsinpseudotime <- cellsinpseudotime[barcode %in% cellsindyn,]
#Define start cell and end cell as minimum pseudotime value
# and maximum pseudotime value in monocle analysis
startcell <- as.character(cellsinpseudotime[pseudotime==min(pseudotime),1])
endcell <- as.character(cellsinpseudotime[pseudotime==max(pseudotime),1])
#Add start and end_id as prior information
dataset_dyn <- add_prior_information(dataset = dataset_dyn, start_id = startcell, end_id = endcell, start_n = 1, end_n = 1)
method_tis <- list(ti_comp1(), ti_mst(), ti_shuffle())
set.seed(seed = 456)
dataset_model <- list()
##Make 3 TI models
##Repeat these methods 3 times and store all 9 models in the list object dataset_model
for (i in seq(3)){
model_method <- list()
for (j in seq(length(method_tis))){
model_method[[j]] <- infer_trajectories(dataset = dataset_dyn, method = method_tis[j], verbose = TRUE, give_priors = c("start_id", "end_id"))
}
dataset_model[[i]] <- model_method
}
#Root All Models by CD14 Expression
for (i in seq(length(dataset_model))){
for (j in seq(length(method_tis))){
if ("dynwrap::with_trajectory" %in% class(dataset_model[[i]][[j]][[6]][[1]])){
dataset_model[[i]][[j]][[6]][[1]] <- add_root_using_expression(dataset_model[[i]][[j]][[6]][[1]], features_oi = "CD14", expression_source = dataset_dyn$expression)
}
}
}
#Calculate Pseudotime for All Models
for (i in seq(dataset_model)){
for (j in seq(length(method_tis))){
if ("dynwrap::with_trajectory" %in% class(dataset_model[[i]][[j]][[6]][[1]])){
dataset_model[[i]][[j]][[6]][[1]]<- add_pseudotime(dataset_model[[i]][[j]][[6]][[1]])
}
}
}
#Add cell waypoints for all models
for (i in seq(length(dataset_model)){
for (j in seq(length(method_tis))){
if (is_wrapper_with_trajectory(dataset_model[[i]][[j]][[6]][[1]]) == TRUE){
dataset_model[[i]][[j]][[6]][[1]] <- add_cell_waypoints(trajectory = dataset_modelcp[[i]][[j]][[6]][[1]])
}
}
}
##Select metrics to calculate: include all metrics except metrics that have the category "average"
metric_ids <- dyneval::metrics %>% filter(category != "average") %>% pull(metric_id)
metrics <- list()
##Calculate Metrics
for (i in seq(length(dataset_model))){
metric_trial <- list()
for (j in seq(length(method_tis))){
if (is_wrapper_with_trajectory(dataset_model[[i]][[j]][[6]][[1]]) == TRUE){
metric_trial[j] <- map_dfr(dataset_model, dyneval::calculate_metrics, dataset = dataset_dyn, metrics = metric_ids)
}
}
metrics[i] <- metric_trial
}
The names of the dataset and (one example) model object are below:
> names(dataset_dyn)
[1] "id" "cell_ids" "cell_info" "counts" "expression" "expression_projected"
[7] "feature_info" "feature_ids" "prior_information"
> names(dataset_model[[1]][[1]][[6]][[1]])
[1] "id" "cell_ids" "cell_info" "milestone_ids"
[5] "milestone_network" "divergence_regions" "milestone_percentages" "progressions"
[9] "trajectory_type" "directed" "pseudotime" "dimred"
[13] "dimred_milestones" "dimred_segment_progressions" "dimred_segment_points" "root_milestone_id"
Hey @cohenp05,
In order to calculate any metrics, your dataset_dyn
needs to have a gold standard trajectory. For this, you need to run something like:
dataset_dyn <-
wrap_expression(
expression = dataset_norm,
counts = dataset_counts
) %>%
add_trajectory(
milestone_network = ...,
milestone_percentages = ...
)
Where you define milestone network and milestone percentage according to the trajectory that is in the dataset. Check ?add_trajectory
for more information of the format of the milestone network and percentages. Do you have this information?
Robrecht
Hey @rcannood,
Thanks so much for your help, and sorry for such a long lapse in response. Looking back I think perhaps my question is a bit naive and I think that I wasn't communicating it well. Ideally, I'd like to compare several trajectories generated using dynverse in an experimental setting where the ground truth is unknown. I was wondering if there are commands in dyneval (or elsewhere) that can be used to determine a "best fit" for the data. Your group has outlined a wealth of tests to run in the supplement, but in an experimental setting without a known ground truth, is there a better way to choose a trajectory (once several have been built)? Thank you again for your help
Phil C
Hello,
Thanks for a great package and publication! I am trying to use dyneval to compare some trajectories that I've built using infer_trajectories, and I looped over it to run the same set of methods n=10 times on the same dataset. However, when I run the dyneval::calculate_metrics call I get the error
Error in UseMethod("rename") : no applicable method for 'rename' applied to an object of class "NULL"
Here is the call I ran. The models are stored in a list object titled Trial 1 containing a list object that houses the output of infer_trajectories (per model run).
dyneval::calculate_metrics(dataset = monocyte_dyn, model = trial1$Method_1$model[[1]], metrics = metric_ids, expression_source = monocyte_dyn$expression)
I checked the class of each input. Which is shown below:
Here is my sessionInfo():
Thank you for your help!