分享是一种态度
这篇文章《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基因的背景知识:
下面来看看控制一个参数变量改变的影响吧:
具体参数:
具体参数:
可以看到在高变基因3000的时候,T细胞和B细胞是分开的。然而在高变基因2000的时候,T细胞和B细胞是连接在一起的。
具体参数:
具体参数:
在PCA维数为110和30的情况下,T细胞和B细胞是连接着的。但是PCA维数为110的时候,T细胞分为了两个部分,而维数为30的时候,T细胞是一个部分。
具体参数:
具体参数:
这个参数的影响不大
这个参数的影响挺大的。min_dist为0.01的时候,可以将T细胞和B细胞分开。
可以看到在dims参数为1:27的时候,T细胞和B细胞是分开的。然而在dims参数为1:15的时候,T细胞和B细胞是连接在一起的。
从这次的结果看,增加高变基因数目会把T细胞和B细胞分开。而PCA维数的增加没有使T细胞和B细胞分开,但是把T细胞分成了两个部分。
之前在这篇推文初探单细胞分析 — 标准化与降维聚类分群的理解提到过:高变基因个数越多,PCs的数目越多,保留的数据信息也就越多,更有可能引入噪音,运行速度也越慢。但是也不能太少,否则会丢失很多数据信息。
所以有可能是增加的高变基因使得T细胞和B细胞之间有更多的差异,可以分开。而PCA维数的增加把T细胞分成了两个部分,可能是由于有噪音引入。
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::fread( file.path(dir,pro),
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 <- JoinLayers(sce.all)
dim(sce.all[["RNA"]]$counts )
as.data.frame(sce.all@assays$RNA$counts[1:10, 1:2])
head(sce.all@meta.data, 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 <- NormalizeData(sce.all.filt,
normalization.method = "LogNormalize",
scale.factor = 1e4)
## 评估参数的影响
variable_nums <- c(2000,3000,4000)
npcs_nums <- c(30,50,70,90,110)
neighbors_nums <- c(20,30,40,50)
min_dist_nums <- c(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 <- FindVariableFeatures(sce.all.filt, selection.method = "vst",
nfeatures = i)
tmp_sce <- ScaleData(tmp_sce)
tmp_sce <- RunPCA(tmp_sce, features = VariableFeatures(object = tmp_sce),
verbose = FALSE, npcs= j)
tmp_sce <- RunHarmony(tmp_sce, "orig.ident")
tmp_sce <- FindNeighbors(tmp_sce, dims = 1:10, verbose = FALSE,reduction = "harmony")
tmp_sce <- FindClusters(tmp_sce, resolution = 0.5, verbose = FALSE)
table(tmp_sce$seurat_clusters)
for (n in neighbors_nums){
for (d in min_dist_nums){
tmp_sce <- RunUMAP(tmp_sce, dims = 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 <- DimPlot(tmp_sce,label = 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 <- DotPlot(tmp_sce , features = genes_to_check ) +
coord_flip() +
theme(axis.text.x=element_text(angle=45,hjust = 1))
print(p2)
p3 <- FeaturePlot(object = tmp_sce, features = T_genes,raster=F)
print(p3)
p4 <- FeaturePlot(object = tmp_sce, features = B_genes,raster=F)
print(p4)
dev.off()
}
}
}
}