ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

单细胞转录组聚类分析R包——SC3的应用

2021-10-14 21:59:37  阅读:378  来源: 互联网

标签:www 单细胞 转录 https 聚类 SC3 聚类分析


前言

在做单细胞相关聚类分析的过程中,为了研究异质性,发现了一种专门用于单细胞转录组的聚类分析R包——SC3(single cell consensus clustering,单细胞一致性聚类)。该聚类算法发表于Nature Methods(2017),能够让单细胞RNA-seq基于转录组特征对细胞类型进行定量表征。这是一种用户友好的无监督聚类工具,它通过一致方法将多个聚类解决方案组合在一起,从而得到高精度和鲁棒性的分群结果。SC3包的官方具体使用教程见:http://www.bioconductor.org/packages/release/bioc/vignettes/SC3/inst/doc/SC3.html
SC3包源代码见Github:https://github.com/hemberg-lab/SC3

原文的SC3聚类流程如图所示:SC3聚类的一般流程
(a)使用SC3框架进行聚类的概述(参见方法)。用Treutlein数据举例说明了一致步骤。
(b)用于设置SC3参数的已发布数据集。N是数据集中的细胞数;k为作者最初确定的簇数;单位:RPKM是每千碱基每百万次读的转录本,RPM是每百万次读的转录本,FPKM是每千碱基每百万次读的转录本片段,TPM是每百万次读的转录本片段。
( c)ARI>处d值的直方图。金标准数据集达到95。黑色竖线表示细胞总数N的d = 4-7%,分类准确率高。
(d) (b)中所示数据集的SC3聚类的100个实现。条对应点的中位数。红色和灰色分别对应有一致步长和无一致步长时的聚类。这条黑线对应的ARI=0。8。黑色虚线分隔了金和银标准数据集。

笔者主要想通过该包对单细胞转录组进行聚类数目的估计。在网上没有直接写SC3关于这个功能的帖子,故我去原文,官方教程以及代码研究了一下,该包主要通过随机矩阵理论(Random Matrix Theory,RMT) 确定聚类的数目。

安装SC3包

if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("SC3")

同时还需要安装如下包:

library(SingleCellExperiment)#安装方法和SC3一样
library(SC3)
library(ggplot2)
library(scater)

利用SC3实现聚类数目的估计

首先看一下源代码:

if (k_estimator) {
        object <- sc3_estimate_k(object)
        # Do not override cluster if user has set a k
        if (is.null(ks))
        {
            ks <- metadata(object)$sc3$k_estimation
        }

把该方法应用于自己的代码:
注意:在跑自己的数据之前,要先准备两个文件:(1)单细胞表达矩阵文件;(2)每个细胞的注释文件。

pbmc <- read.table("SC3/分组表达/exp_SCE_M.csv",stringsAsFactors=FALSE, header=TRUE, check.names=FALSE,row.names=1, sep=",")
ann<-read.table("SC3/ann/SCE_M.csv", header=T, row.names=1, sep=",")
#构建sce对象
sce <- SingleCellExperiment(
    assays = list(
        counts = as.matrix(pbmc),
        logcounts = log2(as.matrix(pbmc) + 1)
    ),
	colData=ann
	)

# define feature names in feature_symbol column
rowData(sce)$feature_symbol <- rownames(sce)
# remove features with duplicated names
sce <- sce[!duplicated(rowData(sce)$feature_symbol), ]
# define spike-ins
#isSpike(sce, "ERCC") <- grepl("ERCC", rowData(sce)$feature_symbol)

#开始聚类(estimate k)
sce<-sc3_estimate_k(sce) #估计k值
ks<-metadata(sce)$sc3$k_estimation #调用k值
ks
a <- ks - 1
b <- ks + 1
sce <- sc3(sce, ks = a:b, biology = TRUE) #这一步会运行很久

#主成分分析PCA
sce <- runPCA(sce)
plotPCA(sce, colour_by = "cell_type")

#画降维后的分布
sc3_n_clusters <- paste0("sc3_", as.character(ks),"_clusters")
sc3_n_log2_outlier_score <- paste0("sc3_", as.character(ks), "_log2_outlier_score")
sc3_n_clusters
sc3_n_log2_outlier_score
plotPCA(
    sce, 
    colour_by = sc3_n_clusters, 
    size_by = sc3_n_log2_outlier_score
)

#计算一致性矩阵
sc3_plot_consensus(
    sce, k = ks, 
    show_pdata = c(
        #"cell_type", 
        #"log10_total_features",
        sc3_n_clusters, 
        sc3_n_log2_outlier_score
    )
)

#计算簇内稳定性
sc3_plot_cluster_stability(sce, k = ks)

绘制的PCA结果和一致性矩阵结果如下图:
在这里插入图片描述
在这里插入图片描述
SC3包后续还可以进行基因表达矩阵热图的绘制; 提供调整p值< 0.01的所有差异表达基因的列表,并绘制p值最低的50个基因的基因表达谱;以及找到标记基因。

参考:
代码:https://blog.csdn.net/weixin_42889479/article/details/108943720
https://www.jianshu.com/p/65094dfd6922
https://www.jianshu.com/p/24f65dd2c88f
https://www.jianshu.com/p/ea81c3635a8b

标签:www,单细胞,转录,https,聚类,SC3,聚类分析
来源: https://blog.csdn.net/m0_51042606/article/details/120772013

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有