library(stats4)
library(splines)
library(VGAM)
library(parallel)
library(irlba)
library(Matrix)
library(DDRTree)
library(BiocGenerics)
library(Biobase)
library(ggplot2)
library(Seurat)
library(monocle)

#Load object
#subset myogenic clusters only
subset = subset(muscle, idents = c("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "12", "15"))
DimPlot(object = subset, reduction = "umap")

#downsampling to 500 cells per cluster
subset.small <- subset(subset, downsample = 500)

#Rename seurat object
seuratX <- subset.small

#Extract data, phenotype data, and feature data from the SeuratObject
data <- as(as.matrix(seuratX@assays$RNA@data), 'sparseMatrix')

pd <- new('AnnotatedDataFrame', data = seuratX@meta.data)

fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data))
fd <- new('AnnotatedDataFrame', data = fData)

#Construct monocle cds
HSMM <- newCellDataSet(data,
                       phenoData = pd,
                       featureData = fd,
                       #lowerDetectionLimit = 0.5,
                       expressionFamily = uninormal())# since I have already normalized, thresholded and scalled in Suerat v3.0.0.9150

#View data
pData(HSMM)
fData(HSMM)


#Run ordering algorithm
var_genes <-seuratX[["RNA"]]@var.features
ordering_genes <- var_genes

HSMM <- setOrderingFilter(HSMM, ordering_genes)
print(dim(exprs(HSMM)))

## reduce dimension - do not normalize or include pseudo count. Use monocle scaling
HSMM <- reduceDimension(HSMM,norm_method="none", 
                        reduction_method="DDRTree",
                        max_components=4,
                        scaling=TRUE,
                        verbose=TRUE,
                        pseudo_expr=0)

# First decide what you want to color your cells by
print(head(pData(HSMM)))

## order cells change colors and theta to match your plot
HSMM <- orderCells(HSMM)
plot_cell_trajectory(HSMM, 
                     color_by = "seurat_clusters",
                     theta = -15,
                     show_branch_points = FALSE,
                     show_tree = TRUE,
                     cell_size = 1))


marker.genes <- c("PAX7", "MYF5", "MYOD1", "MYOG")

marker.genes2 <- c("SPRY1", "HEY1", "DLK1", "CAV1")

plot_genes_branched_pseudotime(HSMM[marker.genes,],
                               branch_point = 3,
                               color_by = "seurat_clusters",
                               cell_size = 0.5,
                               ncol = 2)








#pseudotemporalplot