很多刚刚入门细胞组学数据分析的小伙伴常常会问，单细胞数据分析流程有没有所谓范式？从拿到表达矩阵开始，有没有一个标准的技术路线以供参考？对于这样的问题，不妨从已发表的单细胞转录组文章中寻找答案。
今天开始我们将深入探讨一篇发表于《Cell Stem Cell》的关于人类皮肤创伤愈合的单细胞/空间转录组文章——《Spatiotemporal single-cell roadmap of human skin wound healing》。这篇文章通过单细胞转录组和空间转录组技术，揭示了损伤愈合过程中细胞和分子的动态变化。我们将通过代码复现的方式，一步步解析数据处理的流程，并展示如何从原始数据中提取有价值的信息。无论你是生物信息学新手还是经验丰富的研究者，相信此系列推文都能为你带来新的启发！
文章使用的数据包括人和小鼠的单细胞转录组数据、人的空间转录组数据、以及多个已发表的人皮肤单细胞转录组公开数据，具体如下：

• 损伤愈合组：对5名健康志愿者，使用4毫米活检针在每位志愿者的臀部上方区域创建3个全层皮肤伤口。在创伤后的第1天、第7天和第30天，使用6毫米活检针采集伤口边缘组织。样本用于单细胞RNA测序（10X Genomics）（GSE241132）和/或空间转录组学分析（Visium）（GSE241124）。
• 静脉溃疡组（VU）：对9名捐赠者，包括静脉溃疡患者和年龄、身体部位匹配的健康对照者，使用4毫米活检针采集慢性伤口边缘组织。样本用于单细胞RNA测序（GSE265972）。

• 小鼠：使用8-10周龄的C57BL/6J野生型小鼠，使用4毫米活检针在背部皮肤上创建一个全层伤口，深度穿透至肉膜层。在创伤后的第3天，使用6毫米活检针采集伤口边缘组织以及距离伤口至少1.5厘米的完整皮肤组织。样本用于单细胞RNA测序（GSE218430）。
• 其他公开数据库数据：糖尿病足溃疡（DFU）皮肤单细胞RNA数据（GSE165816）；人皮肤单细胞RNA数据（E-MTAB-8142）；小鼠皮肤损伤单细胞RNA数据（GSA: CRA010641）。

这些数据都可以从公开数据库如GEO、EBI直接下载，下载到的数据基本为原始的表达矩阵，便于我们从头进行代码复现。
根据文章，作者首先对损伤愈合组的单细胞转录组和空间转录组数据进行了伪bulk分析，通过PCA分析和差异基因的GO富集分析验证了样本间的关系，并定义了取样的愈合阶段。

1. 单细胞转录组数据的PCA分析

首先，我们从GEO数据库下载GSE241132（sc）和GSE241124（st）下的所有数据。
将sc的所有数据存放至一个目录下（eg. /sc_data/），每个以样本名命名的目录下为cellranger输出的三个矩阵文件（matrix.mtx.gz、features.tsv.gz、barcodes.tsv.gz）（这里需要对下载的文件进行重命名，这一步可以在Linux下使用for循环实现，这里不做额外说明）。另外，在sc_data/下sample.list文件中写入所有样本名，每个样本一行，格式如下：

PWH26D0
PWH26D1
PWH26D30

使用Seurat包中的AggregateExpression函数将单细胞表达矩阵转换为pseudo bulk表达矩阵，并使用DESeq2包进行PCA分析。

library(Seurat)
library(DESeq2)
library(ggplot2)

###读入样本列表
l <- read.table('./sample.list', header=FALSE)
mlist <- list()

###数据读入
for (i in 1:length(l$V1)){
      count <- Read10X(paste0('/sc_data/',l$V1[i],'/'))
      obj <- CreateSeuratObject(count, assay='RNA', min.cells=3, min.features=200, project=l$V1[i])
      mlist[[i]] <- obj
}
a_merge <- merge(x = mlist[[1]], y = mlist[2:length(mlist)])

###以单样本为单位转为pseudo bulk矩阵
bulk <- AggregateExpression(a_merge, return.seurat = F, slot = "counts", assays = "RNA",
                                                group.by = "orig.ident")
mat <- bulk$RNA

###使用DESeq2进行PCA分析
col <- data.frame(row.names=colnames(mat),
                             donor=c(rep('Donor5',4),rep('Donor4',4),rep('Donor3',4)),
                             stage=c('D0','D1','D30','D7','D0','D1','D30','D7','D0','D1','D30','D7'))
dds <- DESeqDataSetFromMatrix(countData=mat, colData = col, design=~stage)
rlo <- rlog(dds, blind = FALSE)
pcaData <- plotPCA(rlo, intgroup=c("stage","donor"),returnData = T)

ggplot(pcaData, aes(x=PC1,y=PC2,color=stage, shape=donor)) + 
            geom_point(size=4) + theme_bw()

复现图

原文图

复现的PCA结果基本与原文一致。

2. GO富集分析

下面筛选每个取样阶段相对于D0的差异基因，参考文章中的阈值（logFC>1 & p-value < 0.05），取Top500差异基因，以Wound1为例。

dds1 <- DESeq(dds, fitType = 'mean', minReplicatesForReplace = 7, parallel = FALSE)
###D1
res_d1 <- results(dds1, contrast=c('stage', 'D1', 'D0'))
res1 <- data.frame(res_d1, stringsAsFactors = FALSE, check.names = FALSE)
res1 <- res1[order(res1$padj, res1$log2FoldChange, decreasing = c(FALSE, TRUE)), ]
res1_up<- res1[which(res1$log2FoldChange >= 1 & res1$padj < 0.05),]
write.table(head(rownames(res1_up),500), quote=F, row.names=F, col.names=F, file='Wound1.gene.list')

然后对Top500差异基因进行GO富集分析，同样以Wound1为例。

library(clusterProfiler)
library(org.Hs.eg.db)
model <- "org.Hs.eg.db"

genes <- read.table('Wound1.gene.list', header = FALSE)
genes$V1 <- as.character(genes$V1)
full <- bitr(genes$V1, fromType = "SYMBOL" ,toType = c("ENSEMBL","ENTREZID") , OrgDb = model)
go <- enrichGO(gene = full$ENTREZID, 
                          OrgDb = model, 
                          ont = "BP", 
                          pAdjustMethod = "fdr", 
                          pvalueCutoff = 0.05, readable = TRUE)
head(go,5)[,2]
[1] "ribonucleoprotein complex biogenesis"
[2] "cellular respiration"                
[3] "aerobic respiration"                 
[4] "oxidative phosphorylation"           
[5] "ATP metabolic process"

Wound1时期的差异基因富集的Top5 GO条目与核糖体蛋白、ATP代谢等相关，与文中基本一致。
那么如何生成文中的GO富集棒棒糖图呢？我们已对每个时期的DEG进行了GO富集分析并将结果Top5输出至csv文件中。使用ggplot2的geom_segment和geom_point绘制棒棒糖图，并使用grid包指定分面颜色。

library(ggplot2)
library(dplyr)
library(grid)

w1 <- read.csv('Wound1.go.csv')
w7 <- read.csv('Wound7.go.csv')
w30 <- read.csv('Wound30.go.csv')
w1$stage <- 'Wound1'
w7$stage <- 'Wound7'
w30$stage <- 'Wound30'

data <- rbind(w1,w7,w30)
data$logP <- -log10(data$p.adjust)

data$Description <- factor(data$Description, levels=unique(data$Description), ordered=T)
data$stage <- factor(data$stage, levels=c('Wound1', 'Wound7', 'Wound30'), ordered=T)
p <- ggplot(data, aes(x = Count, y = Description)) +
                    geom_segment(aes(yend=Description, x=0, xend=Count, color=logP)) + 
                    geom_point(size=4, aes(color=logP)) +
                    facet_wrap(~ stage, ncol = 3) +
                    scale_color_continuous(low="pink",high="red") +
                    labs(x = "Gene Count", y = "GO Term", title = "Top 5 GO Terms per Condition") +
                    theme_bw()

###设置分面颜色       
g <- ggplotGrob(p)
strip_bg <- g$grobs[grep("strip", g$layout$name)]
colors <- c("salmon", "skyblue", "lightgreen")
for (i in seq_along(strip_bg)) {
    strip_bg[[i]]$grobs[[1]]$children[[1]]$gp$fill <- colors[i]
}
g$grobs[grep("strip", g$layout$name)] <- strip_bg
grid.newpage()
grid.draw(g)

复现图

原文图

3. 空间转录组数据的PCA分析

我们使用类似的方法对空间样本也进行PCA分析。这里只额外说明数据读入，其他与上述一致。
将st的所有数据存放至一个目录下（eg. /st_data/），每个以样本名命名的目录下为spaceranger输出的矩阵文件（filtered_feature_bc_matrix.h5）和图像目录（spatial）。另外，在st_data/下st_sample.list文件中写入所有样本名，每个样本一行。

library(Seurat)

l <- read.table('./st_sample.list', header=FALSE)
mlist <- list()

for (i in 1:length(l$V1)){
    obj <- Load10X_Spatial(paste0('/st_data/',l$V1[i]), filename='/filtered_feature_bc_matrix.h5')
    obj$batch <- l$V1[i]
    mlist[[i]] <- obj
}

a_merge <- merge(x = mlist[[1]], y = mlist[2:length(mlist)])

bulk <- AggregateExpression(a_merge, return.seurat = F, slot = "counts", assays = "Spatial",
                                                group.by = "batch")
mat <- bulk$Spatial

至此，我们就完成了文章中对创伤愈合组的伪bulk分析。在下一篇，我们将复现单细胞转录组和空间转录组数据的预处理过程。