细分亚群后需要使用harmony去除样品差异

生信技能树

发布于 2021-08-25 17:21:17

1.4K0

发布于 2021-08-25 17:21:17

文章被收录于专栏：生信技能树

经过了大量的单细胞转录组数据分析基础讲解，相信大家对第一层次降维聚类分群都不陌生了。参考我们的《明码标价》专栏里面的单细胞内容

如下所示的一个数据集处理，我们会得到如下所示：

基础降维聚类分群

特定的基因名字对应的特定的生物学亚群，如下所示：

'CD68', 'CD163', 'CD14', # myeloids
'TPSAB1' , 'TPSB2',  # mast cells,
'RCVRN','FPR1' , 'ITGAM' ,
'C1QA',  'C1QB',  # mac
'S100A9', 'S100A8', 'MMP19',# monocyte
'LAMP3', 'IDO1','IDO2',## DC3 
'CD1E','CD1C', # DC2
'KLRB1','NCR1', # NK 
'FGF7','MME', 'ACTA2', ## fibo 
'DCN', 'LUM',  'GSN' , ## mouse PDAC fibo 
'Amy1' , 'Amy2a2', # Acinar_cells
'PECAM1', 'VWF',  ## endo 
'EPCAM' , 'KRT19', 'PROM1', 'ALDH1A1'  # epi

简单的肉眼就可以粗略进行各自细胞亚群合并后的生物学命名，如下所示：

生物学命名

但是这样的分析太常规了，远不够！一般来说，需要进行各自细胞亚群独立进行更加细致的降维聚类分群：

#拆分为 多个 seurat子对象
sce.all.list <- SplitObject(sce.all , split.by = "celltype")
sce.all.list 
names(sce.all.list)
for (i in names(sce.all.list)) {
  epi_mat = sce.all.list[[i]]@assays$RNA@counts
  epi_phe = sce.all.list[[i]]@meta.data
  sce=CreateSeuratObject(counts = epi_mat, 
                         meta.data = epi_phe )
  sce
  table(sce@meta.data$orig.ident) 
  save(sce,file = paste0(i,'.Rdata'))
}

得到如下所示各个亚群的rdata文件，都存放在 Rdata-celltype/ 文件夹里面：

$ ls -lh Rdata-celltype/|cut -d" " -f 5-

5.6M Aug 22 09:39 Bcells.Rdata
 14M Aug 22 09:39 CD4.Rdata
 22M Aug 22 09:39 CD8.Rdata
3.2M Aug 22 09:40 NK.Rdata
1.3M Aug 22 09:39 PTGDS.Rdata
2.4M Aug 22 09:40 SMC.Rdata
 26K Aug 22 09:39 check_stromal_marker_by_RNA_snn_res.0.8.pdf
5.8M Aug 22 09:40 endo.Rdata
 63M Aug 22 09:39 epi.Rdata
 13M Aug 22 09:39 fibo.Rdata
 41M Aug 22 09:40 myeloid.Rdata
6.8M Aug 22 09:39 plasma.Rdata
7.2M Aug 22 09:39 treg.Rdata

继续降维聚类分群

这里我们拿fibo亚群举例：

rm(list=ls())
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)

load('../Rdata-celltype/fibo.Rdata')
table(sce@meta.data$orig.ident) 
sce  
library(stringr) 
sce <- NormalizeData(sce, normalization.method =  "LogNormalize",  
                     scale.factor = 1e4)
GetAssay(sce,assay = "RNA")
sce <- FindVariableFeatures(sce, 
                            selection.method = "vst", nfeatures = 2000)  
sce <- ScaleData(sce) 
sce <- RunPCA(object = sce, pc.genes = VariableFeatures(sce)) 
DimHeatmap(sce, dims = 1:12, cells = 100, balanced = TRUE)
ElbowPlot(sce) 

sce <- FindNeighbors(sce, dims = 1:15)
sce <- FindClusters(sce, resolution = 0.8)
table(sce@meta.data$RNA_snn_res.0.8)  

set.seed(123)
sce <- RunTSNE(object = sce, dims = 1:15, do.fast = TRUE)
DimPlot(sce,reduction = "tsne",label=T)

sce <- RunUMAP(object = sce, dims = 1:15, do.fast = TRUE)
DimPlot(sce,reduction = "umap",label=T) 
head(sce@meta.data)
DimPlot(sce,reduction = "umap",label=T,group.by = 'orig.ident') 

sce <- FindClusters(sce, resolution = 0.1)
DimPlot(sce,reduction = "umap",label=T) 
ggsave(filename = 'umap_recluster_by_0.1.pdf') 
DimPlot(sce,reduction = "umap",label=T,
        group.by = 'orig.ident') 
ggsave(filename = 'umap_sce_recluster_by_patients.pdf')

table(sce@meta.data$orig.ident,sce@meta.data$seurat_clusters)
colnames(sce@meta.data)

可以看到，如下所示的非常尴尬的降维聚类分群结果：

每个亚群都是各自的病人特异

也就是说这个时候区分出来的每个亚群都是各自的病人特异的，有点类似于肿瘤细胞的特性。

如下所示：

每个亚群都是各自的病人特异umap

而且常见fibo相关的基因都不特征性表达：

亚群并不是fibo的戏份

这个时候根本就没办法对这些细胞亚群进行生物学命名。

强烈建议这个时候使用harmony进行整合

代码如下所示；


library(harmony)
seuratObj <- RunHarmony(sce, "orig.ident")
names(seuratObj@reductions)
seuratObj <- RunUMAP(seuratObj,  dims = 1:15, 
                     reduction = "harmony")
DimPlot(seuratObj,reduction = "umap",label=T ) 
p1.compare=plot_grid(ncol = 3,
                     DimPlot(sce, reduction = "pca", group.by = "orig.ident")+NoAxes()+ggtitle("PCA"),
                     DimPlot(sce, reduction = "tsne", group.by = "orig.ident")+NoAxes()+ggtitle("tSNE"),
                     DimPlot(sce, reduction = "umap", group.by = "orig.ident")+NoAxes()+ggtitle("UMAP"),
                     DimPlot(seuratObj, reduction = "pca", group.by = "orig.ident")+NoAxes()+ggtitle("PCA"),
                     DimPlot(seuratObj, reduction = "tsne", group.by = "orig.ident")+NoAxes()+ggtitle("tSNE"),
                     DimPlot(seuratObj, reduction = "umap", group.by = "orig.ident")+NoAxes()+ggtitle("UMAP")
)
p1.compare
ggsave(plot=p1.compare,filename="before_and_after_inter.pdf")
 
sce=seuratObj
sce <- FindNeighbors(sce, reduction = "harmony",dims = 1:15)
sce <- FindClusters(sce, resolution = 0.1)
table(sce@meta.data$seurat_clusters)
DimPlot(sce,reduction = "umap",label=T) 
ggsave(filename = 'harmony_umap_recluster_by_0.1.pdf') 
DimPlot(sce,reduction = "umap",label=T,
        group.by = 'orig.ident') 
ggsave(filename = 'harmony_umap_sce_recluster_by_patients.pdf')

整合前后效果很明显：