./steel filtered_feature_bc_matrix  spatial/tissue_positions_list.csv --pca=20  pca20out

library (Seurat)

library (dplyr)

library (ggplot2)

library (magrittr)

library (cowplot)

library (gtools)

library (stringr)

library (Matrix)

library (tidyverse)

library (patchwork)

orc1<- Load10X_Spatial ("Slide_1")

## 直接导入 STEEL 的聚类信息

STEEL<- read.delim ("STEEL/Slide1.map.40",row.names = 1)

orc1 [["seurat_clusters"]] <- NA

clusters<-data.frame (STEEL$Cluster)

rownames (clusters) <- rownames (STEEL)

orc1 [["seurat_clusters"]][rownames (clusters),] <- clusters

## 去除 NA 值只保留有聚类的

orc1_new <- orc1 [,rownames (clusters)]

Idents (orc1_new) <- 'seurat_clusters'

Idents (orc1_new) <- factor (Idents (orc1_new),levels=mixedsort (levels (Idents (orc1_new))))

orc1_new <- SCTransform(orc1_new, assay = "Spatial", return.only.var.genes = FALSE, verbose = FALSE)

orc1_new <- RunPCA(orc1_new, features = VariableFeatures(orc1_new)) 

orc1_new <- RunTSNE(orc1_new, dims = 1:20)

p1 <- DimPlot(orc1_new, reduction = "tsne", label = TRUE)

p2 <- SpatialDimPlot(orc1_new, group.by = "seurat_clusters",label.size = 3, pt.size.factor = 1.3)

pearplot <- plot_grid(p1,p2)

ggsave("STEEL/tsne_Slide1_40.pdf", plot = pearplot, width = 6, height = 6)

p1<- SpatialDimPlot(orc1_new, cells.highlight = CellsByIdentities(object = orc1_new, idents = c(1:40)), facet.highlight = TRUE, ncol = 5)

ggsave("STEEL/cluster_Slide1all.pdf", plot = p1, width = 19, height = 12)

library (monocle)

subdata <- subset (orc1_new, idents = c (19,21,37,38,39,40))

# 选择要分析的亚群

expression_matrix = subdata@assays$Spatial@counts

cell_metadata <- data.frame (group = subdata [['orig.ident']],clusters = Idents (subdata))

gene_annotation <- data.frame (gene_short_name = rownames (expression_matrix), stringsAsFactors = F) 

rownames (gene_annotation) <- rownames (expression_matrix)

pd <- new ("AnnotatedDataFrame", data = cell_metadata)

fd <- new ("AnnotatedDataFrame", data = gene_annotation)

HSMM <- newCellDataSet (expression_matrix,

                       phenoData = pd,

                       featureData = fd,

                       expressionFamily=negbinomial.size ())

HSMM <- detectGenes (HSMM,min_expr = 0.1)

HSMM_myo <- estimateSizeFactors (HSMM)

HSMM_myo <- estimateDispersions (HSMM_myo)

disp_table <- dispersionTable (HSMM_myo)

disp.genes <- subset (disp_table, mean_expression >= 0.1 )

disp.genes <- as.character (disp.genes$gene_id)

HSMM_myo <- reduceDimension (HSMM_myo, max_components = 2, method = 'DDRTree')

HSMM_myo <-orderCells (HSMM_myo,reverse = T)

#State 轨迹分布图

plot1 <- plot_cell_trajectory (HSMM_myo, color_by = "State")

##Cluster 轨迹分布图

plot2 <- plot_cell_trajectory (HSMM_myo, color_by = "clusters")

##Pseudotime 轨迹图

plot3 <- plot_cell_trajectory (HSMM_myo, color_by = "Pseudotime")

plotc <- plot1|plot2|plot3

ggsave ("STEEL/Combination1.pdf", plot = plotc, width = 18, height = 6.2)

# 绘制拟时间

cell_Pseudotime <- data.frame (pData (HSMM_myo)$Pseudotime)

rownames (cell_Pseudotime) <- rownames (cell_metadata) 

# 把拟时间对应到到组织切片位置上

orc1_new [['Pseudotime']] <- NA

orc1_new [['Pseudotime']][rownames (cell_Pseudotime),] <- cell_Pseudotime

p1 <- SpatialFeaturePlot (orc1_new, features = c ("Pseudotime"),pt.size.factor = 1.3)

ggsave ("STEEL/pseudotime_feature1.pdf", plot = p1, width = 8, height = 9)

data_subset <- HSMM_myo ['PAXXG054350',]

p1<-plot_genes_in_pseudotime (data_subset, color_by = "clusters")

data_subset <- HSMM_myo ['PAXXG051950',]

p2<-plot_genes_in_pseudotime (data_subset, color_by = "clusters")

data_subset <- HSMM_myo ['PAXXG086750',]

p3<-plot_genes_in_pseudotime (data_subset, color_by = "clusters")

data_subset <- HSMM_myo ['PAXXG345890',]

p4<-plot_genes_in_pseudotime (data_subset, color_by = "clusters")

data_subset <- HSMM_myo ['PAXXG010560',]

p5<-plot_genes_in_pseudotime (data_subset, color_by = "clusters")

data_subset <- HSMM_myo ['PAXXG074500',]

p6<-plot_genes_in_pseudotime (data_subset, color_by = "clusters") #color_by 可以换成 state 或者 pseudotime

pearplot <- plot_grid (p1,p2,p3,p4,p5,p6,align = "v",axis = "b",ncol = 1)

ggsave ("STEEL/gene_pseudotime1.pdf", plot = pearplot, width = 5, height = 15)

# 拟时相关基因聚类热图   

disp_table <- dispersionTable (HSMM_myo)

disp.genes <- subset (disp_table, mean_expression >= 0.5&dispersion_empirical >= 1*dispersion_fit)

disp.genes <- as.character (disp.genes$gene_id)

mycds_sub <- HSMM_myo [disp.genes,]

diff_test <- differentialGeneTest (HSMM_myo [disp.genes,], cores = 4, 

                                  fullModelFormulaStr = "~sm.ns (Pseudotime)")

sig_gene_names <- row.names (subset (diff_test, qval < 1e-04))

p2 = plot_pseudotime_heatmap (HSMM_myo [sig_gene_names,], num_clusters=6,

                             show_rownames=F, return_heatmap=T)

ggsave ("STEEL/pseudotime_heatmap1.pdf", plot = p2, width = 5, height = 10)

## BEAM 分析

my_pseudotime_de <- differentialGeneTest (HSMM_myo, cores = 5)

BEAM_res <- BEAM (HSMM_myo, branch_point = 1, cores = 4)

BEAM_res <- BEAM_res [order (BEAM_res$qval),]

BEAM_res <- BEAM_res [,c ("gene_short_name", "pval", "qval")]

mycds_sub_beam <- HSMM_myo [row.names (subset (BEAM_res, qval < 1e-4)),]

pdf ("STEEL/BEAM1.pdf", width = 8, height = 12)

plot_genes_branched_heatmap (mycds_sub_beam, branch_point = 1, num_clusters = 6, cores = 4,show_rownames = FALSE)

dev.off ()