dafabet手机黄金版_dafabet黄金手机版

在诸多单细胞数据分析的过程中，对于聚类得到的细胞亚群是否可以进一步细分对于后续细胞类型定义息息相关。本文从背景知识、算法原理、性能评估、应用案例、代码与结果展示五个方面，将理论与实践结合，图文并茂的详细介绍了解决该问题的一个算法工具。

单细胞聚类过程涉及很多参数，诸如PCA，Resolution等等，这些参数的选择对于亚群的分类颗粒度及数目多少息息相关。如果数据是过度聚类的，那么许多类群纯粹是由技术噪音驱动的，并不能反映不同的生物状态。如果数据是欠聚类的，细胞亚型之间的差异信号很可能被忽略而失去细胞亚型异质性的判断。现有的评估聚类质量的工具，如广泛使用的silhouettecoefficient（轮廓系数），不能揭示聚类内的差异性是由于亚群的存在还是随机噪声的干扰。因此，只有合理的参数选择才会使得我们对不同亚群细胞类型的定义更加精准。那么如何判断当前参数下，聚类得到的细胞亚群结果如何呢？

来自荷兰莱顿大学的副教授StefanSemrau课题组开发了名为Sigma（SIGnal-Measurement-Angle）的工具，并于2021年5月11号发表在预印刊bioRxiv上，MirceaM ,  Hochane M ,  Fan X , et al. A clusterability measure forsingle-cell transcriptomics reveals phenotypic subpopulations.  2021.

该作者从随机矩阵理论和扰动理论的视角将我们获得的细胞基因表达量矩阵视为真实细胞基因表达信号对随机噪声信号矩阵的扰动。

接下来，我们来讲解一下该工具的算法原理。在说明其原理之前，我们先来简单介绍一下其中用到的一个常见算法，奇异值分解（SingularValue Decomposition简称SVD），该方法在PCA、异常检测、图像去噪、推荐系统等等领域都有广泛的应用。

图2：不同奇异值选择项下的信号情况

图2中，为日本著名女星上野树里的黑白照片，高宽为450*333，我们将A记为该图片的像素值矩阵。我们可以对A进行奇异值分解如公式（1）所示，将其分解为若干秩一矩阵之和，并根据σ大小进行排序，使得，其中σ即为对应每个对应秩一矩阵的奇异值。

......（1）

我们在图2中由左到右分别取前1、5、20、50项，我们可以看到得到的图像信号将越来越接近原始图像信号。因此奇异值往往对应着矩阵中隐含的重要信息，且重要性和奇异值大小正相关。我们可以理解为，将一张张不同厚度的特定规整（对应秩一的矩阵）的纱窗按照厚度从厚到薄平铺到桌面上，那么我们只需要前几十张就可以清晰的还原出原始图像的样子！而对应越厚的纱窗包含越多越重要原始图像的信息，越薄的纱窗越可有可无，将更被认为是噪声。

也就是说，如果有一副图像包含有噪声，我们会认为那些较小的奇异值就是对应的噪声矩阵，而重要的图像信息矩阵往往对应着更大的奇异值。

图3a中，我们可将测序得到的细胞基因表达量矩阵看做随机噪声矩阵（noise）和信号矩阵(signal)的叠加，根据随机矩阵理论，随机高阶矩阵的奇异值服从Marchenko-Pastur（MP）分布，而在MP分布之外TracyWidom(TW)阈值以上的奇异值则对应于细胞真实的表达信号。那么阈值之上的奇异值向量与MP分布中奇异值向量的最小夹角的余弦值的平方即为SIGMA值，因此SIGMA值在[0,1]范围内取值，越接近1，则说明该亚群能够被进一步分群的可能性越大，越接近于0，则说明该亚群信噪比已经很小，噪声相对信号已经很强烈，越难被各种算法进一步细分。

图3：SIGMA算法理论框架与理论评估

为了评价算法的性能，作者提出利用兰德系数（AdjustedRand Index[ARI]），如公式（2）所示，作为一个很好的评价方法，来评估聚类结果与真实分群之间的一致性。理论上可实现的ARI(TheoreticallyachievableTARI)可以从对贝叶斯错误率推导出来。因此，作者通过检验TARI这种评估方法与SIGMA评估结果的相关性得到如图2d的曲线，显示出两种评估方法的强烈相关性。作者同时也对比了与之前发表的单细胞聚类纯度算法ROGUE（如图4所示）进行了对比，在模拟的数据集上用一系列主流的不同参数的算法进行聚类，并评估理论轮廓系数与ROGEU以及SIGMA的相关性，如图5所示，结果发现ROGUE与SIGMA评价结果并不一致，作者解释为可能ROGUE是从不同的视角来解释聚类纯度的问题。同时，作者使用真实实验数据的不同比例混样对SIGMA进行评估，如图6所示，结果显示其依然明显优于其他指标。

......（2）

图4：ROGUE算法评估亚群纯度

图5：最佳轮廓系数的近似上限和ROGUE与tARI的一致性

图6：SIGMA在实验数据上的评价结果优于其他测量方法

作者将SIGMA算法应用到先前已发表文章的数据集中，尝试发现新的生物学亚群并阐释新的未被挖掘出的生物学意义。

作者首先将其应用在骨髓单核细胞（bonemarrow mononuclearcells）中进行细胞亚群的纯度分析，如图7所示，对每个亚群进行SIGMA值计算，发现有很多亚群都有较高的SIGMA值，即提示有进一步细分的可能性，作者选取了SIGMA值排名最大的前两个亚群进行进一步聚类细分，在红细胞(RBC)祖细胞中，重新鉴定出了4个亚群，它们分别对应于不同的分化阶段，从红细胞前体细胞到高度分化的红细胞，这些细胞亚群也被先前文献的研究证实。在树突状细胞(DC)祖细胞簇中，发现了两个亚群，它们分别对应于经典树突状细胞和浆细胞样树突状细胞的前体细胞。

图7：将SIGMA应用于BMNC数据可以推动发现具有生物学意义的子簇

作者又将SIGMA应用于他们课题组之前发表过的胎儿人类肾脏数据集（fetalhumankidney），如图8所示，并对亚群中SIGMA值得分最高的输尿管芽/集管(UBCD)进行进一步聚类分群，并计算出亚群的标记基因，并通过免疫荧光共定位的实验来进行验证，如图9所示，结果显示免疫荧光染色与计算得到的标记基因表达具有高度的一致性。

基于以上论证，作者认为SIGMA，是一种有助于评估细胞亚群聚类可及性的度量方法，能够在单细胞转录组数据中重新发现可能被忽视的微小的表达信号，同时，还能够提供一个新的视角来解释差异分群的驱动基因。

图8：SIGMA应用于人胎儿肾导致发现有生物学意义的亚簇

图9：免疫荧光染色显示SIGMA值最高的聚类具有一致性

作者提供的Github官方链接为：https://github.com/Siliegia/SIGMA，其中简要介绍了几个主要函数。我们基于此，开发参考代码如下，从外部输入对象project.rdata，该对象包含由细胞基因表达量矩阵构建的标准Seurat对象，主要用到其细胞表达量矩阵project@assays$RNA@counts和对应的分群信息project$Cluster，当然也可以将其替换为其他存储形式的文件直接读入。结果如图10~13所示：

devtools::install_github("Siliegia/SIGMA")

library(SIGMA)

library(getopt)

library(Seurat)

parameter=matrix(c(

"project.rdata","p","1","character","Seuratproject rdata,required",

"celltype","c","1","character","celltypeof each cells,optional"#,

#"cluster_sc","c","1","character","clusters_SC.csv,need"

),byrow = TRUE,ncol = 5)

opt=getopt(parameter)

if (!is.null(opt$project.rdata)){

load(opt$project.rdata)

expr=as.matrix(project@assays$RNA@counts)

expr.norm <- t(t(expr)/colSums(expr))*10000

expr.norm.log <- log(expr.norm + 1)  

order_cells=data.frame(Barcode=colnames(expr.norm.log),Cluster=as.numeric(project$Cluster))

order_cells2=order_cells[order(order_cells$Cluster),]

expr.norm.log2=expr.norm.log[,as.character(order_cells2$Barcode)]

expr2=expr[,as.character(order_cells2$Barcode)]

out2 <- sigma_funct(expr = expr.norm.log2,clusters = as.character(order_cells2$Cluster),

exclude = data.frame(clsm =log(colSums(expr2) + 1)))

pdf("sigma_plot.pdf")

plot_sigma(out2)

dev.off()

pdf("sigma_all_plot.pdf")

plot_all_sigmas(out2)

dev.off()

pdf("g_sigma_all_plot.pdf")

plot_all_g_sigmas(out2)

dev.off()

for (i in unique(project$Cluster)) {

data_temp=get_var_genes(out2, i)

if (data_temp=="No further clusters"){

write.table(file=paste0("Cluster",i,"_variance_driving_genes.xls"),data_temp)

}else{

write.table(file=paste0("Cluster",i,"_variance_driving_genes.xls"),data.frame(Order=rownames(data_temp),data_temp),sep='\t',quote= F,row.names = F)

}

}

pdf("MP_distribution_fits_to_the_data.pdf")

for (i in unique(project$Cluster)) {

plot_MP(out2, i,plot.title =paste0("Cluster",i))

}

dev.off()

pdf("plotting_clusters_with_the_singularvectors.pdf")

for (i in unique(project$Cluster)) {

plot_singular_vectors(out2, i, colour =project$Cluster[project$Cluster == i])

}

dev.off()

}else{

print("lack of project.rdata file")

}

图10：各聚类亚群的SIGMA值

图11：对应亚群的奇异值分布

图12：根据奇异值对亚群可分性的展示

图13：奇异值视角差异信号驱动基因

参考文献

[1] Mircea M ,  Hochane M ,  Fan X , et al. Aclusterability measure for single-cell transcriptomics revealsphenotypic subpopulations.  2021.

[2]https://www.biorxiv.org/content/10.1101/2021.05.11.443685v1.full.pdf ###paper availability

[3] https://github.com/Siliegia/SIGMA  ###Codeavailability

[4]https://mp.weixin.qq.com/s?__biz=MzU2MDgyNzg3NQ==&mid=2247484393&idx=1&sn=ad6f6c7097c586a19b4a1a1d347077aa&scene=19#wechat_redirect####相关报道说明

[5]“An entropy-based metric for assessing thepurity of single cell populations.”  ##文章来源及路径

[6]https://www.ncbi.nlm.nih.gov/pmc/articles/pmid/32572028/