给学徒们收集整理了几套带GitHub源代码的文献图表合辑，让优秀者一点一滴拆解开来分享给大家。（全部的代码复制粘贴即可运行，欢迎尝试以及批评指正）
现在是雪貂支气管肺泡灌洗液单细胞转录组显示SARS-CoV-2感染期间巨噬细胞的顺序变化专辑第4讲：主要是巨噬细胞内部细分亚群的拟时序和RNA速率分析
下面是前年实习生（日行一膳）的分享
1643459432584
本次复现的是于2021年7月27日发表在Nature Communications上的”Single-cell transcriptome of bronchoalveolar lavage fluid reveals sequential change of macrophages during SARS-CoV-2 infection in ferrets“中的Figure4
Fig. 4 Trajectory analysis from MDIM to M1 and M2 macrophages.
1643729136434
代码语言：javascript复制## 安装所需要的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
1644047189899
1644067694276
代码语言：javascript复制## 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()
1644048057620
代码语言：javascript复制## 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
1644049309001
1644068852744
代码语言：javascript复制## 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()
这样的基础认知，也可以看单细胞的基础10讲：
- 01. 上游分析流程
- 02.课题多少个样品，测序数据量如何
- 03. 过滤不合格细胞和基因（数据质控很重要）
- 04. 过滤线粒体核糖体基因
- 05. 去除细胞效应和基因效应
- 06.单细胞转录组数据的降维聚类分群
- 07.单细胞转录组数据处理之细胞亚群注释
- 08.把拿到的亚群进行更细致的分群
- 09.单细胞转录组数据处理之细胞亚群比例比较
