作者：堯小飛
審稿：童蒙
編輯：angelica

引言

隨著單細(xì)胞數(shù)據(jù)的增多，單細(xì)胞亞群的鑒定越來越成為一個單細(xì)胞分析的一個關(guān)鍵的砰蠢、限制性的過程亏栈。今天繼續(xù)分享一下我們究竟如何去“讀取單細(xì)胞這本天書”-單細(xì)胞亞群鑒定岳颇。

單細(xì)胞亞群鑒定--cellassign

下圖為cellassign工作的pipeline氧骤，提供表達(dá)矩陣和細(xì)胞類型及其marker基因，就可以通過機器學(xué)習(xí)預(yù)測細(xì)胞類型几于。

1 cellassign安裝

cellassign（https://github.com/Irrationone/cellassign）是使用輸入的表達(dá)矩陣蕊苗，根據(jù)提供的marker基因，基于TensorFlow的概率模塊進(jìn)行機器學(xué)習(xí)的方法預(yù)測細(xì)胞類型沿彭。

因此安裝的時候需要安裝TensorFlow模塊朽砰，也需要python，此工具最好用annoconda進(jìn)行安裝喉刘。

 1   install.packages("tensorflow") 
 2   library(tensorflow) 
 3   install_tensorflow(extra_packages = "tensorflow-probability")
 4   #Please ensure this installs version 2 of tensorflow. You can check this by calling
 5    tensorflow::tf_config() 
 6   #> TensorFlow v2.1.0 (/usr/local/lib/python3.7/site-packages/tensorflow) 
 7   #> Python v3.7 (~/.virtualenvs/r-reticulate/bin/python) 
 8    sess = tf$Session() 
 9    hello <- tf$constant('Hello, TensorFlow!') 
10    sess$run(hello) 
11    BiocManager::install('cellassign'）
12    或者安裝
13    devtools::install_github("Irrationone/cellassign")
14    ##
15    library(SingleCellExperiment)
16    library(Cellassign)
17    library(Seurat)
18    library(scran)
19    library(dplyr)

2 cellassign使用

cellassign數(shù)據(jù)的準(zhǔn)備：

 1 # 表達(dá)矩陣的準(zhǔn)備
 2 count<-as.matrix(rds@assays$RNA@counts)
 3 # sizeFactors的準(zhǔn)備
 4 cell_anns <- data.frame(Barcode =  names(Idents(rds)),celltype=Idents(rds),samples=rds@meta.data$stim)
 5 rownames(cell_anns) <- names(Idents(rds))
 6 sceset <- SingleCellExperiment(assays = list(counts = as.matrix(count)),colData=cell_anns)
 7 str(sceset)
 8 qclust <- quickCluster(sceset, min.size = 30)
 9 sceset <- computeSumFactors(sceset, sizes = 15, clusters = qclust)
10 sceset <- normalize(sceset)
11 saveRDS(sceset,paste(oudir,paste(pref,'sceset.rds',sep="_"),sep="/"))
12 s <- sizeFactors(sceset)
13 # marker基因的準(zhǔn)備
14 marker_mat <-read.table(markerfile,sep="\t",header=T)
15 marker_gene_list <- list()
16 for (i in marker_mat[,1]){
17    print(unlist(strsplit(as.character(marker_mat[marker_mat[,1]==i,2]),split = ",",fixed=T)))
18    marker_gene_list[[i]]<-  unlist(strsplit(as.character(marker_mat[marker_mat[,1]==i,2]),split = ",",fixed=T))
19}
20 print(marker_list_to_mat(marker_gene_list))
21 marker_mat<-marker_list_to_mat(marker_gene_list)
22 fit <- cellassign(exprs_obj = sceset[as.vector(rownames(marker_mat)),], 
23                  marker_gene_info = marker_mat, 
24                  s = s, 
25                  learning_rate = 1e-2, 
26                  shrinkage = TRUE,
27                  verbose = FALSE)
28

cellassign輸入文件主要是三個瞧柔，一個是表達(dá)矩陣，另外一個是marker基因list睦裳，最后一個是sizeFactors文件造锅。

不過在這里sizeFactors計算直接通過rds文件的信息進(jìn)行計算，因此沒有單獨列出廉邑，可以直接看上述命令即可（在進(jìn)行computeSumFactors時候哥蔚，有時候可能會報錯，可以調(diào)大size參數(shù)即可）蛛蒙。

cellassign根據(jù)marker基因的list進(jìn)行機器學(xué)習(xí)糙箍，預(yù)測細(xì)胞類型，因此此軟件耗時較長牵祟，一般來說倍靡，10000個細(xì)胞的話，差不多24h课舍。

3 結(jié)果展示

其結(jié)果展示也與之前講的singleR一樣塌西，提供了熱圖展示，但是沒有細(xì)胞亞群信息筝尾，因此我們需要把亞群信息添加上去捡需，才更能反映我們數(shù)據(jù)的結(jié)果。

 1#官方展示命令
 2pheatmap::pheatmap(cellprobs(fit)) 
 3
 4cell_cluster_cellprobs<-t(cellprobs(fit))
 5colnames(cell_cluster_cellprobs)<-sceset$Barcode
 6annotation_col<-data.frame(samples=rds@meta.data$stim,Celltype=rds@meta.data$seurat_clusters)
 7   rownames(annotation_col)<-rownames(rds@meta.data)
 8   order_cells<-annotation_col %>% dplyr::mutate(.,barcod=rownames(.)) %>% dplyr::arrange(.,Celltype,samples)
 9   gaps_col<-c() 
10   m<-0
11   for (i in 1:max(as.numeric(annotation_col[,c('Celltype')]))){
12    gaps_col<-c(gaps_col,table(annotation_col[,c('Celltype')])[i]+m)
13    m<-table(annotation_col[,c('Celltype')])[i]+m
14   }
15  #熱圖展示命令
16pheatmap(cell_cluster_cellprobs[,as.vector(order_cells$barcod)],show_rownames=T,show_colnames=F,cluster_cols = F,cluster_rows= F,annotation_col = annotation_col,treeheight_row=0,border=FALSE,gaps_col = gaps_col,main = 'All Single-cell Recognition ',cellheight=16,fontsize=16)

cellassign最終會得到一個細(xì)胞類型的熱圖筹淫，顏色越深代表可能性越高站辉。上圖可以看出不同亞群、不同樣品的細(xì)胞類型預(yù)測的可能性损姜。顏色越紅饰剥，預(yù)測可能性越高，最高為1摧阅。

單細(xì)胞亞群鑒定--clusterprofiler

clusterprofiler是Y叔開發(fā)的目前功能富集最為廣泛的工具汰蓉。

由于有marker基因數(shù)據(jù)庫，因此可以用marker基因數(shù)據(jù)庫對單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)的marker基因list進(jìn)行細(xì)胞類型富集分析棒卷。

根據(jù)富集結(jié)果進(jìn)行細(xì)胞亞群的判斷顾孽，這是做細(xì)胞亞群分析最快的方法，但是此方法有個缺陷比规，那就是如果經(jīng)過sort后的細(xì)胞或者其他純化后的細(xì)胞若厚，關(guān)鍵基因不在marker列表中，其結(jié)果將極其不準(zhǔn)確蜒什。

clusterprofiler：
http://www.bioconductor.org/packages/release/bioc/vignettes/clusterProfiler/inst/doc/clusterProfiler.html

1 clusterprofiler安裝和使用

 1if (!requireNamespace("BiocManager", quietly = TRUE)) 
 2install.packages("BiocManager")
 3BiocManager::install("clusterProfiler") 
 4
 5cell_markers <- vroom::vroom('http://bio-bigdata.hrbmu.edu.cn/CellMarker/download/Human_cell_markers.txt') %>%
 6   tidyr::unite("cellMarker", tissueType, cancerType, cellName, sep=", ") %>% 
 7   dplyr::select(cellMarker, geneID) %>%
 8   dplyr::mutate(geneID = strsplit(geneID, ', '))
 9cell_markers
10markerfile<-"all.markers.csv" #Seurat分析fFindAllMarkers函數(shù)得到所有亞群差異的輸出csv文件
11seurat_marker<-function(markerfile,sep=",",colname="cluster"){
12    marker_list<-list()
13    print(paste("開始讀取markerfile文件测秸，時間為：",Sys.time(),sep=" "))
14    marker<-fread(markerfile,header=T,sep=sep,stringsAsFactors = FALSE,quote = "",data.table = F)
15    for (i in unique(marker[[colname]])){
16        marker_list[[as.character(i)]]<-as.vector(marker[marker[[colname]]==i,]$gene)
17    }
18    #names(marker_list)<-unique(marker[[colname]])
19    print(paste("完成markerfile轉(zhuǎn)換成marker_list，時間為：",Sys.time(),sep=" "))
20    return (marker_list)
21}
22marker_list <-seurat_marker(markerfile)
23y <- compareCluster(marker_list, fun='enricher',TERM2GENE=cell_markers,minGSSize=1,pvalueCutoff = 1, pAdjustMethod = "BH", qvalueCutoff = 1)
24pdf(“compareCluster_cell_markers.pdf, w=12,h=8)
25p1<-dotplot(y, showCategory=as.numeric(showCategory),includeAll=TRUE)
26print(p1)
27dev.off()
28write.table(y@compareClusterResult,paste(outdir,paste(prefix,"compareCluster_cell_markers.xls",sep="_"),sep="/"),quote=F,sep="\t")

2 結(jié)果解讀

這里可以使用enricher或者compareCluster函數(shù)進(jìn)行富集分析灾常，然后可以通過氣泡圖進(jìn)行展示霎冯，最終得到如上所示氣泡圖，每一行為一個細(xì)胞類型岗憋，每一列為一個細(xì)胞亞群肃晚，顏色越紅代表富集越顯著性，氣泡越大仔戈，代表富集程度越高关串。

總結(jié)

目前來看，幾乎所有的工具都是針對人和小鼠兩個模式生物监徘，較少的工具支持其他物種晋修。因此對于除了人、小鼠以外的物種而言凰盔，可能像scMCA/scHCA工具不能適用墓卦，marker基因數(shù)據(jù)庫也有較多的限制，需要從文獻(xiàn)中尋找想要的marker基因户敬。如果有參考數(shù)據(jù)集落剪，可以根據(jù)已有的參考數(shù)據(jù)集通過singleR進(jìn)行細(xì)胞類型預(yù)測或者根據(jù)marker基因通過cellassign和Garnett類似半監(jiān)督的工具進(jìn)行細(xì)胞亞群預(yù)測睁本。

參考文獻(xiàn)

Abdelaal T, Michielsen L, Cats D, et al. A comparison of automatic cell identification methods for single-cell RNA sequencing data[J]. Genome Biology, 2019, 20(1): 1-19.
Han X, Wang R, Zhou Y, et al. Mapping the Mouse Cell Atlas by Microwell-Seq[J]. Cell, 2018, 172(5): 1307-1307.
Han, X. et al. Construction of a human cell landscape at singlecell level. Nature https://doi.org/10.1038/s41586-020-2157-4 (2020).
Park J, Shrestha R, Qiu C, et al. Single-cell transcriptomics of the mouse kidney reveals potential cellular targets of kidney disease[J]. Science, 2018, 360(6390): 758-763.
Zhang X, Lan Y, Xu J, et al. CellMarker: a manually curated resource of cell markers in human and mouse[J]. Nucleic Acids Research, 2019.
Franzen O, Gan L, Bjorkegren J, et al. PanglaoDB: a web server for exploration of mouse and human single-cell RNA sequencing data[J]. Database, 2019.
Aran D, Looney A P, Liu L, et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage.[J]. Nature Immunology, 2019, 20(2): 163-172.
Zhang A W, Oflanagan C H, Chavez E, et al. Probabilistic cell-type assignment of single-cell RNA-seq for tumor microenvironment profiling[J]. Nature Methods, 2019, 16(10): 1007-1015.
Yu G, Wang L G, Han Y, et al. clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters[J]. Omics A Journal of Integrative Biology, 2012, 16(5): 284-287.
Faget J, Groeneveld S, Boivin G, et al. Neutrophils and Snail Orchestrate the Establishment of a Pro-tumor Microenvironment in Lung Cancer[J]. Cell Reports, 2017, 21(11): 3190-3204.

關(guān)注“生信阿拉丁”，第一時間查收“新款”生信學(xué)習(xí)干貨忠怖。

單細(xì)胞亞群鑒定知多少（三）

單細(xì)胞亞群鑒定知多少（三）

引言

單細(xì)胞亞群鑒定--cellassign

1 cellassign安裝

2 cellassign使用

3 結(jié)果展示

單細(xì)胞亞群鑒定--clusterprofiler

1 clusterprofiler安裝和使用

2 結(jié)果解讀

總結(jié)

參考文獻(xiàn)

單細(xì)胞亞群鑒定知多少（三）

引言

單細(xì)胞亞群鑒定--cellassign

1 cellassign安裝

2 cellassign使用

3 結(jié)果展示

單細(xì)胞亞群鑒定--clusterprofiler

1 clusterprofiler安裝和使用

2 結(jié)果解讀

總 結(jié)

參 考 文 獻(xiàn)

總結(jié)

參考文獻(xiàn)