jackbibby1 / SCPA

R package for pathway analysis in scRNA-seq data
https://jackbibby1.github.io/SCPA/
GNU General Public License v3.0
57 stars 6 forks source link

Odd result in pseudotime #73

Closed Umaarasu closed 1 month ago

Umaarasu commented 3 months ago

@jackbibby1 ...I have been trying to dun your pseudotime codes in my own data. The data entered as per the instruction, but when I run the plot command this is the figure that comes

plot_dimred(model_n4, "pseudotime", pseudotime = calculate_pseudotime(model_n4), hex_cells = F, plot_trajectory = T, size_cells = 1, alpha_cells = 0.8) + theme(aspect.ratio = 1) root cell or milestone not provided, trying first outgoing milestone_id Using '3' as root Warning message: In calculate_pseudotime(model_n4) : Trajectory is not rooted. Add a root to the trajectory using dynwrap::add_root(). This will result in an error in future releases.

Rplot02

When I added a root cell root_cell <- "celltype1" model_n4 <- add_root(model_n4, root_cell) Error in add_root(model_n4, root_cell) : Invalid root_cell_id

I notice that the "dataset_n4" has cell id as barcodes AAACCAACACAGCCAT-1.

Could you help me fix this? Thanks :)

jackbibby1 commented 3 months ago

Hi,

Is the code for the trajectory calculation exactly as per here?

Jack

jackbibby1 commented 3 months ago

I've just re-ran this code on a fresh installation of dyno and everything looked fine

library(Seurat)
library(tidyverse)
library(dyno)

naive_cd4 <- readRDS("~/Downloads/naive_cd4.rds")
naive_cd4 <- UpdateSeuratObject(naive_cd4)

naive_cd4 <- subset(naive_cd4, idents = "Treg", invert = T)

df <- as.matrix(naive_cd4[["RNA"]]@data)
var_genes <- names(sort(apply(df, 1, var), decreasing = TRUE))[1:1000]

counts <- Matrix::t(as(as.matrix(naive_cd4@assays$RNA@counts[var_genes,]), 'sparseMatrix'))
expression <- Matrix::t(as(as.matrix(naive_cd4@assays$RNA@data[var_genes,]), 'sparseMatrix'))

dataset_n4 <- wrap_expression(expression = expression,
                              counts = counts)

model_n4 <- infer_trajectory(dataset_n4, method = ti_slingshot(), verbose = T)

plot_dimred(model_n4, 
            "pseudotime", 
            pseudotime = calculate_pseudotime(model_n4), 
            hex_cells = F,
            plot_trajectory = T, 
            size_cells = 1, alpha_cells = 0.8) + 
  theme(aspect.ratio = 1)

--------

root cell or milestone not provided, trying first outgoing milestone_id
Using '1' as root
Warning message:
In calculate_pseudotime(model_n4) :
  Trajectory is not rooted. Add a root to the trajectory using dynwrap::add_root(). This will result in an error in future releases.

image

Umaarasu commented 2 months ago

Hi,

Is the code for the trajectory calculation exactly as per here?

Jack

@jackbibby1 I ran the code same like the example. One difference between the example data set and mine is that I cannot split the data by anything. It is from a bunch of diseased samples. I know that there is a trajectory as I have checked it with monocle.

In your example code when you run the plot command it says "root cell or milestone not provided, trying first outgoing milestone_id Using '1' as root". For me it says using '3'. I realised maybe I should define the root cell for a clean trajectory. But the problem is the object containing trajectory info "model_n4" does not have the celltype metadata. It has the barcodes stored as the cell_id's.

Could you let me know how to procees?

Thanks!

jackbibby1 commented 2 months ago

If you've ran monocle already then you can split the cells across that trajectory and use it as an input to SCPA -- it doesn't have to be slingshot.

Otherwise, this is more likely a question for the Slingshot authors and asking what they would advise

Umaarasu commented 2 months ago

If you've ran monocle already then you can split the cells across that trajectory and use it as an input to SCPA -- it doesn't have to be slingshot.

Otherwise, this is more likely a question for the Slingshot authors and asking what they would advise @jackbibby1 Would it be possible for you to give a example code of how you can feed a monocle output into SCPA trajectory. Thanks :)

jackbibby1 commented 2 months ago

This is likely not the best or most elegant way to do it, but it should work. You can split the trajectory into x number of parts across all values, then use these to annotate your Seurat object or whatever you're using e.g.

pseudotime <- seq(0, 7.5, 0.1) ## this is just a distribution of values representing your actual trajectory values
trajectory_groups <- case_when(pseudotime <= 2.5 ~ "one",
                               pseudotime > 2.5 & pseudotime <= 5 ~ "two",
                               pseudotime > 5 & pseudotime <= 7.5 ~ "three")
seurat_obect$trajectory_groups <- trajectory_groups

... then use the trajectory groups in SCPA