当你想要单细胞分出特定的群数时，试试FindCluster2

获取代码和更佳阅读体验获取，请移步：当你想要单细胞分出特定的群数时，试试FindCluster2

前言

在做单细胞分群时，有时候我们根据文献或者偏好已经有的了预想的细胞群数，那就需要对分辨率一点点地调整，为了更简便地自动化地实现这一过程，小编开发了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)
}