实战探究五个参数对UMAP图可视化的影响

这篇文章《 A comprehensive single-cell map of T cell exhaustion-associated immune environ- ments in human breast cancer 》的UMAP图中T细胞和B细胞是分开的，但之前复现的时候发现T细胞和B细胞是连接在一起的。有很多小伙伴纠结于为什么我复现的可视化图和原文献很不一样呢。其实也不用太紧张，有的参数只是影响了肉眼但不会影响细胞的本质的属性（也就是分群）。但是如果想较真还是可以探索一下参数的影响滴。

这篇推文的目的是探索一些重要参数对后续分群UMAP可视化的影响。参数主要考虑： 高变基因个数；pca维数；UMAP中的n_neighbors，min_dist和dims参数 。影响主要看T细胞和B细胞是否分开。

marker基因的背景知识：

• T细胞（CD3D,CD3E）
• B细胞（MS4A1,CD79A）

参数定义

• FindVariableFeatures 的 nfeatures 参数（高变基因个数）
• RunPCA 的 npcs 参数（PCA的维数）
• RunUMAP 的 n.neighbors 和 min.dist 参数
• n_neighbors -- 用于构建初始高维图的近似最近邻点的数量。它有效地控制了UMAP如何平衡局部结构和全局结构--低值会促使UMAP在分析高维数据时，通过约束考虑的邻点数量，更加关注局部结构，而高值则会促使UMAP倾向于表示大局结构，同时失去精细的细节。
• min_dist -- 是低维空间中点之间的最小距离。这个参数控制了UMAP将点紧密地聚集在一起的程度，数值越低，嵌入的点就越紧密。更大的min_dist值会使UMAP将点更松散地包装在一起，而专注于保存广泛的拓扑结构。
• dims -- 考虑的特征维数（Which dimensions to use as input features, used only if features is NULL）

下面来看看控制一个参数变量改变的影响吧：

高变基因的影响

高变基因个数为3000：

具体参数：

• 高变基因3000；
• pca维数110；
• UMAP参数：n_neighbors为50，min_dist为0.1，dims为1:15

高变基因个数为2000：

具体参数：

• 高变基因2000；
• pca维数110；
• UMAP参数：n_neighbors为50，min_dist为0.1，dims为1:15

可以看到在高变基因3000的时候，T细胞和B细胞是分开的。然而在高变基因2000的时候，T细胞和B细胞是连接在一起的。

PCA维数的影响

PCA维数（npcs）为110：

具体参数：

• 高变基因3000；
• pca维数110；
• UMAP参数：n_neighbors为30，min_dist为0.3，dims为1:15

PCA维数（npcs）为30：

具体参数：

• 高变基因3000；
• pca维数30；
• UMAP参数：n_neighbors为30，min_dist为0.3，dims为1:15

在PCA维数为110和30的情况下，T细胞和B细胞是连接着的。但是PCA维数为110的时候，T细胞分为了两个部分，而维数为30的时候，T细胞是一个部分。

UMAP参数中的n_neighbors的影响

n_neighbors为50：

具体参数：

• 高变基因3000；
• pca维数50；
• UMAP参数：n_neighbors为50，min_dist为0.3，dims为1:15

n_neighbors为20：

具体参数：

• 高变基因3000；
• pca维数50；
• UMAP参数：n_neighbors为20，min_dist为0.3，dims为1:15

这个参数的影响不大

UMAP参数中的min_dist的影响

min_dist为0.01：

min_dist为0.5：

这个参数的影响挺大的。min_dist为0.01的时候，可以将T细胞和B细胞分开。

UMAP中dims参数的影响

dims参数为1:27：

dims参数为1:15：

可以看到在dims参数为1:27的时候，T细胞和B细胞是分开的。然而在dims参数为1:15的时候，T细胞和B细胞是连接在一起的。

小结

会影响分群和UMAP可视化的参数：

从这次的结果看，增加高变基因数目会把T细胞和B细胞分开。而PCA维数的增加没有使T细胞和B细胞分开，但是把T细胞分成了两个部分。

之前在这篇推文 初探单细胞分析 — 标准化与降维聚类分群的理解 提到过：高变基因个数越多，PCs的数目越多，保留的数据信息也就越多，更有可能引入噪音，运行速度也越慢。但是也不能太少，否则会丢失很多数据信息。

所以有可能是增加的高变基因使得T细胞和B细胞之间有更多的差异，可以分开。而PCA维数的增加把T细胞分成了两个部分，可能是由于有噪音引入。

不会影响分群，只会影响UMAP可视化的参数：

UMAP参数中的n_neighbors，min_dist和dims这三个参数，是不会影响细胞的本质属性的，只是影响UMAP可视化的图。可以看到减小的min_dist和增加的dims会把T细胞和B细胞分开。n_neighbors的影响不大。

希望大家看完之后对参数改变对UMAP图的影响有一些大致的理解。

附实例代码

rm(list=ls())
options(stringsAsFactors = F) 
source('../scRNA_scripts/lib.R')
getwd()

###### step1: 导入数据 ######   
# 付费环节 800 元人民币
# 参考：https://mp.weixin.qq.com/s/tw7lygmGDAbpzMTx57VvFw
dir='../inputs/'
samples=list.files( dir )
samples = samples[str_detect(samples,"singlecell_count_matrix.txt")]
sceList = lapply(samples,function(pro){ 
  # pro=samples[7] 
  print(pro)  
  ct                            data.table = F)
  ct[1:4,1:4]
  rownames(ct)=ct[,1]
  ct=ct[,-1]
  #ct=t(ct) 
  sce =CreateSeuratObject(counts =  ct ,
                          project =  gsub('_singlecell_count_matrix.txt.gz','',pro)  ,
                          min.cells = 5,
                          min.features = 300 )
  return(sce)
}) 
do.call(rbind,lapply(sceList, dim))
sce.all=merge(x=sceList[[1]],
              y=sceList[ -1 ],
              add.cell.ids = samples  ) 
names(sce.all@assays$RNA@layers)
sce.all[["RNA"]]$counts 
# Alternate accessor function with the same result
LayerData(sce.all, assay = "RNA", layer = "counts")
sce.all dim(sce.all[["RNA"]]$counts )

as.data.frame(sce.all@assays$RNA$counts[1:10, 1:2])
head([email protected], 10)
table(sce.all$orig.ident) 

library(stringr)
sce.all$orig.ident=str_split(colnames(sce.all),'[-_]',simplify = T)[,1]
table(sce.all$orig.ident) 
sce.all

# 如果为了控制代码复杂度和行数 
# 可以省略了质量控制环节
###### step2: QC质控 ######
dir.create("./1-QC")
setwd("./1-QC")
# 如果过滤的太狠，就需要去修改这个过滤代码
source('../../scRNA_scripts/qc.R')
sce.all.filt = basic_qc(sce.all)
print(dim(sce.all))
print(dim(sce.all.filt))
setwd('../')
sp='human'
getwd()

sce.all.filt                            normalization.method = "LogNormalize",
                           scale.factor = 1e4) 

## 评估参数的影响
variable_nums 2000,3000,4000)
npcs_nums 30,50,70,90,110)
neighbors_nums 20,30,40,50)
min_dist_nums 0.001,0.01,0.1,0.3,0.5) 
# https://zhuanlan.zhihu.com/p/352461768
# RunUMAP: n.neighbors default is 30[5-50]; min.dist default is 0.3[0.001-0.5]
if(F){
variable_nums=2000
npcs_nums=50
neighbors_nums=30
min_dist_nums=0.001
}

for (i in variable_nums){
  for (j in npcs_nums){
    tmp_sce "vst",
                                    nfeatures = i)
    tmp_sce     tmp_sce                       verbose = FALSE, npcs= j)
    tmp_sce "orig.ident")
    tmp_sce 1:10, verbose = FALSE,reduction = "harmony")
    tmp_sce 0.5, verbose = FALSE)
    table(tmp_sce$seurat_clusters)
    
    for (n in neighbors_nums){
      for (d in min_dist_nums){
        tmp_sce 1:15, umap.method = "uwot", metric = "cosine", reduction = "harmony", n.neighbors=n, min.dist=d)
        
        pdf(paste0("pca_object_v",i,".pca_",j,"neighbors_",n,"dist_",d,".pdf"))
        p1 T,reduction = "umap",raster=F)
        print(p1)
        genes_to_check = c("CD3E","CD3D","MS4A1","CD79A")
        T_genes = c("CD3E","CD3D")
        B_genes = c("MS4A1","CD79A")
        p2         coord_flip()  +
        theme(axis.text.x=element_text(angle=45,hjust = 1))
        print(p2)
        p3 F) 
        print(p3)