单细胞亚群鉴定知多少（二）
作者：尧小飞
审稿：童蒙
编辑：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)
这一期我们介绍两个自动化进行单细胞亚群鉴定的工具，下期我们继续。
关注“生信阿拉丁”，第一时间查收“新款”生信学习干货。
