​harmony、不harmony，这是个问题

harmony基本原理

harmony应用主成分分析，将转录组表达谱嵌入到低维空间中，然后应用迭代过程去除数据集特有的影响

我在跑单细胞前单独学习harmony算法时就有这样的感悟：

感觉这种方法默认需要批次效应不是很大，大到让同一类型细胞PCA却分的很开

使用harmony等算法去除单细胞样品差异得谨慎

这篇推文评论区也有类似的声音

使用harmony算法与否对GZMA在细胞分群中表达量无影响

原推文已有代码略

这周的推文是对一篇文献 《Systemic Blockade of Clever-1 Elicits Lymphocyte Activation Alongside Checkpoint Molecule Downregulation in Patients with Solid Tumors: Results from a Phase I/II Clinical Trial》中的数据集进行复现，该篇文献只对Tcell进行单细胞测序。整合C1D0 (C1D1),C4D7 (C4D8),C4D14 (C4D15)样本数据进行第一层次降维分群，之后再对CD8+ Tcell进行二次降维分群，我们拿到的数据是 CD8A-expressing CD8+ T-cell的csv文件。

所以这篇推文就是要看看harmony与否对作者文章中的结果有无影响

样本名称:

导入数据

标准化和降维分群

特征需要跟生物背景一致

这里很明显看到C1D1与其他两种样本存在差异

harmony前：

可以发现GZMA在cluster3很明显表达量很低，与作者的图一致

目的：

根据前面的图，我们把GZMA人为分为GZMA+和GZMA-来比对一下harmony前后对该基因在细胞中的表达量是否会不同

原文部分代码有误

counts2 <- sce
counts2@meta.data$celltype = "NA"
table(counts2@meta.data$RNA_snn_res.0.6 )

for(i in 1:nrow(celltype)){
  counts2@meta.data[which(counts2@meta.data$RNA_snn_res.0.6 == celltype$ClusterID[i]),'celltype'] <- celltype$celltype[i]}
table(counts2@meta.data$celltype)

#harmony前
tab.1=table(counts2@meta.data$celltype,counts2@meta.data$seurat_clusters)
colnames(counts2@meta.data)   
library(gplots)
tab.1 
pro='before_harmony'
pdf(file = paste0(pro,'before_GZMA.pdf'),
    width = 9,height  = 8)
balloonplot(tab.1, main ="GZMA- VS  GZMA+ ", xlab ="", ylab="",
            label = T, show.margins = F)
dev.off() 

harmony

seuratObj <- RunHarmony(sce, "orig.ident")
names(seuratObj@reductions)
seuratObj <- RunUMAP(seuratObj,  dims = 1:15, 
                     reduction = "harmony")
DimPlot(seuratObj,reduction = "umap",label=T )  

sce=seuratObj
sce <- FindNeighbors(sce, reduction = "harmony",
                     dims = 1:15)
sce <- FindClusters(sce, resolution = 0.6)
table(sce@meta.data$seurat_clusters)
DimPlot(sce,reduction = "umap",label=T)   

+这里harmony后代码所画的两幅Dimplot有什么区别？

（harmony前）

sce <- FindVariableFeatures(sce, 
                            selection.method = "vst", nfeatures = 2000)  
sce <- ScaleData(sce) 
sce <- RunPCA(object = sce, pc.genes = VariableFeatures(sce)) # 默认输出5个PC
DimHeatmap(sce, dims = 1:12, cells = 100, balanced = TRUE)  # 画12个PC
ElbowPlot(sce) 

sce <- FindNeighbors(sce, dims = 1:15)  # 输入15个维度
sce <- FindClusters(sce, resolution = 0.6)  # 作者0.6
table(sce@meta.data$RNA_snn_res.0.6)  

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

（harmony后）

seuratObj <- RunUMAP(seuratObj,  dims = 1:15, 
                     reduction = "harmony")
DimPlot(seuratObj,reduction = "umap",label=T )  

sce=seuratObj
sce <- FindNeighbors(sce, reduction = "harmony",
                     dims = 1:15)
sce <- FindClusters(sce, resolution = 0.6)
table(sce@meta.data$seurat_clusters)
DimPlot(sce,reduction = "umap",label=T) 

我们上期也提到，harmony前后相同分辨率下UMAP降维分组会发生变化

这两幅Dimplot中

第一张，区别于harmony前的UMAP图，reduction参数被设置为了"harmony"

绘制UMAP前要有基本降维

绘制UMAP的一般流程：

初探单细胞下游

• 数据标准化、归一化
• 识别高变基因
• PCA降维，默认使用前面找到的高变基因的scale矩阵用于降维
• UMAP聚类
• • FindNeighbors函数通过计算数据集中每个细胞之间的相似性，找到每个细胞的k个最近邻，并可选地计算共享最近邻图。这些最近邻关系可以用来构建细胞之间的连接，用于后续的聚类分析、可视化和其他细胞间关系的研究
• • 通过调用 FindClusters函数，可以根据细胞之间的共享最近邻关系，在数据集中识别出具有相似性的细胞聚类

所以第一张Dimplot属于，在harmony将PCA降维好的细胞分布迭代优化后的聚类分布图（无cluster id）上，UMAP使用的仍是harmony前计算的cluster id的label

体现出来就是分布图形形状变化

而第二张Dimplot使用 FindNeighbors函数和 FindClusters函数重新计算并给细胞聚类打上cluster id label

这也导致了seurat对象中seurat_clusters水平发生变化

使用的是harmony优化后的聚类分布图重新计算后的label

这里其实可以看到原推文harmony前后降维分组 GZMA-cluster 所含细胞数量 vs GZMA+cluster 所含细胞数量图 放反了

两次数据分析对比结果可视化

#我们将两次数据未进行harmony和 进行harmony的对比
tab.2=table(counts2@meta.data$celltype,sce@meta.data$celltype)
balloonplot(tab.2, main =" before harmony VS after harmony", xlab ="", ylab="",
            label = T, show.margins = F)

最后通过table得到的balloonplot发现其实做不做harmony对基因GZMA在细胞分群中表达量基本上没有影响。

其实我个人看来，从整个scRNAseq分析来说，上述只能说明对这一步”基因GZMA在细胞分群中表达量“没有影响

harmony究竟该不该做，还需下面进一步讨论

两个问题

harmony后没有：

①单细胞中什么时候样本间差异需要去除？

根据生物背景知识

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

第一次降维分组粗略合并命名... 但是这样的分析太常规了，远不够！一般来说， 需要进行各自细胞亚群独立进行更加细致的降维聚类分群... 这里我们拿fibo亚群举例：

但是常见fibo相关的基因都不特征性表达，所以很可能是样本差异 强烈建议这个时候使用harmony进行整合

抹去了样本差异

使用harmony等算法去除单细胞样品差异得谨慎

15个鼻咽癌样品，加上1个正常人样品。全部的样品的单细胞转录组数据整合后，如果不使用harmony等算法去除样品差异，默认的降维聚类分群，如下所示：

可以看到，效果还不错，很有意思， 给大家的感觉是 harmony等算法去除样品差异并不是必须的。 但是如果我们具体到每个样品，有如下所示的现象：

上皮细胞大的亚群里面，每个病人独立成为小亚群，泾渭分明，这个符合预期，因为每个肿瘤病人都有自己的特异性 但是免疫细胞各个亚群里面，病人之间的界限就模糊很多。值得注意的是P07这个病人的样品，它主要是T细胞和髓系细胞，而且是独立成为一个亚群了，这就是单细胞转录组的样品差异，理论上是需要去除的！

fibo、肿瘤上皮细胞、免疫细胞...

所以我们需要根据不同细胞类型和样本降维聚类分组图来判断需不需要去除批次效应

②harmony后可能导致的问题

使用harmony等算法去除单细胞样品差异得谨慎

... 可以看到，首先上皮细胞大的亚群里面，每个病人独立成为小亚群，泾渭分明，这个符合预期，因为每个肿瘤病人都有自己的特异性。但是免疫细胞各个亚群里面，病人之间的界限就模糊很多。值得注意的是P07这个病人的样品，它主要是T细胞和髓系细胞，而且是独立成为一个亚群了，这就是单细胞转录组的样品差异，理论上是需要去除的！ 有意思的事情就来了

第一层次降维聚类分群最好是分辨率调大一点

我也检测了这些上皮细胞看看它们是否具有病人异质性，因为我前面使用了harmony进行多个样品整合，所以才能造成多个病人的正常细胞（成纤维，内皮细胞，免疫细胞）能很好的混合在一起，但是副作用就是多个病人的本来就应该是千差万别的肿瘤上皮细胞也被迫糅合在一起啦。

根据这两个问题我们整理的资料

探索本文数据集harmony前后对已定义的细胞亚型聚类的影响：

这里进行降维分群需要注意的是，因为实验只是对Tcell测序，所以这里第一次定义细胞亚群其实已经需要精细的marker

如果不经过思考当前研究的背景，直接拿我们前面的代码走，可以发现：

实验只是对Tcell测序可以看到Tcell的相关marker表达量各个cluster中都较高

免疫治疗前后T细胞亚群比例变化真的重要吗？准确吗？ 'PTPRC', 'CD3D', 'CD3E', 'CD4','CD8A'

> genes_check_after_harmony <- DotPlot(sce, features = c("EPCAM"),
+                                      assay='RNA'  )  + coord_flip()
Error in FetchData.Seurat(object = object, vars = features, cells = cells) : 
  None of the requested variables were found: EPCAM

所以虽然是第一次降维分组，实际上已经需要精细定义

T细胞 CD4和CD8的T细胞在单细胞转录组水平本来就很难确定亚群和名字 各个单细胞亚群的top特异性高表达基因也不会跟第一层次那样的泾渭分明。通常我们拿到了肿瘤相关的单细胞转录组的表达量矩阵后的第一层次降维聚类分群通常是：

• immune (CD45+,PTPRC),
• epithelial/cancer (EpCAM+,EPCAM),
• stromal (CD10+,MME,fibro or CD31+,PECAM1,endo)

参考我前面介绍过 CNS图表复现08—肿瘤单细胞数据第一次分群通用规则，这3大单细胞亚群构成了肿瘤免疫微环境的复杂。绝大部分文章都是抓住免疫细胞亚群进行细分，包括淋巴系（T,B,NK细胞）和髓系（单核，树突，巨噬，粒细胞）的两大类作为第二次细分亚群。但是也有不少文章是抓住stromal 里面的 fibro 和endo进行细分，并且编造生物学故事的。 分享一些用得着的CD4和CD8的T细胞的基因吧 比如按照功能进行划分，naive, memory ，effector，cytotoxic，Exhaustion：

• naive (LEF1, SELL, TCF7),
• effector (IFNG),
• cytotoxicity (GZMB, PRF1),
• early and general exhaustion (PDCD1, CTLA4, ENTPD1 ) .
• antigen presentation (CD74, HLA-DRB1/5, HLA-DQA2)

如果你使用上面的基因列表，你会发现主要的naive状态的是CD4的T细胞，其它主要是CD8的T细胞。 还有：

• naive T cell, such as CCR7, LTB, and SELL
• memory CD4+ T cells (TEM) CXCR4, TNFRSF4, and CCR6, RORA, IL6ST, and IL17RA, suggesting a TH17 cell identity
• cytotoxicity associated genes (GZMH, GZMA, and GZMB) were expressed in cluster T4
• T cell receptor alpha and beta chain variable genes with TRAV13-2 and TRBV7-9 were enriched in cluster T5
• Cluster T6 was characterized as naive CD8+ T cells with the high expression of CCR7, LEF1, and CD27.

（harmony前）

c("LEF1", "Sell", "TCF7", "IFNG", "GZMB", "PRF1", "PDCD1", "CTLA4", "ENTPD1", "HLA-DRB1/5", "HLA", "HLA-DQA2")

> DotPlot(counts2, 
+         features = genes_to_check1,
+         assay='RNA'  )  + coord_flip()
Warning message:
In FetchData.Seurat(object = object, vars = features, cells = cells) :
  The following requested variables were not found: Sell, ENTPD1, HLA-DRB1/5, HLA, HLA-DQA2

c("CCR7", "LTB", "CXCR4", "TNFRSF4", "CCR6", "RORA", "IL6ST", "IL17RA", "GZMH", "GZMA", "GZMB", "TRAV13-2", "TRBV7-9", "CCR7", "LEF1", "CD27")

很明显cluster3为naive T cell

其它其实分的也不是很开

和我们前面的GZMB+ vs GZMB- 分组是一样的

更多的时候其实是打分，并不能完全是看某个基因或者某些基因在某个单细胞亚群里面的排他性的特异性高表达

（harmony后）

如果从目前的简单分组来看harmony没有把定义好的细胞亚型打散，这也在意料之中，因为我们这一步已经全部是T细胞了，不会像前面提到的那些推文中那样epi、fibo、immu...

但我们这里还是整理一下

针对harmony打散已定义细胞亚型情况解决方法

就像本文开头提到我一开始学的时候预感的那样

感觉这种方法默认需要批次效应不是很大，大到让同一类型细胞分的很开

• 当检查各批次之间降维聚类分组无明显区分时，在进行第一次降维分组（这里默认只进行两次分组的话，第一次粗略，第二次精细）时我们并不进行harmony，因为这时可能会去掉细胞的异质性

此外在第一次降维分组时还需注意：

第一层次降维聚类分群最好是分辨率调大一点

第一层次降维聚类分群最好是分辨率调大一点，然后根据我们给大家的基因列表背后的生物学意义去给它们亚群进行生物学命名。然后针对不同的单细胞亚群进行细分即可，这样的单细胞转录组数据分析的标准降维聚类分群，并且进行生物学注释后的结果。

第一次初步定义细胞亚群往往用的是通用规则：

听说你还缺单细胞亚群标记基因

免疫治疗前后T细胞亚群比例变化真的重要吗？准确吗？

# 一般来说肿瘤样品的单细胞就应该是首先是按照如下所示的标记基因进行第一次分群 ：
# 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) 
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' ,
                   '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' )

• 在进行精细分组进一步定义亚型时，选取对应子集，才进行harmony，因为这时批次效应不会那么大，而且不会影响到粗略分组定义好的其它亚型

就我这个想法问了曾老师 曾老师说根据经验，一般是有肿瘤上皮细胞第一步不harmony，其它基本上都可以 我的思考：确实，肿瘤细胞病人异质性很大，会造成较大的批次效应

针对本文这个数据集，虽然只进行一次降维分组，但因为已经是只对Tcell测序，所以已经不是初步粗略分组了，鉴于样本间的批次效应，感觉是要harmony为好。

往期回顾

NK细胞的单细胞层面细分亚群

被忽视的真相 | 细胞亚群内互作

单细胞分辨率下鉴定中胚层诱导的 ESC 中的转录组学、调控网络和增强子

跟学单细胞周更（二）

单细胞测序最好的教程（七）: 数据整合与批次效应校正