单细胞经典文章复现（二）

天意生信云

发布于 2025-03-06 21:53:30

6300

代码可运行

运行总次数：0

代码可运行

今天复现的文章是2022年发表在《Nature》子刊《Communications Biology》上，题为“Single-cell transcriptomics reveals cellular heterogeneity and molecular stratification of cervical cancer”的宫颈癌单细胞研究。

下载数据

本文中共有6个子宫样本，其中3个来自于健康对照，另外3个则来自于宫颈癌。编号为E-MTAB-11948，点击文章提供的链接可跳转文章数据库，下载数据。

作者上传到数据库中的数据已经是经过cellRanger处理过后得到的表达矩阵，每个样本有三个文件barcodes.tsv.gz 、features.tsv.gz、 matrix.mtx.gz。

数据读取，构建seurat对象

在基于R的seurat包中，我们可以完成单细胞数据的降维聚类和细胞注释等上游分析。本文中的单细胞分析中没有涉及到批次效应的问题，只需要将样本合并处理即可。

# 加载ArrayExpress包
library(ArrayExpress)

#options(timeout=600)
# 使用getAE()函数下载E-MTAB-11948数据
expression_data <- getAE("E-MTAB-11948")

# 查看数据
head(expression_data)#可能会因为超时，下载不完全，自检

在读取数据之前，需要更改一下文件的命名，并单独为每一个样本创建一个独立的文件夹。

## 创建每个样本文件夹

# 运行
ls *barcodes* |perl -ne 'chomp;/(.*)_(.*.gz)/;print"mkdir $1\n";' |sh

文件重命名
# 运行
ls *.gz |perl -ne 'chomp;/(.*)_.*L(.*gz)/;print"mv $_ $1/$2$3\n";' |sh

每个样本的文件夹都整理成下图这种：

#加载必须的包
library(Seurat)
library(dplyr)
library(patchwork)
library(RColorBrewer)


#读取每个样本的表达矩阵并构建seurat对象
Samples <- c("sample1", "sample2", "sample3","sample4","sample6")
Samples_list <- list()
for (i in Samples) {
  input_dir <- paste0("/mnt/data/home/sunting/2-single_cell_rna_seq/3-fuxianwenzhang/1/0-data/", i , "/") 
  Samples_list[[i]] <- Read10X(data.dir = input_dir) 
  Samples_list[[i]] <- CreateSeuratObject(Samples_list[[i]], project = i)  
}

在运行过程中，读取sample5这个样本的数据文件时发现作者误传了两边barcodes.tsv.gz。Read10x在读取时不能识别对应基因名。所以在后续的分析中先删掉了sample5这个样本。

正式开始单细胞的上游分析

#将每个样本的seurat对象合并成一个整体
#在本次单细胞分析中没有涉及到批次效应的问题
sc_Merge <- merge(Samples_list[[1]], Samples_list[-1], add.cell.ids = names(Samples_list))

#计算每个细胞中线粒体基因的占比
sc_Merge[["percent.mt"]] <- PercentageFeatureSet(sc_Merge, pattern = "^MT-")
#可以查看都有哪些线粒体基因
rownames(sc_Merge)[grep("^MT-", rownames(sc_Merge))]
#在处理自己的数据的时候，线粒体基因可能没有共同特征，那我们可以去gtf或者gff文件中查看，然后把pattern展开写

#质控，过滤掉count数小于200，基因数小于200，线粒体基因占比大于10的细胞
sc_Merge <- subset(sc_Merge, subset = nFeature_RNA > 200 & nCount_RNA > 200 & percent.mt < 10)


#对数据进行标准化处理
#因为细胞间的异质性除了可能来自于细胞自身基因的差异表达，还有可能是由测序技术误差造成的。数据标准化的目的便是排除技术误差，让测序深度和文库大小不同的细胞的基因表达量具有可比性，进而更好地呈现基因水平的差异。
sc_Merge <- NormalizeData(sc_Merge, normalization.method = "LogNormalize", scale.factor = 10000)

#寻找高变基因
sc_Merge <- FindVariableFeatures(sc_Merge, selection.method = "vst", nfeatures = 2000)


all.genes <- rownames(sc_Merge)

#对所有基因的表达量进行归一化处理
sc_Merge <- ScaleData(sc_Merge, features = all.genes)
#-----------节省时间的做法---------
#ScaleData()默认使用的是先前鉴定的2000个高变基因。所以将features参数去掉，将all.genes变成高变基因即可减少运行时间
#sc_Merge <- ScaleData(sc_Merge, features = scale.genes)


#根据前面计算的2000个高变基因进行主成分分析
sc_Merge <- RunPCA(sc_Merge, features = VariableFeatures(object = sc_Merge))


#Seurat 应用了一种基于图表的聚类方法。重要的是，驱动聚类分析的距离度量(基于先前识别的 PC)保持不变。将细胞嵌入到图形结构中-例如 K 最近邻(KNN)图，在具有相似特征表达模式的细胞之间绘制边界，然后将该图划分为高度相互关联的不同簇。
#选取前10个主成分进行降维聚类分析
sc_Merge <- FindNeighbors(sc_Merge, dims = 1:10)
sc_Merge <- FindClusters(sc_Merge, resolution = 0.5)
sc_Merge <- RunTSNE(sc_Merge, dims = 1:10)

DimPlot(sc_Merge, reduction = "tsne")
#根据文章中给定的markers基因进行细胞类型注释， 现在单细胞分析各个组织中的分析都已经很丰富，
#一些典型的marker基因都可以经过文献查阅找到，因此在做自己的单细胞分析的时候不用担心无法确定marker基因
FeaturePlot(sc_Merge, features = c("CDH1","PECAM1","COL1A2","ACTA2","CD3E","CD79A","LYZ","CSF3R"),reduction = "tsne")

#知道哪些细胞是属于什么类型后，我们开始对这些细胞进行重命名
new.cluster.ids = c("Endothelial", "Epithelial", "Epithelial", "Smooth muscle cells", "Smooth muscle cells",
                    "Lymphocytes", "Fibroblast", "Fibroblast", "Epithelial", "Epithelial", "Epithelial",
                    "Fibroblast", "Neutrophils", "Macrophage", "Epithelial", "Lymphocytes", "Lymphocytes",
                    "Fibroblast", "Endothelial", "Epithelial", "Fibroblast", "Endothelial","Myocyte","Stem cells")


names(new.cluster.ids)=levels(sc_Merge)
sc_Merge <- RenameIdents(sc_Merge,new.cluster.ids)
sc_Merge$celltype=Idents( sc_Merge)
cols = c("#80BD7D", "#BB436F", "#9D4D3C", "#3666A7", "#C5B7D9", "#CF5D15", "#2E8339", "#F0A202", "#D84B16")
DimPlot(sc_Merge, reduction = "tsne", label = TRUE, pt.size = 0.5,cols = cols)

#原文中已经给出的celltype的markers基因
Idents(sc_Merge)=sc_Merge$seurat_clusters

Markers=c("CDKN2A","EPCAM","CD24","CDH1",#epithelial
          "PECAM1","CDH5","ENG",#endothelial
          "COL1A2","DCN","APOD",#fibroblast
          "ACTA2","ACTG2",#smmoth muscle cells
          "CD3E","CD3D","CD2",#T cells
          "CD79A","CD19",#B  cells
          "CD68","CD163","LYZ",#macrophage          
          "CSF3R")#neutrophils                                        

Markers.avg=AverageExpression(sc_Merge, features = Markers)

Markers.avg$RNA=Markers.avg$RNA[,c(13,14,16,6,17,4,5,7,8,21,12,18,1,22,19,10,3,9,15,20,2,11)]
group=sc_Merge@meta.data[,c(6,7)] %>% unique()


rownames(group)=group$seurat_clusters
group=group[c(13,14,16,6,17,4,5,7,8,21,12,18,1,22,19,10,3,9,15,20,2,11),]


colnames(Markers.avg$RNA) <- gsub("g", "", colnames(Markers.avg$RNA))
bk <- seq(-2, 2, length.out = 100)

# 使用pheatmap绘制热图，确保color参数正确生成
pheatmap::pheatmap(Markers.avg$RNA, scale = "row",  
                   cluster_rows = F, cluster_cols = F,  
                   show_colnames = T, breaks = bk, annotation_col = group, 
                   color = colorRampPalette(rev(brewer.pal(n = 11, name = "PiYG")))(length(bk)))

#画细胞数量在不同组中的分布，需要将样本的分组信息添加近seurat对象中的metadata中
samples.group=data.frame(orig.ident=c("sample1","sample2","sample3","sample4","sample6"),                        
                         group=c(rep("Tumor",3), rep("Control",2)))

sc_Merge@meta.data$Barcode=rownames(sc_Merge@meta.data)
sc_Merge@meta.data=left_join(sc_Merge@meta.data, samples.group, by="orig.ident")
rownames(sc_Merge@meta.data)=sc_Merge@meta.data$Barcode


#统计细胞类型在Tumor组和Control组中的细胞数量
celltype.number=sc_Merge@meta.data[,c(1,7,9)] %>% 
  group_by(orig.ident, celltype, group) %>%                         
  summarize(count=n())
celltype.number$group=factor(celltype.number$group, levels = c("Tumor","Control"))

#利用ggplot2画柱状图展示
p4 <- ggplot2::ggplot(celltype.number,aes(x=celltype, y=count, fill=group))+  
  geom_bar(stat="identity",position = "stack")+  
  theme_bw()+coord_flip()+  
  scale_fill_manual(values = c("#9562A6","#91BB63"))+  
  ylab("Cell number")+xlab("")
ggsave("celltype_bar.pdf", plot = p4, path = "/mnt/data/home/sunting/2-single_cell_rna_seq/3-fuxianwenzhang/1/1-result", width = 10, height = 8)

#提取数据集中的Epithelial, 后面针对Epithelial的分析其实和前面的步骤相同
Idents(sc_Merge)=sc_Merge$celltype
Epithelial=subset(sc_Merge, idents = "Epithelial")
Epithelial <- FindVariableFeatures(Epithelial, selection.method = "vst", nfeatures = 2000)
Epithelial <- RunPCA(Epithelial, features = VariableFeatures(object = Epithelial))
ElbowPlot(Epithelial)
#选取前10个主成分进行降维聚类分析
Epithelial <- FindNeighbors(Epithelial, dims = 1:8)
Epithelial <- FindClusters(Epithelial, resolution = 0.2)
Epithelial <- RunTSNE(Epithelial, dims = 1:8)
DimPlot(Epithelial, reduction = "tsne",label = T, cols=cols)
DimPlot(Epithelial, reduction = "tsne",label = T, group.by = "group", cols = c("#91BB63","#9463A7"))

All.Markers=FindAllMarkers(Epithelial, only.pos = T, logfc.threshold = 0.25)

TOP.genes=All.Markers %>% group_by(cluster) %>% top_n(n=20, wt=avg_log2FC)

#画TOP20的特异性表达的差异基因
TOP.genes.avg=AverageExpression(Epithelial, features = unique(TOP.genes$gene))
bk <- c(seq(-2,2,by=0.01))
pheatmap::pheatmap(TOP.genes.avg$RNA, scale = "row",                    
                   cluster_rows = F, cluster_cols = F,                    
                   show_colnames = T,breaks=bk,                   
                   color =colorRampPalette(rev(brewer.pal(n = 11, name = "PiYG")))(length(bk)))