单细胞亚群鉴定知多少 第二篇

单细胞亚群鉴定

知多少（二）

上一期我们详细介绍了几种细胞marker基因的获取方法，为细胞亚群鉴定打好了坚实的基础。在这里我们接着看如何进行单细胞亚群鉴定，读取“单细胞这本天书”。

01

单细胞亚群鉴定

scMCA/scHCA

单细胞亚群鉴定：

scMCA/scHCA(http://bis.zju.edu.cn/MCA/blast.html)

在上一篇（单细胞亚群鉴定知多少 第一篇）介绍marker基因数据库的时候提到过scHCA数据库，这不仅仅可以用于单细胞marker基因的查找，而且郭国冀老师课题组还提供了单细胞类型预测的R工具包和在线工具。

在线的这里就不介绍了，这里介绍一下R工具包的安装和使用（scHCA）。由于小鼠的和人的R工具包使用类似，小鼠的就不再赘述。

1. scMCA/scHCA的安装

scMCA/scHCA的安装比较简单，由于该工具是GitHub上面的包，因此需要使用devtools进行安装，具体安装方法如下所示：

1#This require devtools
2install.packages('devtools') 
3library(devtools) # scHCL requires ggplot2/reshape2/plotly/shiny/shinythemes/shiny 
4install_github("ggjlab/scHCL") 
5#scMCA工具包的安装，小鼠图谱
6install_github("ggjlab/scMCA")

2. scMCA/scHCA的使用

 1library(Seurat)
 2library(scHCL)
 3###这里直接使用Seurat分析结果后的对象保存的rds文件
 4rds<-readRDS('cluster.ifnb.rds')
 5counts<-as.matrix(rds@assays$RNA@counts) ##官网建议的是Rawcounts
 6dim(counts)
 7#[1] 2884 80 
 8mca_result <- scHCL(scdata = counts, numbers_plot = 3)  # scMCA
 9str(mca_result)
10#List of 4
11#$ cors_matrix    : num [1:894, 1:80] 0.133 0.119 0.188 0.14 0.136 ...
12# ..- attr(*, "dimnames")=List of 2
13# .. ..$ : chr [1:894] "Oligodendrocyte Lineage .. ..$ : chr [1:80] "E18_2_C06" "E18_2_C07" "E18_2_C11" "E18_2_C12" ...
14#$ top_cors       : num 3
15# $ scMCA          : Named chr [1:80] "AT1 Cell(Lung)" "AT1 Cell(Lung)" "AT2 Cell(Lung)" "AT1 Cell(Lung)" ...
16# ..- attr(*, "names")= chr [1:80] "E18_2_C06" "E18_2_C07" "E18_2_C11" "E18_2_C12" ...
17# $ scMCA_probility:'data.frame':    240 obs. of  3 variables:
18# ..$ Cell     : Factor w/ 80 levels "E18_2_C06","E18_2_C07",..: 1 1 1 2 2 2 3 3 3 4 ...
19# ..$ Score    : num [1:240] 0.306 0.45 0.316 0.377 0.452 .

需要注意的是：我们这里使用的是Seurat包进行细胞亚群分析以后、pbmc对象的rds文件，因此直接读取rds文件，从rds文件中提取表达矩阵即可，后续的结果展示也直接从rds文件中提取相应的细胞亚群信息、样品信息。

下表格为输入的表达矩阵文件：

Gene	cell1	cell2	…	cellN
Gene1	1.00	0.00	…	0.00
………	………	………	…	………
GeneN	0.00	0.00	…	3.00

3. 结果展示

其实官方给了一个shinny的展现形式，但是由于我们一般10xGenomics单细胞转录组分析的时候是有很多细胞，我们关注的重点在于细胞亚群的细胞类型，因此我们可以使用另外一种展现形式。

1#shinny展现形式
2scHCL_vis(mca_result)

根据rds文件的细胞亚群注释信息，我们可以提供另外的一种展现形式，如下图所示。

 1#' @export
 2gettissue <- function(x,Num=3){
 3  top_markers <-order(x,decreasing = T)[1:Num]
 4  return(top_markers)
 5}
 6cors <- mca_result$cors_matrix
 7cors_index <- apply(cors,2,gettissue,numbers_plot)
 8cors_index <- sort(unique(as.integer(cors_index)))
 9data = cors[cors_index,]
10annotation_col<-data.frame(stim=[e-mail]$stim,Celltype=[e-mail]$seurat_clusters)
11rownames(annotation_col)<-rownames([e-mail])
12order_cells<-annotation_col %>% dplyr::mutate(.,barcod=rownames(.)) %>% dplyr::arrange(.,Celltype,stim)
13gaps_col<-c()
14m<-0
15for (i in 1:max(as.numeric(annotation_col[,c('Celltype')]))){
16  gaps_col<-c(gaps_col,table(annotation_col[,c('Celltype')])[i]+m)
17  m<-table(annotation_col[,c('Celltype')])[i]+m
18}
19library(pheatmap)
20library(RColorBrewer)
21pheatmap(data[79:95,as.vector(order_cells$barcod)],show_rownames=T,show_colnames=F,cluster_cols = F,cluster_rows= F,annotation_col = annotation_col,treeheight_row=0,cellheight=12,fontsize=8,cellwidth=520/length(rownames(annotation_col)),border=FALSE,gaps_col = gaps_col,color = colorRampPalette(c("grey", "white", "red"))(100),main = " ")

02

单细胞亚群鉴定

SingleR

SingleR V1.0.0之前的版本，只需要选择human或者mouse两个物种进行分析即可，无需设置参考数据集，参考数据集都是内置在包中，但是运行极慢，极其占用资源。V1.0.0版本参考数据集可以自己选择，参考数据集更友好，运行速度极快。

SingleR 为biocondutor的包，安装方便，直接使用BiocManager安装即可，其预测需要参考数据集和未知数据集，参考数据集官方给了7个数据集，可以直接使用。

SingleR：

http://www.bioconductor.org/packages/release/bioc/vignettes/SingleR/inst/doc/SingleR.html

1. SingleR的安装以及使用

SingleR安装比较方便，直接使用BiocManager安装即可，使用也比较方便，参考数据集直接使用官方给的参考数据集即可，这里的表达量需要log一下。

 1if (!requireNamespace("BiocManager", quietly = TRUE)) 
 2install.packages("BiocManager") 
 3BiocManager::install("SingleR")
 4
 5library(SingleR) 
 6ref_data <- HumanPrimaryCellAtlasData()  #参考数据集
 7ref_data 
 8# ref_data<-readRDS(‘human.hpca.rds’)    #对象 SummarizedExperiment
 9library(scRNAseq) 
10querry <- as.matrix(rds@assays$RNA@data) 
11# SingleR() expects log-counts, but the function will also happily take raw # counts for the test dataset. The reference, however, must have log-values. 
12pred.hesc <- SingleR(test = querry, ref =ref_data, labels = ref_data$label.fine) 
13###这里选择的是用label.fine，而没有用label.main，主要是考虑结果的精细化
14pred.hesc

2. SingleR的结果

SingleR的具体结果如下表格所示，其中labels就是我们最终所需要的结果。

3. SingleR的参考数据集介绍

下面表格为官方提供的7种数据集的相关信息，人有5个数据集、小鼠2个。这些数据集有不同的特征，如果是做免疫相关的研究，则选择免疫相关的数据集，细胞亚群鉴定的结果可能会更细致，比如可以鉴定到CD4 naïve或者effector细胞亚群。

由于在测试的过程中，官方的函数下载参考数据一直不成功，因此直接下载官方参考数据集，并进行处理。内容较多，已上传到了百度网盘（感兴趣的小伙伴，后台回复小编“SingleR”）。该文件夹下有较多的文件，可以直接使用该文件夹下的rds文件，具体数据处理过程也一并提供了，见文件夹下的r脚本。

4. SingleR的结果展示

SingleR工具也提供官方的结果展示，但是其没有细胞亚群信息，因此在我们真正作分析的时候，还是需要具有样品信息和亚群信息的结果才能做下游分析，为了得到具有亚群注释信息的结果，根据如下代码进行结果展示（分析是基于Seurat的rds文件）。

 1annotation_col<-data.frame(clusters=[e-mail]$seurat_clusters,stim=[e-mail]$stim) #
 2rownames(annotation_col)<-rownames([e-mail]) 
 3order_cells<-annotation_col %>% dplyr::mutate(.,barcod=rownames(.)) %>% dplyr::arrange(.,clusters,stim) #
 4pred.hesc$clusters<-annotation_col[as.vector(rownames(pred.hesc)),]$clusters 
 5pred.hesc$stim<-annotation_col[as.vector(rownames(pred.hesc)),]$stim 
 6pred.hesc<-pred.hesc[as.vector(order_cells$barcod),] 
 7pdf('stim_celltype.pdf',w=12,h=8) 
 8p<-plotScoreHeatmap(pred.hesc,annotation_col=annotation_col,cells.order=order(pred.hesc$clusters))
 9print(p)
10dev.off()
11write.table(dplyr :: bind_cols(Barcode=rownames(pred.hesc),as.data.frame(pred.hesc)),‘singleR_celltype.xls，sep="\t",row.names=F) 

这一期我们介绍两个自动化进行单细胞亚群鉴定的工具，下期我们继续。

作者：尧小飞

审稿：童蒙

编辑：angelica

猜你喜欢的文章
► Denovo mutation 分析► 三代组装软件简介► 如何利用NR库快速进行物种鉴定► udp -葡萄糖6-脱氢酶功能缺失突变可引起隐性发展性癫痫性脑病► 如何拆分pacbio的数据