Monocle3 script

BiocManager::install(c('BiocGenerics', 'DelayedArray', 'DelayedMatrixStats',
                       'limma', 'lme4', 'S4Vectors', 'SingleCellExperiment',
                       'SummarizedExperiment', 'batchelor', 'HDF5Array',
                       'terra', 'ggrastr'))

devtools::install_github('cole-trapnell-lab/monocle3')

install.packages('Seurat')

devtools::install_github("satijalab/seurat-wrappers")

library(monocle3)
library(Seurat)

files <- list.files('data')
files <- grep("2i|Dox", files, value = TRUE)
files <- grep("h5", files, value = TRUE)
files

input_files <- function(h5_path){
  day <- sub("^.*_D([^_]*)_Dox.*$", "\\1", h5_path)
  day <- sub("^.*_D([^_]*)_2i.*$", "\\1", day)
  
  data <- Read10X_h5(paste0('data/', h5_path))
  data <- CreateSeuratObject(data, min.cells = 0, min.features = 200)
  data[["percent.mt"]] <- PercentageFeatureSet(data, pattern = "mt-")
  lb <- quantile(data[["nFeature_RNA"]]$nFeature_RNA, probs = 0.02)
  ub <- quantile(data[["nFeature_RNA"]]$nFeature_RNA, probs = 0.97)
  data <- data[, data[["nFeature_RNA"]] > lb & data[["nFeature_RNA"]] < ub & data[["percent.mt"]] < 15]
  
  data$day <- day
  return(data)
}

data_list <- sapply(files, input_files)
data <- merge(data_list[1]$gene_bc_mat.h5, y = data_list[2:length(data_list)])

data <- merge(data_list[1]$gene_bc_mat.h5, y = data_list[2:length(data_list)])

data <- NormalizeData(object = data, verbose = FALSE)
data <- FindVariableFeatures(object = data, nfeatures = 2000, verbose = FALSE, selection.method = 'vst')
data <- ScaleData(data, verbose = FALSE)
data <- RunPCA(data, npcs = 30, verbose = FALSE)
data <- FindNeighbors(data, dims = 1:30)

data <- RunUMAP(data, reduction = "pca", dims = 1:30)


data@active.assay = 'RNA'

DimPlot(data, group.by = c("day"))

data$dayint <- data[[]]$day
data$dayint <- ifelse(data$dayint == "iPSC", 20, data$dayint)
data$dayint <- as.numeric(data$dayint)
FeaturePlot(data, "dayint")

FeaturePlot(data, "dayint")

cds <- SeuratWrappers::as.cell_data_set(data) #change to cds here

cds <- cluster_cells(cds)

plot_cells(cds, show_trajectory_graph = FALSE,
           color_cells_by = "partition")

cds <- learn_graph(cds, use_partition = FALSE) #graph learned across all partitions

cds <- order_cells(cds)

plot_cells(cds, color_cells_by = "pseudotime", label_branch_points=FALSE, label_leaves=FALSE)

rowData(cds)$gene_name <- rownames(cds)
rowData(cds)$gene_short_name <- rowData(cds)$gene_name

plot_cells(cds,
           genes=c('Sox2', 'Nanog', 'Col6a2'),
           label_cell_groups=FALSE,
           show_trajectory_graph=FALSE, 
           min_expr = 3)

cds_pt_res <- graph_test(cds, neighbor_graph="principal_graph", cores=8)

cds_pt_res <- readRDS("cds_pt_res.rds")

cds_pt_res

cds_pt_res <- na.omit(cds_pt_res)
cds_pt_res <- cds_pt_res[cds_pt_res$p_value < 0.05 & cds_pt_res$status == "OK", ]
cds_pt_res

cds_pt_res[order(-cds_pt_res$morans_test_statistic),] 

plot_cells(cds, genes=c("Col1a2", "Uba52", "Serpine1", "Dppa5a"),
           show_trajectory_graph=FALSE,
           label_cell_groups=FALSE,
           label_leaves=FALSE)

cds_subset <- choose_cells(cds)

cds_subset

cds_subset_pt_res <- graph_test(cds_subset, neighbor_graph="principal_graph", cores=8)
cds_subset_pt_res <- na.omit(cds_subset_pt_res)
cds_subset_pt_res <- cds_subset_pt_res[cds_subset_pt_res$p_value < 0.05 & cds_subset_pt_res$status == "OK", ]
cds_subset_pt_res

cds_subset_pt_res[order(-cds_subset_pt_res$morans_test_statistic),] 

plot_cells(cds_subset, genes=c("Rpl7a", "Eef1a1", "Mgst1", "Lgals1"),
           show_trajectory_graph=FALSE,
           label_cell_groups=FALSE,
           label_leaves=FALSE)

cds_subset_subset <- cds_subset[rowData(cds_subset)$gene_short_name %in% c("Rpl7a", "Eef1a1", "Mgst1", "Lgals1")]

plot_genes_in_pseudotime(cds_subset_subset,
                         color_cells_by="dayint",
                         min_expr=0.5)

Integration Example

features <- SelectIntegrationFeatures(object.list = data_list)

scale_pca <- function(x){
  x <- ScaleData(x, features = features, verbose = FALSE)
  x <- RunPCA(x, features = features, verbose = FALSE)
  return(x)
}


data_list <- lapply(X = data_list, scale_pca)

anchors <- FindIntegrationAnchors(object.list = data_list, anchor.features = features, reduction = "rpca")
saveRDS(anchors, file = "integration_anchors.rds")


data <- IntegrateData(anchorset = anchors)