# ## 安装所需要的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)
#设置镜像
#options(repos="http://mirrors.tuna.tsinghua.edu.cn/CRAN/")
#options(BioC_mirror="http://mirrors.tuna.tsinghua.edu.cn/bioconductor/")
##载入巨噬细胞
setwd("~/Rdata/sars_scRNA/")
MP <- readRDS("GSE171828_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)))
data <- as(as.matrix(SO.traj1@assays$SCT@data), 'sparseMatrix')
cds <- newCellDataSet(data, phenoData = pd, featureData = fd, expressionFamily = negbinomial.size())
#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)
> head(MP@reductions$umap@cell.embeddings)
UMAP_1 UMAP_2
Ctrl-2:AAACCCACAATAGTAGx -4.5643387 4.824413
Ctrl-2:AAACCCACACCAGCTGx -5.9036785 3.614603
Ctrl-2:AAACCCACATGACTGTx 3.2230615 -11.148943
Ctrl-2:AAACCCAGTAGAGATTx 0.8020717 1.683502
Ctrl-2:AAACCCAGTCATAAAGx 3.4056079 4.289331
Ctrl-2:AAACCCATCAAGGACGx -5.1570769 4.094346
> head(cds@phenoData@data)
orig.ident nCount_RNA nFeature_RNA Sample_name Experimental_group nCount_SCT
Ctrl-2:AAACCCAGTCATAAAGx SeuratProject 7550.153 2237 Ctrl-2 Ct 8394
Ctrl-2:AAAGTGAAGAAGTGTTx SeuratProject 7832.050 2136 Ctrl-2 Ct 8429
Ctrl-2:AACAAGAAGTTGAAGTx SeuratProject 2246.858 991 Ctrl-2 Ct 7883
Ctrl-2:AACAGGGCATATCTCTx SeuratProject 12854.146 3141 Ctrl-2 Ct 9512
Ctrl-2:AACGAAACAAGACCTTx SeuratProject 4119.166 1723 Ctrl-2 Ct 7721
Ctrl-2:AAGAACAGTGAGTGACx SeuratProject 6715.242 1945 Ctrl-2 Ct 8183
nFeature_SCT SCT_snn_res.0.3 seurat_clusters SCT_snn_res.0.5 SCT_snn_res.0.4
Ctrl-2:AAACCCAGTCATAAAGx 2225 2 9 8 4
Ctrl-2:AAAGTGAAGAAGTGTTx 2135 2 9 8 4
Ctrl-2:AACAAGAAGTTGAAGTx 1552 5 8 6 4
Ctrl-2:AACAGGGCATATCTCTx 3106 2 9 8 0
Ctrl-2:AACGAAACAAGACCTTx 1813 2 8 6 4
Ctrl-2:AAGAACAGTGAGTGACx 1945 5 8 6 4
Annotation SCT_snn_res.1 SCT_snn_res.0.6
Ctrl-2:AAACCCAGTCATAAAGx Monocyte-derived infiltrating macrophage 10 9
Ctrl-2:AAAGTGAAGAAGTGTTx Monocyte-derived infiltrating macrophage 5 9
Ctrl-2:AACAAGAAGTTGAAGTx Weakly activated M1 9 8
Ctrl-2:AACAGGGCATATCTCTx Monocyte-derived infiltrating macrophage 16 9
Ctrl-2:AACGAAACAAGACCTTx Weakly activated M1 9 8
Ctrl-2:AAGAACAGTGAGTGACx Weakly activated M1 9 8
Annotation_2 Size_Factor clusters
Ctrl-2:AAACCCAGTCATAAAGx Infiltrating macrophage 1.0953478 Monocyte-derived infiltrating macrophage
Ctrl-2:AAAGTGAAGAAGTGTTx Infiltrating macrophage 1.0700055 Monocyte-derived infiltrating macrophage
Ctrl-2:AACAAGAAGTTGAAGTx Transitional M1 0.8270419 Weakly activated M1
Ctrl-2:AACAGGGCATATCTCTx Infiltrating macrophage 1.4714588 Monocyte-derived infiltrating macrophage
Ctrl-2:AACGAAACAAGACCTTx Transitional M1 0.9654185 Weakly activated M1
Ctrl-2:AAGAACAGTGAGTGACx Transitional M1 1.0523062 Weakly activated M1
num_genes_expressed Pseudotime State
Ctrl-2:AAACCCAGTCATAAAGx 2225 28.72614 5
Ctrl-2:AAAGTGAAGAAGTGTTx 2135 27.32175 5
Ctrl-2:AACAAGAAGTTGAAGTx 1552 16.39644 5
Ctrl-2:AACAGGGCATATCTCTx 3106 31.99751 5
Ctrl-2:AACGAAACAAGACCTTx 1813 26.84087 5
Ctrl-2:AAGAACAGTGAGTGACx 1945 27.36084 5
MP@meta.data<-MP@meta.data[rownames(cds@phenoData@data),]
MP@meta.data$Pseudotime <- cds@phenoData@data$Pseudotime
head(MP@meta.data)
head(MP@meta.data)
orig.ident nCount_RNA nFeature_RNA Sample_name Experimental_group nCount_SCT
Ctrl-2:AAACCCAGTCATAAAGx SeuratProject 7550.153 2237 Ctrl-2 Ct 8394
Ctrl-2:AAAGTGAAGAAGTGTTx SeuratProject 7832.050 2136 Ctrl-2 Ct 8429
Ctrl-2:AACAAGAAGTTGAAGTx SeuratProject 2246.858 991 Ctrl-2 Ct 7883
Ctrl-2:AACAGGGCATATCTCTx SeuratProject 12854.146 3141 Ctrl-2 Ct 9512
Ctrl-2:AACGAAACAAGACCTTx SeuratProject 4119.166 1723 Ctrl-2 Ct 7721
Ctrl-2:AAGAACAGTGAGTGACx SeuratProject 6715.242 1945 Ctrl-2 Ct 8183
nFeature_SCT SCT_snn_res.0.3 seurat_clusters SCT_snn_res.0.5 SCT_snn_res.0.4
Ctrl-2:AAACCCAGTCATAAAGx 2225 2 9 8 4
Ctrl-2:AAAGTGAAGAAGTGTTx 2135 2 9 8 4
Ctrl-2:AACAAGAAGTTGAAGTx 1552 5 8 6 4
Ctrl-2:AACAGGGCATATCTCTx 3106 2 9 8 0
Ctrl-2:AACGAAACAAGACCTTx 1813 2 8 6 4
Ctrl-2:AAGAACAGTGAGTGACx 1945 5 8 6 4
Annotation SCT_snn_res.1 SCT_snn_res.0.6
Ctrl-2:AAACCCAGTCATAAAGx Monocyte-derived infiltrating macrophage 10 9
Ctrl-2:AAAGTGAAGAAGTGTTx Monocyte-derived infiltrating macrophage 5 9
Ctrl-2:AACAAGAAGTTGAAGTx Weakly activated M1 9 8
Ctrl-2:AACAGGGCATATCTCTx Monocyte-derived infiltrating macrophage 16 9
Ctrl-2:AACGAAACAAGACCTTx Weakly activated M1 9 8
Ctrl-2:AAGAACAGTGAGTGACx Weakly activated M1 9 8
Annotation_2 Pseudotime
Ctrl-2:AAACCCAGTCATAAAGx Infiltrating macrophage 28.72614
Ctrl-2:AAAGTGAAGAAGTGTTx Infiltrating macrophage 27.32175
Ctrl-2:AACAAGAAGTTGAAGTx Transitional M1 16.39644
Ctrl-2:AACAGGGCATATCTCTx Infiltrating macrophage 31.99751
Ctrl-2:AACGAAACAAGACCTTx Transitional M1 26.84087
Ctrl-2:AAGAACAGTGAGTGACx Transitional M1 27.36084
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()
