在做单细胞分群时,有时候我们根据文献或者偏好已经有的了预想的细胞群数,那就需要对分辨率一点点地调整,为了更简便地自动化地实现这一过程,小编开发了FindCluster2
函数。一起来学习下开发过程吧。
首先需要指定想要的cluster数或范围,在默认起始分辨率下运行FindCluster
计算当前分辨率下的细胞群数,并于指定的范围进行比较,决定下一步是增加还是减少分辨率还是退出循环,这种数值渐变且低于/高于某值的循环当然就是使用while
循环了,简单写了一下:
FindClusters2 <- function(obj, cluster = c(35, 40), verbose = FALSE){
if(verbose)
cat("Find suitable resolution, start with 1.\n")
nCluster = res = 1
while(nCluster < cluster[1] | nCluster > cluster[2]){
if(verbose)
cat("resolution", res, "... ")
obj <- FindClusters(obj, resolution = res, verbose = FALSE)
nCluster = length(unique(obj$seurat_clusters))
if(nCluster < cluster[1]){
res = res + 0.1
} else if (nCluster > cluster[2]){
res = res - 0.1
} else{
break
}
if(verbose)
cat(nCluster, "clusters.\n")
}
obj[["best.res"]] = res
if(verbose)
cat("Final resolution:", res, "with", nCluster, "clusters.\n")
resL = list("obj" = obj, "best.res" = res)
return(resL)
}
默认分辨率设为1,默认步长为0.1。
经过测试,对付一般情况下,足够使用了。但仔细想一想这里代码有两个大的漏洞。
第一,分辨率取值有误。我们知道分辨率的取值范围是大于0的,但是我们代码每个循环都减去固定的一个值,那当指定的细胞群数很少时,需要的分辨率小于0.1时,则分辨率将继续减去0.1,就出bug了。
因此应该控制分辨率的取值范围要大于0,这让我想到了逻辑斯蒂方程,其取值范围时0-1,那我们再乘以10就得到了0-10取值范围的分辨率值,足够我们使用了。
"定义逻辑斯蒂函数"
sigmoid <- function(x) 1 / (1 + exp(-x))
"根据输入值计算分辨率"
res = signif(sigmoid(x) * 10, 3)
第二,有死循环的风险。当指定的细胞群数范围较小或步长较大时,指定的范围有可能被跳过,这将会造成左右无限蹦迪的死循环现象,所以要增加个判断。首先想想,正常情况下,在判断当前细胞群数与指定细胞群数时,大于或小于的情况永远只会出现一种,如果都曾经出现则说明有跳过折返的情况,因此只要判断大于和小于的情况如果都出现过,则抛出错误,提示指定的范围被跳过,并建议扩宽细胞群数范围或减小步长。给每个判断语句下面加个计数器即可。
还有两点可以优化。第一,上面代码最终最佳分辨率是以列表的形式和原始对象输出,这是因为我开始没用找到将最佳分辨率加到对象里的方法。又搜了以下发现:搜索add slot to seurat
找到https://github.com/satijalab/seurat/discussions/5617
, 第一个方法是定义一个新类,这跟用列表封装也没啥区别了;第二个建议是加入到misc
这个slot里面,那么这是slot是干什么的?搜索misc slot seurat
得到https://mojaveazure.github.io/seurat-object/reference/Misc.html
,miscellaneous意思是杂七杂八的东西,所以这里你可以放任何东西。
类内元素的赋值提取也有两种方法,一种直接类似于列表的提取方法,另一种是使用类成员函数
。
"赋值"
obj@misc[["best.resolution"]] = res
Misc(obj, slot = "best.resolution") <- res
"提取"
obj@misc[["best.resolution"]]
Misc(obj, slot = "best.resolution")
死去的C++突然开始攻击我。。。
第二,可指定初始分辨率。对于经验丰富的人或者已经经过几轮筛选的情况下,可能已经有了大概的分辨率的取值范围,那从这个值开始计算的话,速度会快很多。因为我们是通过逻辑斯蒂方程计算的分辨率,那得到特定分辨率时的x值就要使用其反函数了,也很容易计算:
x = -log(10/res - 1)
FindClusters2 <- function(obj, cluster.range = c(35, 40), by = 0.1, res = 1, verbose = FALSE){
if(verbose)
cat("Find suitable resolution, start with", res, "\n")
if(length(cluster.range) == 1){
cluster.range = c(cluster.range, cluster.range)
}
sigmoid <- function(x) 1 / (1 + exp(-x))
x = -log(10/res - 1)
plusCounter = minusCounter = 0
nCluster = 1
while(nCluster < cluster.range[1] | nCluster > cluster.range[2]){
if(verbose)
cat("resolution", res, "... ")
obj <- FindClusters(obj, resolution = res, verbose = FALSE)
nCluster = length(unique(obj$seurat_clusters))
if(nCluster < cluster.range[1]){
x = x + by
plusCounter = plusCounter + 1
} else if (nCluster > cluster.range[2]){
x = x - by
minusCounter = minusCounter + 1
} else{
break
}
res = signif(sigmoid(x) * 10, 3)
if(plusCounter & minusCounter){
cat("\n")
stop("Specific cluster ranger was skipped! Try expanding the cluster range or reducing the resolution step size.")
}
if(verbose)
cat(nCluster, "clusters.\n")
}
obj@misc[["best.resolution"]] = res
if(verbose)
cat("Final resolution:", res, "with", nCluster, "clusters.\n")
return(obj)
}
原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。
如有侵权,请联系 cloudcommunity@tencent.com 删除。
原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。
如有侵权,请联系 cloudcommunity@tencent.com 删除。