作者：尧小飞
审稿：童蒙
编辑：angelica

引言

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

1 单细胞亚群鉴定-scMCA/scHCA

单细胞亚群鉴定：scMCA/scHCA(http://bis.zju.edu.cn/MCA/blast.html)

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

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

1.1 scMCA/scHCA的安装

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

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

1.2 scMCA/scHCA的使用

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

1.3 结果展示

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

#shinny展现形式
scHCL_vis(mca_result)

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

 #' @export
 gettissue <- function(x,Num=3){
   top_markers <-order(x,decreasing = T)[1:Num]
   return(top_markers)
 }
 cors <- mca_result$cors_matrix
 cors_index <- apply(cors,2,gettissue,numbers_plot)
 cors_index <- sort(unique(as.integer(cors_index)))
 data = cors[cors_index,]
annotation_col<-data.frame(stim=rds@meta.data$stim,Celltype=rds@meta.data$seurat_clusters)
rownames(annotation_col)<-rownames(rds@meta.data)
order_cells<-annotation_col %>% dplyr::mutate(.,barcod=rownames(.)) %>% dplyr::arrange(.,Celltype,stim)
gaps_col<-c()
m<-0
for (i in 1:max(as.numeric(annotation_col[,c('Celltype')]))){
  gaps_col<-c(gaps_col,table(annotation_col[,c('Celltype')])[i]+m)
  m<-table(annotation_col[,c('Celltype')])[i]+m
}
library(pheatmap)
library(RColorBrewer)
pheatmap(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 = " ")

2 单细胞亚群鉴定-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

2.1 SingleR的安装以及使用

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

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

2.2 SingleR的结果

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

2.3 SingleR的参考数据集介绍

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

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

2.4 SingleR的结果展示

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

annotation_col<-data.frame(clusters=rds@meta.data$seurat_clusters,stim=rds@meta.data$stim) # 
rownames(annotation_col)<-rownames(rds@meta.data)  
order_cells<-annotation_col %>% dplyr::mutate(.,barcod=rownames(.)) %>% dplyr::arrange(.,clusters,stim) # 
pred.hesc$clusters<-annotation_col[as.vector(rownames(pred.hesc)),]$clusters  
pred.hesc$stim<-annotation_col[as.vector(rownames(pred.hesc)),]$stim  pred.hesc<-
pred.hesc[as.vector(order_cells$barcod),]  
pdf('stim_celltype.pdf',w=12,h=8)  
p<-plotScoreHeatmap(pred.hesc,annotation_col=annotation_col,cells.order=order(pred.hesc$clusters)) 
print(p)
dev.off()
write.table(dplyr :: bind_cols(Barcode=rownames(pred.hesc),as.data.frame(pred.hesc)),‘singleR_celltype.xls，sep="\t",row.names=F)

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

关注“生信阿拉丁”，第一时间查收“新款”生信学习干货。

单细胞亚群鉴定知多少（二）