harmony整合“失败”才是正常的

我们的马拉松授课的单细胞环节给了大家一套适用于几千个单细胞转录组数据集的完整降维聚类分群代码，但是有一些小伙伴在应用于自己感兴趣的研究的时候就出现了一些看起来是bug的结果。

比如咱们九月份的学员就反馈了这个数据集：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE125688

GSM3580367 BEC_ctr1
GSM3580368 BEC_ctr2
GSM3580369 BEC_ctr3
GSM3580370 BEC_DDC1
GSM3580371 HEP_ctr1
GSM3580372 HEP_DDC1

对应的文章是：Single-Cell Analysis of the Liver Epithelium Reveals Dynamic Heterogeneity and an Essential Role for YAP in Homeostasis and Regeneration. Cell Stem Cell 2019 Jul 3;25(1):23-38.e8. PMID: 31080134

这个研究关心的是这两个不同的上皮细胞： hepatocytes and biliary epithelial cells (BECs)，我们很容易下载到不同样品的作者提供的表达量矩阵文件：

 ls -lh *gz|cut -d" " -f 6-
  1.7M  1 26  2019 GSM3580367_bil-adult1.coutb.csv.gz
  2.3M  1 26  2019 GSM3580367_bil-adult1.coutc.csv.gz
  1.8M  1 26  2019 GSM3580367_bil-adult1.coutt.csv.gz
  1.6M  1 26  2019 GSM3580368_bil-adult2.coutb.csv.gz
  2.1M  1 26  2019 GSM3580368_bil-adult2.coutc.csv.gz
  1.7M  1 26  2019 GSM3580368_bil-adult2.coutt.csv.gz
  2.3M  1 26  2019 GSM3580369_bil-adult3.coutb.csv.gz
  3.2M  1 26  2019 GSM3580369_bil-adult3.coutc.csv.gz
  2.5M  1 26  2019 GSM3580369_bil-adult3.coutt.csv.gz
  3.5M  1 26  2019 GSM3580370_bil-DDC1.coutb.csv.gz
  4.7M  1 26  2019 GSM3580370_bil-DDC1.coutc.csv.gz
  3.9M  1 26  2019 GSM3580370_bil-DDC1.coutt.csv.gz
  6.8M  1 26  2019 GSM3580371_hep-adult1.coutb.csv.gz
  8.7M  1 26  2019 GSM3580371_hep-adult1.coutc.csv.gz
  8.5M  1 26  2019 GSM3580371_hep-adult1.coutt.csv.gz
  5.6M  1 26  2019 GSM3580372_hep-DDC1.coutb.csv.gz
  7.5M  1 26  2019 GSM3580372_hep-DDC1.coutc.csv.gz
  6.8M  1 26  2019 GSM3580372_hep-DDC1.coutt.csv.gz

很简单的批量读取代码：

dir='GSE125688_RAW/'
samples=list.files( dir  ,pattern = 'c.csv.gz')
samples
sceList = lapply(samples,function(pro){ 
  # pro=samples[3] 
  print(pro)   
  tmp=fread( file.path(dir,pro )   ,
              header = F,data.table = F ) 
  head(tmp)
  gid=str_split(tmp[,1],'_',simplify = T)[,1]
  kp=!duplicated(gid);table(kp)
  tmp=tmp[kp,]
  gid=gid[kp]
  rownames(tmp)=gid
  ct=tmp[,-1] 
  sce =CreateSeuratObject(counts =  ct ,
                          project =  gsub('.coutc.csv.gz','',
                                          gsub('GSM[0-9]*_','',pro) )  ,
                          min.cells = 5,
                          min.features = 300 )
  return(sce)
}) 
do.call(rbind,lapply(sceList, dim)) 
samples=gsub('.coutc.csv.gz','',
     gsub('GSM[0-9]*_','',samples) ) 
samples
sce.all=merge(x=sceList[[1]],
              y=sceList[ -1 ] ,   
              add.cell.ids = samples  )  
sce.all <- JoinLayers(sce.all) 

降维聚类分群后可以看到：

               0    1    2    3
  bil-adult1    5 1236    0    0
  bil-adult2    0 1369    0    0
  bil-adult3   12 1491   11    0
  bil-DDC1     38 1688    3   23
  hep-adult1 4354    0   21    4
  hep-DDC1   2563    8   54   58

如果是仅仅是看上面的表格，初学者确实是很容易错误的判断了harmony整合“失败”。比如我们看2023-GSE181919-头颈癌-疾病进展单细胞，成功的整合后应该是如下所示 ：

其实是因为，初学者对生物学背景的把握不行，这个时候harmony整合“失败”才是正常的，既然我们的样品本来就是纯粹的胆管上皮细胞和肝细胞的，那么样品之间就不能被harmony混合， 合理的结果，说明了我们的harmony可以很好的处理样品差异和真实的生物学差异问题！

因为很容易根据特异性基因去给这些亚群生物学名字，很明显1是胆管上皮细胞，而0是肝细胞，所以它们有样品特异性 ：