其中第一层次降维聚类分群是：clustering and cell type identification was based on scRNA-seq of 212,404 cells from primary mononuclear bone marrow samples of RRMM patients (n = 20) as shown by an UMAP embedding colored by sample without batch-effect correction.
然后呢，效果如下所示，感觉就有点奇怪，因为恰好是肿瘤细胞和正常细胞在umap图上面就泾渭分明成为了两个部分（下图的上半部分是肿瘤细胞，然后下半部分就是正常细胞）：
肿瘤细胞和正常细胞在umap图上面就泾渭分明
起初我是不太相信的，然后我就下载了这个数据集的GSE161801_RAW.tar压缩包然后解压后如下所示的文本文件：
代码语言：javascript复制$ ls -lh *gz|cut -d" " -f 5-
248K Nov 18 2020 GSM4914649_K43R_2FTMUU_T1.csv.gz
3.3M Nov 18 2020 GSM4914650_K43R_2FTMUU_N1.csv.gz
4.7M Nov 18 2020 GSM4914651_K43R_2FTMUU_T2.csv.gz
5.3M Nov 18 2020 GSM4914652_K43R_2FTMUU_N2.csv.gz
2.6M Nov 18 2020 GSM4914653_K43R_32ARQL_T2.csv.gz
3.6M Nov 18 2020 GSM4914654_K43R_32ARQL_N2.csv.gz
7.4M Nov 18 2020 GSM4914655_K43R_37A66E_T1.csv.gz
9.2M Nov 18 2020 GSM4914656_K43R_37A66E_N1.csv.gz
3.6M Nov 18 2020 GSM4914657_K43R_37A66E_T2.csv.gz
1.6M Nov 18 2020 GSM4914658_K43R_37A66E_N2.csv.gz
1.1M Nov 18 2020 GSM4914659_K43R_6K8YTY_T2.csv.gz
。。。。。。。。。。。。。。。。。
所以就批量读取它们，代码如下所示：
代码语言：javascript复制library(data.table)
dir='../GSE161801_RAW/'
samples=list.files( dir ,pattern = 'gz')
samples
library(data.table)
sceList = lapply(samples,function(pro){
# pro=samples[1]
print(pro)
ct=fread(file.path( dir ,pro),data.table = F)
ct[1:4,1:4]
rownames(ct)=ct[,1]
ct=ct[,-1]
sce=CreateSeuratObject(counts = ct ,
project = gsub('matrix.txt.gz','',
strsplit(pro,'_')[[1]][3]) ,
min.cells = 5,
min.features = 300 )
table(sce$orig.ident)
sce$orig.ident = gsub('matrix.txt.gz','',
strsplit(pro,'_')[[1]][3])
return(sce)
})
names(sceList)
library(stringr)
samples=gsub('matrix.txt.gz','',str_split(samples,'_',simplify = T)[,3])
samples
names(sceList) = samples
sce.all=merge(x=sceList[[1]],
y=sceList[ -1 ],
add.cell.ids = samples)
as.data.frame(sce.all@assays$RNA@counts[1:10, 1:2])
head(sce.all@meta.data, 10)
table(sce.all$orig.ident)
一般情况下，大家会对上面的全部的样品的Seurat对象进行harmony整合，示例代码如下所示：
代码语言：javascript复制 sce=sce.all
sce <- NormalizeData(sce,
normalization.method = "LogNormalize",
scale.factor = 1e4)
sce <- FindVariableFeatures(sce)
sce <- ScaleData(sce)
sce <- RunPCA(sce, features = VariableFeatures(object = sce))
if(F){
library(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 <- RunUMAP(sce, dims = 1:15 )
sce <- FindNeighbors(sce, dims = 1:15)
sce.all=sce
FeaturePlot(sce.all,raster = F,features = 'TNFRSF17')
但是我为了还原文章的结果，选择了放弃harmony整合后看降维聚类分群，然后可视化多发性骨髓瘤领域的特殊基因TNFRSF17的表达量，如下所示：
特殊基因TNFRSF17的表达量
居然真的是恶性肿瘤细胞和正常的血液免疫细胞如此的泾渭分明，不是文章瞎搞，上面的UMAP图里面的左边部分全部是肿瘤细胞然后右边部分是正常细胞。而且通过常规基因我们可以给右边的正常的血液免疫细胞进行简单的命名，比如t细胞和b细胞还有髓系免疫细胞。
右边的正常的血液免疫细胞进行简单的命名
但是，如果我们这最开始的降维聚类分群流程里面的选择了对这些样品进行了harmony整合，就会出现很诡异的结果：
出现很诡异的结果
这个时候不要说是肿瘤病人异质性问题了，连肿瘤细胞和正常细胞都会被混合在一起，因为这个是血液肿瘤所以都是免疫细胞，它们这被harmony抹除样品个体差异的时候就顺便把肿瘤细胞和正常细胞的差异也给抹除了。
那么肿瘤病人应该是第一层次降维聚类分群就不harmony吗？
可以看到其它类似的肿瘤病人的单细胞转录组数据如果第一层次降维聚类分群就不harmony那么结果也是病人独立成群，比如《Single-cell analyses of renal cell cancers reveal insights into tumor microenvironment, cell of origin, and therapy response》也是如此
病人独立成群
如果多个样品进行了整合，也就是说第一层次降维聚类分群采用了harmony，一般来说就不会有病人特异性独立成一个群的现象，比如《Evolution of core archetypal phenotypes in progressive high grade serous ovarian cancer》：
不会有病人特异性独立成一个群的现象
也就是说harmony在抹除病人异质性这方面的效果是杠杠的，比如2021的NC文章里面的25个病人也是被整合的非常干净：《An N-Cadherin 2 expressing epithelial cell subpopulation predicts response to surgery, chemotherapy and immunotherapy in bladder cancer》
25个病人也是被整合的非常干净
上皮细胞里面的有一些细分亚群是恶性的肿瘤细胞，需要跑inferCNV流程才有可能区分哦，详见：肿瘤单细胞转录组拷贝数分析结果解读和应用，而且很明显这里面的降维聚类分群理论上是跟inferCNV流程的拷贝数冲突的，并不是说上面的5个亚群要么是恶性要么是正常二倍体，而是5个亚群本身在拷贝数变化层面仍然是有异质性！类似于下面的效果：
亚群本身在拷贝数变化层面仍然是有异质性
上面的画圈的两个亚群里面都是部分细胞是恶性的，而且有意思的是这两个降维聚类分群后的完全不同的亚群里面的部分细胞居然是拷贝数可以融合，这两个维度的分群是可以冲突的。
