hello，下半年開始，我們也要開始新的工作了，今天我們需要解決的問題有兩個(gè)。

（1）細(xì)胞亞群細(xì)分的時(shí)候仍然要分開樣本去除批次效應(yīng)？？？harmony

（2）harmony算法是不是剛才是整合分析就用，還是說再分群再用？？？

其實(shí)就是涉及到harmony整合分析的問題。

首先第一個(gè)問題，再分群批次去除的問題。

我們做腫瘤研究的單細(xì)胞數(shù)據(jù)，一般來說會選擇初步很粗狂的定義大的細(xì)胞亞群，比如我常用的 第一次分群是通用規(guī)則是：

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

然后絕大部分文章都是抓住免疫細(xì)胞亞群進(jìn)行細(xì)分，包括淋巴系（T,B,NK細(xì)胞）和髓系（單核，樹突，巨噬，粒細(xì)胞）的兩大類作為第二次細(xì)分亞群。說起來很簡單，但是實(shí)際上每次做到單細(xì)胞數(shù)據(jù)集的細(xì)分亞群就非常的頭疼，尤其是myeloid的髓系，（單核，樹突，巨噬，粒細(xì)胞）有時(shí)候根本就分不清楚，而且分完之后仍然是可以繼續(xù)細(xì)分。

我嘗試對它這個(gè)數(shù)據(jù)集進(jìn)行數(shù)據(jù)分析圖表復(fù)現(xiàn)，比如單獨(dú)把髓系拿出來進(jìn)行重新降維聚類分群，然后可視化如下所示：

load(file = 'sce_recluster.Rdata')
p1=DimPlot(sce,reduction = "umap",label=T
        ,group.by = 'Cell_type') 
p2=DimPlot(sce,reduction = "umap",label=T
        ,group.by = 'orig.ident') 
table(sce$orig.ident,sce$seurat_clusters)
p3=DimPlot(sce,reduction = "umap",label=T) 
library(patchwork)
p1/p2/p3

可以看到pDC這個(gè)細(xì)胞亞群，由 4和8群組成，而且包含多個(gè)病人！

但是p07這個(gè)病人就非常的詭異，這個(gè)樣品里面的多種髓系細(xì)胞居然是與其它病人樣品的髓系細(xì)胞距離超級遠(yuǎn)！

圖片

但是如果你harmony處理一下，然后再降維聚類分群，代碼如下所示：

load(file = 'main_sce_recluster.Rdata')

sce.all.filt=sce
library(harmony) 
sce.all.int <- RunHarmony(sce.all.filt,
                          c( "orig.ident" ))
names(sce.all.int@reductions)
harmony_embeddings <- Embeddings(sce.all.int, 'harmony')
harmony_embeddings[1:5, 1:5]

sce.all.int=RunTSNE(sce.all.int,reduction = "harmony", dims = 1:30)
sce.all.int=RunUMAP(sce.all.int,reduction = "harmony",dims = 1:10)

sce=sce.all.int
sce <- FindNeighbors(sce, reduction = "harmony",dims = 1:15)
sce <- FindClusters(sce, resolution = 0.8)

load(file = 'after_harmony/main_sce_recluster.Rdata')
p1=DimPlot(sce,reduction = "umap",label=T
           ,group.by = 'Cell_type') 
p2=DimPlot(sce,reduction = "umap",label=T
           ,group.by = 'orig.ident') 
table(sce$orig.ident,sce$seurat_clusters)
p3=DimPlot(sce,reduction = "umap",label=T) 
library(patchwork)
p1/p2/p3

出圖如下：

圖片

可以看到這個(gè)時(shí)候的p07這個(gè)病人不會在獨(dú)立成為一個(gè)亞群啦，而且呢，pDC這個(gè)細(xì)胞亞群也比較純粹一點(diǎn)了,說明再分群的時(shí)候，需要去除批次效應(yīng)。

第二個(gè)問題：harmony是不是一開始就用。

我們以這個(gè)單細(xì)胞轉(zhuǎn)錄組文獻(xiàn)，《Single-cell transcriptomics reveals regulators underlying immune cell diversity and immune subtypes associated with prognosis in nasopharyngeal carcinoma》為例子，15個(gè)鼻咽癌樣品，加上1個(gè)正常人樣品。全部的樣品的單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)整合后，如果不使用harmony等算法去除樣品差異，默認(rèn)的降維聚類分群，如下所示：

圖片

我們根據(jù)左邊的標(biāo)記基因以及生物學(xué)背景知識，可以進(jìn)行如下所示的命名：

圖片

<figcaption style="margin-top: 5px;text-align: center;color: #888;font-size: 14px;"> </figcaption>

可以看到，效果還不錯，很有意思， 給大家的感覺是 harmony等算法去除樣品差異并不是必須的。但是如果我們具體到每個(gè)樣品，有如下所示的現(xiàn)象：

圖片

可以看到，首先上皮細(xì)胞大的亞群里面，每個(gè)病人獨(dú)立成為小亞群，涇渭分明，這個(gè)符合預(yù)期，因?yàn)槊總€(gè)腫瘤病人都有自己的特異性。但是免疫細(xì)胞各個(gè)亞群里面，病人之間的界限就模糊很多。值得注意的是P07這個(gè)病人的樣品，它主要是T細(xì)胞和髓系細(xì)胞，而且是獨(dú)立成為一個(gè)亞群了，這就是單細(xì)胞轉(zhuǎn)錄組的樣品差異，理論上是需要去除的！

有意思的事情就來了

如果我們在樣品層面就開始使用harmony等算法去除樣品差異，又會導(dǎo)致另外一個(gè)可怕的事情發(fā)生，如下所示：

圖片

<figcaption style="margin-top: 5px;text-align: center;color: #888;font-size: 14px;"> </figcaption>

就是本來是應(yīng)該是具備病人特異性的上皮細(xì)胞，這個(gè)時(shí)候被抹除了樣品差異。

好好的上皮細(xì)胞，被拆分的七零八落，如下所示：

圖片.png

我們也可以以病人樣品視角來看：

圖片

這個(gè)算法真的是太可怕了，樣品差異被抹除的干干凈凈了！這不是最可怕的，真正的問題是，這個(gè)上皮細(xì)胞被打散到了其它免疫細(xì)胞里面，因?yàn)檫@個(gè)harmony算法！我們可以對上皮細(xì)胞的最重要的marker基因EPCAM進(jìn)行如下所示可視化，并且使用harmony等算法去除樣品差異前后可以對比看看。

harmony雖好，但也不要貪杯啊