SARS-CoV-2感染的雪貂支气管肺泡灌洗液单细胞转录组数据挖掘（4）巨噬细胞内部细分亚群的拟时序和RNA速率分析

## 安装所需要的R包
chooseCRANmirror(graphics=FALSE)
install.packages('Seurat')
install.packages("tidyverse")
install.packages("ggpubr")
install.packages("cowplot")
install.packages("dplyr")
install.packages("ggplot2")
install.packages("tidyr")
install.packages("msigdbr")

if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("piano")
BiocManager::install("HSMMSingleCell")
BiocManager::install("monocle")


## 载入R包
library(Seurat)
library(tidyverse)
library(ggpubr)
library(cowplot)
library(dplyr)
library(ggplot2)
library(tidyr)
library(piano)
library(msigdbr)
library(monocle)


##载入巨噬细胞
MP <- readRDS("Seurat_object_monocyte_macrophage_DC.Rds")
## define cell cluster name 
MP$Annotation <- as.character(MP$Annotation)
MP$Annotation_2 <- as.character(MP$Annotation)

MP$Annotation[MP$Annotation_2 == "APOE pos FABP4 high tissue M2"] <- "APOE+ tissue macrophage"
MP$Annotation[MP$Annotation_2 == "SPP1 high fibrogenic M2"] <- "SPP1hi CHIT1int profibrogenic M2"
MP$Annotation[MP$Annotation_2 == "Transitional M1"] <- "Weakly activated M1"
MP$Annotation[MP$Annotation_2 == "Interferon stimulated M1"] <- "Highly activated M1"
MP$Annotation[MP$Annotation_2 == "Infiltrating macrophage"] <- "Monocyte-derived infiltrating macrophage"
MP$Annotation[MP$Annotation_2 == "APOE pos FABP4 high tissue M2"] <- "APOE+ tissue macrophage"

## Fig5 b

#伪时序分析计算

SO.traj1 <- MP[,MP$Annotation %in% c("Monocyte-derived infiltrating macrophage", "Weakly activated M1", "Highly activated M1")]
pd <- new('AnnotatedDataFrame', data = SO.traj1@meta.data)
fd <- new('AnnotatedDataFrame', data = data.frame(gene_short_name = row.names(SO.traj1), row.names = row.names(SO.traj1)))
cds <- newCellDataSet(as(SO.traj1@assays$SCT@data, "matrix"), phenoData = pd, featureData = fd, expressionFamily = negbinomial.size())
cds$clusters <- SO.traj1$Annotation

#估计规模因素和离散度
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)

#过滤低质量的细胞
set.seed(123)
cds <- detectGenes(cds, min_expr = 0.1)
markers <- FindAllMarkers(object = SO.traj1, min.pct = 0.25, thresh.use = 0.25)
markers <- subset(markers, p_val_adj < 0.05)
order.genes <- unique(as.character(markers$gene))

#构建单细胞轨迹
cds <- setOrderingFilter(cds, order.genes)
cds <- reduceDimension(cds = cds, max_components = 3,method = 'DDRTree')
cds <- orderCells(cds)
pdf("Fig5d.MP_M1_pseudotime_trajectory.pdf", 400/72,200/72)
plot_cell_trajectory(cds, color_by = "clusters", show_branch_points = F, cell_size = 0.5) + scale_x_reverse()
dev.off()

#构建单细胞UMAP轨迹
head(MP@reductions$umap@cell.embeddings)
head(cds@phenoData@data)
MP@meta.data<-MP@meta.data[rownames(cds@phenoData@data),]
MP@meta.data$Pseudotime <- cds@phenoData@data$Pseudotime
head(MP@meta.data)
mydata<- FetchData(MP,vars = c("UMAP_1","UMAP_2","Pseudotime"))

p <- ggplot(mydata,aes(x = UMAP_1,y =UMAP_2,colour = Pseudotime))+
geom_point(size = 1)+
scale_color_gradientn(values = seq(0,1,0.2),
colours = c('blue','green','yellow'))

p1 <- p + theme_bw() + theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"))
p1

## Fig 5c伪时序分析热图
order.genes <- order.genes[!grepl("ENSMPUG", order.genes)]
pdf("Fig5e.MP_M1_pseudotime_heatmap.pdf", 5/2.54*2,9/2.54*1.5)
plot_pseudotime_heatmap(cds[order.genes,], num_clusters = 4, cores = 1, show_rownames = F, return_heatmap = T)
dev.off()

## Fig5 d,e

#伪时序分析计算

SO.traj2 <- MP[,MP$Annotation %in% c("Monocyte-derived infiltrating macrophage", "SPP1hi CHIT1int profibrogenic M2")]
pd <- new('AnnotatedDataFrame', data = SO.traj2@meta.data)
fd <- new('AnnotatedDataFrame', data = data.frame(gene_short_name = row.names(SO.traj2), row.names = row.names(SO.traj2)))
cds <- newCellDataSet(as(SO.traj2@assays$SCT@data, "matrix"), phenoData = pd, featureData = fd, expressionFamily = negbinomial.size())

cds$clusters <- SO.traj2$Annotation


#估计规模因素和离散度
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)

#过滤低质量的细胞
set.seed(123)
cds <- detectGenes(cds, min_expr = 0.1)
markers <- FindAllMarkers(object = SO.traj2, min.pct = 0.25, thresh.use = 0.25)
markers <- subset(markers, p_val_adj < 0.05)
order.genes <- unique(as.character(markers$gene))

##构建单细胞轨迹
cds <- setOrderingFilter(cds, order.genes)
cds <- reduceDimension(cds = cds, max_components = 3,method = 'DDRTree')
cds <- orderCells(cds)
pdf("Fig5d.MP_M2_pseudotime_trajectory.pdf", 400/72,200/72)
plot_cell_trajectory(cds, color_by = "clusters", show_branch_points = F, cell_size = 0.5) + scale_x_reverse()
dev.off()

#构建单细胞UMAP轨迹
head(MP@reductions$umap@cell.embeddings)
head(cds@phenoData@data)
MP@meta.data<-MP@meta.data[rownames(cds@phenoData@data),]
MP@meta.data$Pseudotime <- cds@phenoData@data$Pseudotime
head(MP@meta.data)
mydata<- FetchData(MP,vars = c("UMAP_1","UMAP_2","Pseudotime"))
p <- ggplot(mydata,aes(x = UMAP_1,y =UMAP_2,colour = Pseudotime))+
geom_point(size = 1)+
scale_color_gradientn(values = seq(0,1,0.2),
colours = c('blue','green','yellow'))
p1 <- p + theme_bw() + theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"))
p1

## Fig 5e伪时序分析热图
order.genes <- order.genes[!grepl("ENSMPUG", order.genes)]
pdf("Fig5e.MP_M2_pseudotime_heatmap.pdf", 5/2.54*2,9/2.54*1.5)
plot_pseudotime_heatmap(cds[order.genes,], num_clusters = 4, cores = 1, show_rownames = F, return_heatmap = T)
dev.off()