小样本量单细胞肿瘤数据分析模范 / 四六文摘

考虑到咱们生信技能树粉丝对单细胞数据挖掘的需求，我开通了一个专栏《100个单细胞转录组数据降维聚类分群图表复现》，也亲自示范了几个，不过自己带娃，读博，时间精力有限，所以把剩余的90多个任务安排了学徒，实习生，学员。真的是太棒了，群策群力！

下面是2021第2期生信学习班优秀学员的投稿

Step0:加载需要的R包
step1:导入数据
step2:QC质控
step3:合并2个数据集(CCA)
step4:分cluster
step5:细胞类型注释(根据marker gene)
step6:细胞类型注释(基于singleR)
step7: 最终细胞类型的差异基因
step8: 看不同细胞亚群比例变化

这次我们要复现的单细胞数据来自2020年10月8日发表的一篇单细胞相关的文章——单细胞测序首次应用于早期前列腺癌的诊断和分层

代码来自于jimmy老师的课程，该文章数据链接GSE157703

Step0:加载需要的R包

library(Seurat) library(ggplot2) library(clustree) library(cowplot) library(dplyr)

step1:导入数据

######## step1:导入数据 ######## rm(list=ls()) options(stringsAsFactors = F) ### 读取txt格式的矩阵 library(data.table) pca1 <- fread("GSM4773521_PCa1_gene_counts_matrix.txt.gz", data.table = F) #这里的Warning message不要怕，是因为读取的时候把行名当成了一列，所以多出来了一列'V1’,后面会删掉这里一列再给数据重新添加行名 pca1[1:4,1:4]

pca2 <- fread("GSM4773522_PCa2_gene_counts_matrix.txt.gz", data.table = F) pca2[1:4,1:4] d1=pca1[,-1] rownames(d1)=pca1[,1]

创建SeuratObject

#分别创建SeuratObject scRNA1 <- CreateSeuratObject(counts = d1, #和文章一样，初步质控 min.cells = 3, #在不少于三个细胞中表达 min.features = 200, #基因数不少于200 project = "pca1") #24576 features & 2896 samples d2=pca2[,-1] rownames(d2)=pca2[,1] scRNA2 <- CreateSeuratObject(counts = d2, min.cells = 3, min.features = 200, project = "pca2") #24399 features & 4931 samples

merge函数整合两样本

#利用merge函数把两个样本合在一起，并且给每个样本都加上各自的样本名，方便后续分析 sce.all = merge(scRNA1, y = scRNA2, add.cell.ids = c("pca1", "pca2"), project = "PCA", merge.data = TRUE) #25948 features & 7827 samples


as.data.frame(sce.all@assays$RNA@counts[1:10, 1:2])

#可以看到每个细胞都被写上样本名了

#head(sce.all@meta.data, 10)

#orig.ident nCount_RNA nFeature_RNA

#pca1_AAACCTGAGCGTTCCG_1       pca1       5608         1823

#pca1_AAACCTGCACACTGCG_1       pca1       7860         2336

#pca1_AAACCTGGTCGCTTTC_1       pca1        607          382

table(sce.all@meta.data$orig.ident) head(sce.all@meta.data) saveRDS(sce.all,file="sce.all_raw.rds")

step2:QC质控

过滤低质量细胞/基因

rm(list=ls()) ######## step2:QC质控 ###### dir.create("./1-QC") setwd("./1-QC") sce.all=readRDS("../sce.all_raw.rds") #计算线粒体基因比例 ##人和鼠的基因名字稍微不一样。（鼠是mt） mito_genes=rownames(sce.all)[grep("^MT-", rownames(sce.all))] mito_genes #13个线粒体基因 sce.all=PercentageFeatureSet(sce.all, "^MT-", col.name = "percent_mito") fivenum(sce.all@meta.data$percent_mito) #计算核糖体基因比例 ribo_genes=rownames(sce.all)[grep("^Rp[sl]", rownames(sce.all),ignore.case = T)] ribo_genes sce.all=PercentageFeatureSet(sce.all, "^RP[SL]", col.name = "percent_ribo") fivenum(sce.all@meta.data$percent_ribo) #计算红血细胞基因比例 rownames(sce.all)[grep("^Hb[^(p)]", rownames(sce.all),ignore.case = T)] sce.all=PercentageFeatureSet(sce.all, "^HB[^(P)]", col.name = "percent_hb") fivenum(sce.all@meta.data$percent_hb)


#可视化细胞的上述比例情况

feats <- c("nFeature_RNA", "nCount_RNA", "percent_mito", "percent_ribo", "percent_hb")

feats <- c("nFeature_RNA", "nCount_RNA")

p1=VlnPlot(sce.all, group.by = "orig.ident", features = feats, pt.size = 0.01, ncol = 2) + 

    NoLegend()

p1

library(ggplot2)

ggsave(filename="Vlnplot1.pdf",plot=p1)

ggsave(filename="Vlnplot1.png",plot=p1)
feats <- c("percent_mito", "percent_ribo", "percent_hb")

p2=VlnPlot(sce.all, group.by = "orig.ident", features = feats, pt.size = 0.01, ncol = 3, same.y.lims=T) + 

    scale_y_continuous(breaks=seq(0, 100, 5)) +

 NoLegend()

p2 

ggsave(filename="Vlnplot2.pdf",plot=p2)

ggsave(filename="Vlnplot2.png",plot=p2)

#用FeatureScatter函数来可视化特征-特征之间的关系 p3=FeatureScatter(sce.all, "nCount_RNA", "nFeature_RNA", group.by = "orig.ident", pt.size = 0.5) ggsave(filename="Scatterplot.pdf",plot=p3) ggsave(filename="Scatterplot.png",plot=p3)

过滤前的"nFeature_RNA", "nCount_RNA", "percent_mito", "percent_ribo", "percent_hb"比例情况

过滤指标1:最少表达基因数的细胞&最少表达细胞数的基因

#过滤指标1:最少表达基因数的细胞&最少表达细胞数的基因 selected_c <- WhichCells(sce.all, expression = nFeature_RNA > 300) selected_f <- rownames(sce.all)[Matrix::rowSums(sce.all@assays$RNA@counts > 0 ) > 3] sce.all.filt <- subset(sce.all, features = selected_f, cells = selected_c) dim(sce.all) #过滤前 25948 7827 dim(sce.all.filt) #过滤后 24653 7459 #可以看到，主要是过滤了基因，其次才是细胞

过滤指标2:线粒体/核糖体基因比例

#过滤指标2:线粒体/核糖体基因比例(根据上面的violin图) selected_mito <- WhichCells(sce.all.filt, expression = percent_mito < 15) selected_ribo <- WhichCells(sce.all.filt, expression = percent_ribo > 3) selected_hb <- WhichCells(sce.all.filt, expression = percent_hb < 0.1) length(selected_hb) length(selected_ribo) length(selected_mito) sce.all.filt <- subset(sce.all.filt, cells = selected_mito) sce.all.filt <- subset(sce.all.filt, cells = selected_ribo) #sce.all.filt <- subset(sce.all.filt, cells = selected_hb) dim(sce.all.filt) table(sce.all.filt$orig.ident)


#可视化过滤后的情况

feats <- c("nFeature_RNA", "nCount_RNA")

p1_filtered=VlnPlot(sce.all.filt, group.by = "orig.ident", features = feats, pt.size = 0.1, ncol = 2) + 

    NoLegend()

ggsave(filename="Vlnplot1_filtered.pdf",plot=p1_filtered)

ggsave(filename="Vlnplot1_filtered.png",plot=p1_filtered)

feats <- c("percent_mito", "percent_ribo", "percent_hb") p2_filtered=VlnPlot(sce.all.filt, group.by = "orig.ident", features = feats, pt.size = 0.1, ncol = 3) + NoLegend() ggsave(filename="Vlnplot2_filtered.pdf",plot=p2_filtered) ggsave(filename="Vlnplot2_filtered.png",plot=p2_filtered) dim(sce.all.filt) ### [1] 24653 7086

过滤后的"nFeature_RNA", "nCount_RNA", "percent_mito", "percent_ribo", "percent_hb"比例情况

过滤指标3:过滤特定基因

#过滤指标3:过滤特定基因 ### Filter MALAT1 管家基因 sce.all.filt <- sce.all.filt[!grepl("MALAT1", rownames(sce.all.filt),ignore.case = T), ] ### Filter Mitocondrial 线粒体基因 sce.all.filt <- sce.all.filt[!grepl("^MT-", rownames(sce.all.filt),ignore.case = T), ] ### 当然，还可以过滤更多 dim(sce.all.filt) ### [1] 24639 7086

细胞周期评分

#细胞周期评分 ##标准化数据 sce.all.filt = NormalizeData(sce.all.filt) s.genes=Seurat::cc.genes.updated.2019$s.genes g2m.genes=Seurat::cc.genes.updated.2019$g2m.genes sce.all.filt=CellCycleScoring(object = sce.all.filt, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE) p4=VlnPlot(sce.all.filt, features = c("S.Score", "G2M.Score"), group.by = "orig.ident", ncol = 2, pt.size = 0.1) ggsave(filename="Vlnplot4_cycle.pdf",plot=p4) ggsave(filename="Vlnplot4_cycle.png",plot=p4)

sce.all.filt@meta.data %>% ggplot(aes(S.Score,G2M.Score))+geom_point(aes(color=Phase))+ theme_minimal() ggsave(filename="cycle_details.pdf" ) ggsave(filename="cycle_details.png" ) ### S.Score较高的为S期，G2M.Score较高的为G2M期，都比较低的为G1期

过滤doublets

###### 检测doublets #Doublets被定义为在相同细胞barcode下测序的两个细胞（例如被捕获在同一液滴中） sce.all.filt = FindVariableFeatures(sce.all.filt) sce.all.filt = ScaleData(sce.all.filt, vars.to.regress = c("nFeature_RNA", "percent_mito")) sce.all.filt = RunPCA(sce.all.filt, npcs = 20) sce.all.filt = RunTSNE(sce.all.filt, npcs = 20) sce.all.filt = RunUMAP(sce.all.filt, dims = 1:10) ### define the expected number of doublet cellscells. nExp <- round(ncol(sce.all.filt) * 0.04) ### expect 4% doublets ### remotes::install_github('chris-mcginnis-ucsf/DoubletFinder') library(DoubletFinder) #这个DoubletFinder包的输入是经过预处理（包括归一化、降维，但不一定要聚类）的 Seurat 对象 sce.all.filt <- doubletFinder_v3(sce.all.filt, pN = 0.25, pK = 0.09, nExp = nExp, PCs = 1:10)


###### 找Doublets
### DF的名字是不固定的，因此从sce.all.filt@meta.data列名中提取比较保险

DF.name = colnames(sce.all.filt@meta.data)[grepl("DF.classification", colnames(sce.all.filt@meta.data))]

#(图右)

p5.dimplot=cowplot::plot_grid(ncol = 2, DimPlot(sce.all.filt, group.by = "orig.ident") + NoAxes(), 

    DimPlot(sce.all.filt, group.by = DF.name) + NoAxes())

ggsave(filename="doublet_dimplot.pdf",plot=p5.dimplot)

ggsave(filename="doublet_dimplot.png",plot=p5.dimplot)
#(图左)红色部分上的黑点就是doublets

p5.vlnplot=VlnPlot(sce.all.filt, features = "nFeature_RNA", group.by = DF.name, pt.size = 0.1)

ggsave(filename="doublet_vlnplot.pdf",plot=p5.vlnplot)

ggsave(filename="doublet_vlnplot.png",plot=p5.vlnplot)
###### 过滤doublet

sce.all.filt=sce.all.filt[, sce.all.filt@meta.data[, DF.name] == "Singlet"]

#过滤到此结束

#最后过滤结果

dim(sce.all.filt) 

#[1] 24639  6803

dim(sce.all) 

#[1] 25948  7827

saveRDS(sce.all.filt, "sce.all_qc.rds")

### rm(list=ls())

#dev.off()

#上面我们用merge函数合并两个数据集姐，批次效应不明显的话是可以的，如果样本量比较大，并且有批次效应，或者有跨测序平台的情况，可以使用CCA进行整合

在看到jimmy老师代码之前，我只知道有doublet cells 这回事，但是很少在普通文章中见到过滤它们的，每天都在涨知识，开心 0v0

step3:合并2个数据集(CCA)

######## step3:合并2个数据集######## setwd('../') dir.create("2-int") getwd() setwd("2-int") #sce.all=readRDS("../1-QC/sce.all_qc.rds") sce.all=sce.all.filt sce.all


### 拆分为2 个seurat子对象(分开pca1和pca2)

table(sce.all@meta.data$orig.ident)

sce.all.list <- SplitObject(sce.all, split.by = "orig.ident")

sce.all.list
### 对两个数据进行处理

#鉴定高度变化的基因:用FindVariableFeatures函数来执行。默认情况下，每个dataset返回2,000个基因，这些基因将用于下游分析，如PCA。

for (i in 1:length(sce.all.list)) {

 print(i)

    sce.all.list[[i]] <- NormalizeData(sce.all.list[[i]], verbose = FALSE)

    sce.all.list[[i]] <- FindVariableFeatures(sce.all.list[[i]], selection.method = "vst", 

        nfeatures = 2000, verbose = FALSE)

}

#接下来的代码耗时较久，耐心等待哦。 #利用Seurat v3里面的CCA对数据进行整合,（里面的具体原理对我来说有难度，所以我打算遵循小洁老师的教诲：会用就行）。单样本做了均一化后，进行多样本的整合,首先是用 FindIntegrationAnchors 函数找到数据集间的锚点，然后用 IntegrateData 函数整合多个数据集。注意这两个函数的参数要一致. alldata.anchors <- FindIntegrationAnchors(object.list = sce.all.list, dims = 1:30, reduction = "cca") sce.all.int <- IntegrateData(anchorset = alldata.anchors, dims = 1:30, new.assay.name = "CCA") #运行IntegrateData之后，Seurat对象将包含一个具有整合（或“批次校正”）表达矩阵的新Assay。 #请注意，原始值（未校正的值）仍存储在“RNA”分析的对象中，因此可以来回切换。然后可以使用这个新的整合矩阵进行下游分析和可视化 names(sce.all.int@assays) #[1] "RNA" "CCA" sce.all.int@active.assay #[1] "CCA"

CCA数据整合——可视化（整合前&整合后）

#### 运行标准流程并进行可视化（合并前和合并后） p1.compare=plot_grid(ncol = 3, DimPlot(sce.all, reduction = "pca", group.by = "orig.ident")+NoAxes()+ggtitle("PCA raw_data"), DimPlot(sce.all, reduction = "tsne", group.by = "orig.ident")+NoAxes()+ggtitle("tSNE raw_data"), DimPlot(sce.all, reduction = "umap", group.by = "orig.ident")+NoAxes()+ggtitle("UMAP raw_data"), DimPlot(sce.all.int, reduction = "pca", group.by = "orig.ident")+NoAxes()+ggtitle("PCA integrated"), DimPlot(sce.all.int, reduction = "tsne", group.by = "orig.ident")+NoAxes()+ggtitle("tSNE integrated"), DimPlot(sce.all.int, reduction = "umap", group.by = "orig.ident")+NoAxes()+ggtitle("UMAP integrated") ) p1.compare ggsave(plot=p1.compare,filename="Before&After_int.pdf",width = 10, height = 7) ggsave(plot=p1.compare,filename="Before&After_int.png",width = 10, height = 7)

虽然整合前的样本之间的差异也不明显，但是也是存在一些的，使用CCA整合后可以看到样本间差异基本可以忽略了。

step4:分cluster

######## step4:分cluster###### sce.all.int #After sce.all #Before sce.all=sce.all.int sce.all@active.assay


#使用Seurat包中的FindNeighbors函数计算构建SNN图。

sce.all=FindNeighbors(sce.all, dims = 1:30, k.param = 60, prune.SNN = 1/15)
#先执行不同resolution 下的分群

#设置不同的分辨率，观察分群效果(选择哪一个？)

#在这里，jimmy老师把常用参数全部循环了一遍，再从中挑选最合适的

for (res in c(0.01, 0.05, 0.1, 0.2, 0.3, 0.5,0.8,1)) {

    sce.all=FindClusters(sce.all, graph.name = "CCA_snn", resolution = res, algorithm = 1)

}

apply(sce.all@meta.data[,grep("CCA_snn_res",colnames(sce.all@meta.data))],2,table)
#可以看出不同分辨率设置对于分群数量来说影响挺大，之前这里的分辨率我都没改过。（表情：0_0）
#低分辨率的分群情况。（0.01，0.1，0.2）

p1_dim=plot_grid(ncol = 3, DimPlot(sce.all, reduction = "umap", group.by = "CCA_snn_res.0.01") + 

    ggtitle("louvain_0.01"), DimPlot(sce.all, reduction = "umap", group.by = "CCA_snn_res.0.1") + 

    ggtitle("louvain_0.1"), DimPlot(sce.all, reduction = "umap", group.by = "CCA_snn_res.0.2") + 

    ggtitle("louvain_0.2"))

ggsave(plot=p1_dim, filename="Dimplot_diff_resolution_low.pdf",width = 14)

ggsave(plot=p1_dim, filename="Dimplot_diff_resolution_low.png",width = 14)
#高分辨率的分群情况。（0.3，0.8，1）

p1_dim=plot_grid(ncol = 3, DimPlot(sce.all, reduction = "umap", group.by = "CCA_snn_res.0.8") + 

                   ggtitle("louvain_0.8"), DimPlot(sce.all, reduction = "umap", group.by = "CCA_snn_res.1") + 

                   ggtitle("louvain_1"), DimPlot(sce.all, reduction = "umap", group.by = "CCA_snn_res.0.3") + 

                   ggtitle("louvain_0.3"))

ggsave(plot=p1_dim, filename="Dimplot_diff_resolution_high.pdf",width = 18)

ggsave(plot=p1_dim, filename="Dimplot_diff_resolution_high.png",width = 18)
#使用clustree包来可视化不同分辨率下细胞在聚类群之间的分配。

### 参考：https://mp.weixin.qq.com/s/WRhMC3Ojy1GWYfLS_4vSeA

p2_tree=clustree(sce.all@meta.data, prefix = "CCA_snn_res.")

ggsave(plot=p2_tree, filename="Tree_diff_resolution.pdf",width = 8)

ggsave(plot=p2_tree, filename="Tree_diff_resolution.png",width = 8)

低分辨率的分群情况。（0.01，0.1，0.2）

高分辨率的分群情况。（0.3，0.8，1）

`clustree`包可视化不同分辨率下细胞在聚类群之间的分配

上图不同的颜色表示不同的 resolution ，点越大表示亚群中包含的细胞数目越多。我们只知道随着 resolution 值越来越大，分的亚群的也来越多，但是通上图，我们可以观察到随着resolution 增大，具体是哪些细胞亚群发生了变化。也可以看到参数设置为0.8或1时，都分成了16个cluster，不会再随着分辨率的增大而分出新的cluster了。所以接下来分析，按照分辨率为0.8进行。

#接下来分析，按照分辨率为0.8进行 sel.clust = "CCA_snn_res.0.8" sce.all <- SetIdent(sce.all, value = sel.clust) table(sce.all@active.ident) saveRDS(sce.all, "sce.all_int.rds") rm(list=ls())

小结：

一分耕耘一分收获，跑一遍老师的代码，感觉自己这个井底之蛙又往井口挪了挪。

写这个笔记前，我连在Typora里面插入图片都不会操作，而现在，我完成了这个笔记，小小的进步就是每天的兴奋剂啦。

感恩老师们的辛苦付出，冲冲冲！！！

上半部分我们对数据进行了一顿操作，接下来该到了细胞类型注释这一步啦

step5:细胞类型注释(根据marker gene)

######## step5:细胞类型注释(根据marker gene)###### ### setwd("/home/data/gma49/scRNA") dir.create("../3-cell") setwd("../3-cell")


sce.all=readRDS( "sce.all_int.rds")

table(sce.all@meta.data$seurat_clusters)

table(sce.all@meta.data$CCA_snn_res.0.8)

#展示按照不同标准进行分群的结果，并没有什么区别 p1=DimPlot(sce.all, reduction = "umap", group.by = "seurat_clusters",label = T) p2=DimPlot(sce.all, reduction = "umap", group.by = "CCA_snn_res.0.8",label = T) ggsave("umap_by_CCA_snn_res.0.8.pdf", plot = p1, width = 8, height = 7) ggsave("umap_by_CCA_snn_res.0.8.png", plot = p1, width = 8, height = 7) ggsave("umap_by_seurat_clusters.pdf", plot = p2, width = 8, height = 7) ggsave("umap_by_seurat_clusters.png", plot = p2, width = 8, height = 7)

各个细胞类型的marker

#各个细胞类型的marker。（这里就太考验生物学背景知识了，而我的生物学背景知识比纸薄0.0，之前都是无脑用singleR包注释） #jimmy老师非常贴心的给列出来了，感恩！ ### T Cells (CD3D, CD3E, CD8A), ### B cells (CD19, CD79A, MS4A1 [CD20]), ### Plasma cells (IGHG1, MZB1, SDC1, CD79A), ### Monocytes and macrophages (CD68, CD163, CD14), ### NK Cells (FGFBP2, FCG3RA, CX3CR1), ### Photoreceptor cells (RCVRN), ### Fibroblasts (FGF7, MME), ### Endothelial cells (PECAM1, VWF). ### epi or tumor (EPCAM, KRT19, PROM1, ALDH1A1, CD24). ### immune (CD45+,PTPRC), epithelial/cancer (EpCAM+,EPCAM), ### stromal (CD10+,MME,fibo or CD31+,PECAM1,endo)


library(ggplot2) 

#查看这些marker基因在各个cluster里的表达情况

genes_to_check = c('PTPRC', 'CD3D', 'CD3E', 'CD4','CD8A','CD19', 'CD79A', 'MS4A1' ,

                   'IGHG1', 'MZB1', 'SDC1',

                   'CD68', 'CD163', 'CD14', 

                   'TPSAB1' , 'TPSB2',  ### mast cells,

                   'RCVRN','FPR1' , 'ITGAM' ,

                   'FGF7','MME', 'ACTA2',

                   'PECAM1', 'VWF', 

                   'EPCAM' , 'KRT19', 'PROM1', 'ALDH1A1' )

library(stringr)  

#sce.all@assays$RNA

p_all_markers <- DotPlot(sce.all, features = genes_to_check,

             assay='RNA'  )  + coord_flip()

p_all_markers

ggsave(plot=p_all_markers, filename="check_all_marker_by_seurat_cluster.pdf",width = 8, height = 6)

ggsave(plot=p_all_markers, filename="check_all_marker_by_seurat_cluster.png",width = 8, height= 6)

用常见的marker先区分一下，大概能知道cluster对应什么类型。这里用到的marker都是在文献里面经常出现的。

#初步定义后，心里大概有一个了解，再利用自己的相关知识开始细分


### 1、Tcells（T细胞）

genes_to_check = c('PTPRC', 'CD3D', 'CD3E', 'CD4','CD8A',

                   'CCR7', 'SELL' , 'TCF7','CXCR6' , 'ITGA1',

                   'FOXP3', 'IL2RA',  'CTLA4','GZMB', 'GZMK','CCL5',

                   'IFNG', 'CCL4', 'CCL3' ,

                   'PRF1' , 'NKG7') 

library(stringr)  

pT <- DotPlot(sce.all, features = genes_to_check,

              assay='RNA'  )  + coord_flip()+ggtitle("Tcells")

pT

ggsave(plot=pT, filename="check_Tcells_marker_by_seurat_cluster.pdf")
### mast cells, TPSAB1 and TPSB2 

### B cell,  CD79A  and MS4A1 (CD20) 

### naive B cells, such as MS4A1 (CD20), CD19, CD22, TCL1A, and CD83, 

### plasma B cells, such as CD38, TNFRSF17 (BCMA), and IGHG1/IGHG4
### 2、Bcells

genes_to_check = c('CD3D','MS4A1','CD79A',

                   'CD19', 'CD22', 'TCL1A',  'CD83', ###  naive B cells

                   'CD38','TNFRSF17','IGHG1','IGHG4', ### plasma B cells,

                   'TPSAB1' , 'TPSB2',  ### mast cells,

                   'PTPRC' ) 

pB <- DotPlot(sce.all, features = genes_to_check,

              assay='RNA'  )  + coord_flip()+ggtitle("Bcells")

pB

ggsave(plot=pB, filename="check_Bcells_marker_by_seurat_cluster.pdf")
### 3、Myeloid

genes_to_check = c('CD68', 'CD163', 'CD14',  'CD86', 'LAMP3', #### DC 

                   'CD68',  'CD163','MRC1','MSR1','ITGAE','ITGAM','ITGAX','SIGLEC7', 

                   'MAF','APOE','FOLR2','RELB','BST2','BATF3')

pMyeloi <- DotPlot(sce.all, features = unique(genes_to_check),

                   assay='RNA'  )  + coord_flip()+ggtitle("Myeloid")

pMyeloi

ggsave(plot=pMyeloi, filename="check_Myeloid_marker_by_seurat_cluster.pdf")
### epi or tumor (EPCAM, KRT19, PROM1, ALDH1A1, CD24).

### - alveolar type I cell (AT1; AGER+)

### - alveolar type II cell (AT2; SFTPA1)

### - secretory club cell (Club; SCGB1A1+)

### - basal airway epithelial cells (Basal; KRT17+)

### - ciliated airway epithelial cells (Ciliated; TPPP3+)
### 4、Epi

genes_to_check = c(  'EPCAM' , 'KRT19', 'PROM1', 'ALDH1A1' ,

                     'AGER','SFTPA1','SCGB1A1','KRT17','TPPP3',

                     'KRT4','KRT14','KRT8','KRT18',

                     'CD3D','PTPRC' ) 

pEpi <- DotPlot(sce.all, features = unique(genes_to_check),

                assay='RNA'  )  + coord_flip()+ggtitle("Epi")

pEpi

ggsave(plot=pEpi, filename="check_epi_marker_by_seurat_cluster.pdf")
### 5、stromal

genes_to_check = c('TEK',"PTPRC","EPCAM","PDPN","PECAM1",'PDGFRB',

                   'CSPG4','GJB2', 'RGS5','ITGA7',

                   'ACTA2','RBP1','CD36', 'ADGRE5','COL11A1','FGF7', 'MME')

pStromal <- DotPlot(sce.all, features = unique(genes_to_check),

                    assay='RNA'  )  + coord_flip()+ggtitle("stromal")
pStromal

ggsave(plot=pStromal, filename="check_stromal_marker_by_seurat_cluster.pdf")

p5=pT+pB+pMyeloi+pEpi+pStromal ggsave(plot=p5, filename="pT+pB+pMyeloi+pEpi+pStromal.png",width=16 ,height=8) ggsave(plot=p5, filename="pT+pB+pMyeloi+pEpi+pStromal.pdf",width=16 ,height=8)

结合对应的关系，初步确认细胞类型

### 需要自行看图，定细胞亚群： celltype=data.frame(ClusterID=0:15, celltype='unkown') celltype[celltype$ClusterID %in% c(0,3,7),2]='Tcells' celltype[celltype$ClusterID %in% c(15),2]='NK' celltype[celltype$ClusterID %in% c(13),2]='Bcells' celltype[celltype$ClusterID %in% c(11),2]='myeloid' celltype[celltype$ClusterID %in% c(2,5,6,14),2]='epithelial' celltype[celltype$ClusterID %in% c( 1,8),2]='endo' celltype[celltype$ClusterID %in% c(4,9),2]='SMC' celltype[celltype$ClusterID %in% c(12),2]='fibo' celltype[celltype$ClusterID %in% c(10),2]='mast'


head(celltype)

celltype 

table(celltype$celltype)#9种细胞类型，基本和文章中一致
sce.all@meta.data$celltype = "NA"

for(i in 1:nrow(celltype)){

  sce.all@meta.data[which(sce.all@meta.data$CCA_snn_res.0.8 == celltype$ClusterID[i]),'celltype'] <- celltype$celltype[i]}

table(sce.all@meta.data$celltype)
genes_to_check = c('PTPRC', 'CD3D', 'CD3E', 'CD4','CD8A','NKG7',

                   'CD19', 'CD79A', 'MS4A1' ,

                   'IGHG1', 'MZB1', 'SDC1',

                   'CD68', 'CD163', 'CD14', 

                   'TPSAB1' , 'TPSB2',  ### mast cells,

                   'RCVRN','FPR1' , 'ITGAM' ,

                   'FGF7','MME', 'ACTA2',

                   'PECAM1', 'VWF', 

                   'EPCAM' , 'KRT19', 'PROM1', 'ALDH1A1' )

p <- DotPlot(sce.all, features = genes_to_check,

             assay='RNA' ,group.by = 'celltype' )  + coord_flip()

p ggsave(plot=p, filename="check_marker_by_celltype.pdf",width = 8, height = 5) ggsave(plot=p, filename="check_marker_by_celltype.png",width = 8, height = 5)

确定好细胞类型后，再利用这几种细胞类型的marker gene进行检查。

确定细胞类型，并可视化

#根据确定的细胞类型，给细胞加上注释，并对每个cluster进行注释 table(sce.all@meta.data$celltype,sce.all@meta.data$CCA_snn_res.0.8)

DimPlot(sce.all, reduction = "umap", group.by = "celltype",label = T) ggsave('umap_by_celltype.pdf') ggsave('umap_by_celltype.png') saveRDS(sce.all, "sce.all_cell.rds")

step6:细胞类型注释(基于singleR)

######## step6:细胞类型注释(基于singleR) ######## rm(list=ls()) library(Seurat) library(ggplot2) library(clustree) library(cowplot) library(dplyr)


getwd() 

#setwd("3-cell") 

sce=readRDS( "sce.all_cell.rds")
### sce=sce.all

sce

colnames(sce@meta.data)
library(SingleR) 

sce_for_SingleR <- GetAssayData(sce, slot="data") 

 

#使用HumanPrimaryCellAtlasData()函数从Human Primary Cell Atlas中获取参考数据，该函数返回一个SummarizedExperiment对象

### library(celldex)

### hpca.se <- HumanPrimaryCellAtlasData()

#load(file='../hpca.se.Rdata')

hpca.se

sel.clust = "CCA_snn_res.0.8"

sce <- SetIdent(sce, value = sel.clust)

table(sce@active.ident) 

table(sce@meta.data$CCA_snn_res.0.8)

clusters=as.character(sce@meta.data$CCA_snn_res.0.8)

table(clusters)  

colnames(sce@meta.data)

DimPlot(sce, reduction = "umap", group.by = "celltype",label = T) 

DimPlot(sce, reduction = "umap", group.by = "CCA_snn_res.0.8",label = T)
pred.hesc <- SingleR(test = sce_for_SingleR, ref = hpca.se, 

                     labels = hpca.se$label.main,

                     method = "cluster", clusters = clusters, 

                     assay.type.test = "logcounts", assay.type.ref = "logcounts")

table(pred.hesc$labels)

celltype = data.frame(ClusterID=rownames(pred.hesc), 

                      celltype=pred.hesc$labels, 

                      stringsAsFactors = F)

celltype

sce@meta.data$singleR=celltype[match(as.character(clusters),

                                     celltype$ClusterID),'celltype']
table(sce@meta.data$singleR,sce@meta.data$celltype)

colnames(sce@meta.data)
#CCA_snn_res.0.8  Tsne

p.dim.cell=DimPlot(sce, reduction = "tsne", group.by = "CCA_snn_res.0.8",label = T) 

p.dim.cell

ggsave(plot=p.dim.cell,filename="CCA_snn_res.0.8_tsne_cell.pdf",width=9, height=7)

ggsave(plot=p.dim.cell,filename="CCA_snn_res.0.8_tsne_cell.png",width=9, height=7)

#CCA_snn_res.0.8 Umap

p.dim.cell=DimPlot(sce, reduction = "umap", group.by = "CCA_snn_res.0.8",label = T) 

p.dim.cell

ggsave(plot=p.dim.cell,filename="CCA_snn_res.0.8_umap_cell.pdf",width=9, height=7)

ggsave(plot=p.dim.cell,filename="CCA_snn_res.0.8_umap_cell.png",width=9, height=7)

#singleR Tsne

p.dim.cell=DimPlot(sce, reduction = "tsne", group.by = "singleR",label = T) 

p.dim.cell

ggsave(plot=p.dim.cell,filename="singleR_tsne_cell.pdf",width=9, height=7)

ggsave(plot=p.dim.cell,filename="singleR_tsne_cell.png",width=9, height=7)

#singleR Umap

p.dim.cell=DimPlot(sce, reduction = "umap", group.by = "singleR",label = T) 

p.dim.cell

ggsave(plot=p.dim.cell,filename="singleR_umap_cell.pdf",width=9, height=7)

ggsave(plot=p.dim.cell,filename="singleR_umap_cell.png",width=9, height=7)

第一种方法需要人工收集marker基因比对各个cluster的显著高表达基因综合分析，第二种方法可以使用SingleR包自动识别细胞类型。

通过人工注释（celltype）和singleR包自动注释（singleR）的结果来看，鉴定结果有一定出入，但是差距不算太大，因此：建议两种方法结合起来进行细胞鉴定。

分别对pca1,pca2进行注释

#分别对pca1,pca2进行注释 #TSNE p.dim.cell=DimPlot(sce, reduction = "tsne", group.by = "singleR",split.by = "orig.ident", label = T) p.dim.cell ggsave(plot=p.dim.cell,filename="singleR_Dim_tsne_cell_split.pdf", width=13, height=7) ggsave(plot=p.dim.cell,filename="singleR_Dim_tsne_cell_split.png", width=13, height=7) #UMAP p.dim.cell=DimPlot(sce, reduction = "umap", group.by = "singleR",split.by = "orig.ident",label = T) p.dim.cell ggsave(plot=p.dim.cell,filename="singleR_Dim_umap_cell_split.pdf", width=13, height=7) ggsave(plot=p.dim.cell,filename="singleR_Dim_umap_cell_split.png", width=13, height=7)

分别对pca1和pca2进行聚类，注释，结果差别不大。

step7: 最终细胞类型的差异基因

######## step7: 最终细胞类型的差异基因 ######## library(gplots) tab.1=table(sce@meta.data$CCA_snn_res.0.8,sce@meta.data$celltype) tab.1 pro='cluster'#每个cluster的细胞类型以及数目 pdf(file = paste0(pro,'_last_celltype.pdf'),width = 12) balloonplot(tab.1, main =" last_celltype VS CCA_snn_res.0.8 ", xlab ="", ylab="", label = T, show.margins = F) dev.off()

#保存meta.data phe=sce@meta.data save(phe,file = 'last_phe.Rdata')

1、可视化每个cluster有多少个细胞，以及他们属于什么细胞类型。

2、可视化每种细胞类型的基因数和Count数

#可视化每种细胞类型的基因数和Count数 feats <- c("nFeature_RNA", "nCount_RNA") p1_filtered=VlnPlot(sce, group.by = "celltype", features = feats, pt.size = 0.1, ncol = 2) + NoLegend() p1_filtered ggsave(filename="Vlnplot1_filtered.pdf",plot=p1_filtered) ggsave(filename="Vlnplot1_filtered.png",plot=p1_filtered)

3、分别可视化"percent_mito", "percent_ribo", "percent_hb"在不同类型细胞中的比例

#分别可视化"percent_mito", "percent_ribo", "percent_hb"在不同类型细胞中的比例 feats <- c("percent_mito", "percent_ribo", "percent_hb") p2_filtered=VlnPlot(sce, group.by = "celltype", features = feats, pt.size = 0.1, ncol = 3) + NoLegend() p2_filtered ggsave(filename="Vlnplot2_filtered.pdf",plot=p2_filtered) ggsave(filename="Vlnplot2_filtered.png",plot=p2_filtered)

step8: 看不同细胞亚群比例变化

######## step8: 看不同细胞亚群比例变化 ######## library(ggplot2) #setwd("3-cell") ### sce.all=readRDS( "sce.all_int.rds") load(file = 'last_phe.Rdata') head(phe)

可视化所有细胞注释分类的饼图、柱状图、玫瑰花图

#可视化2：所有细胞注释分类的饼图 tb=as.data.frame(table(phe$celltype),stringsAsFactors = F) tb=tb[order(tb$Freq,decreasing=F),] colnames(tb) <- c('Celltype', 'Freq') percentage <- scales::percent(tb$Freq / sum(tb$Freq)) percentage labs <- paste0(tb$Celltype,'(', percentage, ')')#设置标签名 labs tb$cell=factor(tb$Celltype,levels=tb$Celltype) p.pie = ggplot(tb, aes(x = "", y = Freq, fill = cell)) + geom_bar(stat = "identity",width = 1) + scale_fill_discrete(breaks = tb$cell, labels = labs) + labs(x = "", y = "", title = "") + coord_polar(theta = "y") + theme(axis.ticks = element_blank()) p.pie ggsave(plot=p.pie,filename="Total_pie.pdf",width=14, height=7) ggsave(plot=p.pie,filename="Total_pie.png",width=14, height=7) #柱状图 head(tb) tb$percentage <-tb$Freq / sum(tb$Freq) tb$Celltype = factor(tb$Celltype,tb$Celltype) p <- ggplot(tb, aes(x=Celltype, y=percentage, fill=Celltype)) + geom_bar(stat="identity") +theme_bw() + theme(legend.position ="none") + geom_text( aes(x=Celltype, y=percentage, label=labs ), #check_overlap = T, size=4, vjust = -1,fontface="bold") p ggsave(plot=p,filename="Total_pie-barplot.pdf",width=14, height=7) ggsave(plot=p,filename="Total_pie-barplot.png",width=14, height=7)


#绘制玫瑰花图：

p2 <- p+coord_polar()+labs(x = "", y = "", title = "Cell Types") + 

  theme(axis.text.y = element_blank()) +     

  theme(axis.ticks = element_blank()) +     

  theme(panel.border = element_blank()) + 

  theme(plot.title = element_text(hjust=0.5,size=14,face = "bold") )+

  theme(axis.text.x=  element_blank() )

p2

ggsave(plot=p2,filename="Total_pie2.pdf",width=14, height=7)

ggsave(plot=p2,filename="Total_pie2.png",width=14, height=7)

这三种方式画出来的图太漂亮了吧！

可视化每个sample/group的细胞类型组成柱状图

#可视化3：每个sample/group的细胞类型组成柱状图 colnames(phe) table(phe$orig.ident) #test=phe[,c("orig.ident" , "patient","group","celltype")] test=phe[,c("orig.ident" ,"celltype")] test$cell=factor(test$celltype,levels=rev(tb$Celltype)) head(test) p.bar.sample=ggplot(test, aes(x=orig.ident, fill=cell)) + geom_bar(position="fill") + labs(title='Cell proportion of each sample') +theme_bw() p.bar.sample ggsave(plot=p.bar.sample,filename="Bar.sample.pdf",width=14, height=7) ggsave(plot=p.bar.sample,filename="Bar.sample.png",width=14, height=7) ###

总结：

用规范的数据来复现文章真的很有必要，既能加深对代码的理解，又能拓宽自己的思路，并且利用比较规范的数据来练手不会那么容易踩坑。

加油鸭！

写在文末

咱们现在这个专栏《100个单细胞转录组数据降维聚类分群图表复现》分享的代码是到此为止，但是一般来说单细胞文章数据分析还有很多进阶图表制作，比如inferCNV看肿瘤拷贝数变异，monocle看拟时序等等。如果你也需要，可以加入我们这个专栏《100个单细胞转录组数据降维聚类分群图表复现》创作团队，获取进阶指引哦！

文末友情推荐

与十万人一起学生信，你值得拥有下面的学习班：

小样本量单细胞肿瘤数据分析模范