scRNAseq可谓是目前科研界研究细胞异质性的有效手段，正处于如火如荼的阶段。单细胞分析一个很重要目的就是为了确定细胞的类型。说到单细胞分析，大家第一时间想到的肯定是三大R包Seurat、monocle、scater，但是今天我准备给大家介绍一个新的R包metacell，可以用来聚类和注释细胞类型，功能堪比Seurat，但实现方法却很不一样。
接下来就说说metacell包与其他包不同的地方，这个R包聚类不用建立在PCA的基础上，而是将细胞分成若干的元胞(MetaCell)；最重要的是这个包可以系统性地注释细胞类型，先用层次聚类将元胞大体分成若干组，然后可以根据mark基因的lfp(相当于基因在元胞中的富集情况)值进一步细分分组。
话不多说，下面来介绍这个包使用方法，流程如下图所示：
具体步骤
1、初始化环境
library("metacell")
if(!dir.exists("testdb")){
dir.create("testdb/")
}
scdb_init("testdb/", force_reinit=T) #会使用testdb文件夹作为database用来存储生成的中间对象
if(!dir.exists("figs")){
dir.create("figs/")
}
scfigs_init("figs/") #会使用figs文件夹作为默认路径存储生成的图片
2、导入数据
#向database里导入表达矩阵（会生database中生成mat.pbmc.Rda）
mcell_import_scmat_10x("pbmc", base_dir="http://www.wisdom.weizmann.ac.il/~atanay/metac_data/pbmc_8k/") #'pbmc'为矩阵对象的id,可将base_dir换成包含10x的矩阵文件的本地路径。
mat <- scdb_mat("pbmc") #从database获取矩阵对象
#或者从 sparse Matrix of class "dgCMatrix"(行为基因，列为细胞) 导入数据，该方式需要自己添加cell_metadata，matrix格式如下：
AAACCCAAGGAGAGTA AAACGCTTCAGCCCAG AAAGAACAGACGACTG
7SK . . .
A1BG . . .
A1CF . . .
#获取cell_metadata信息
cellid <- colnames(matrix)
cell_metadata <- data.frame(cellid,type='10x',batch_set_id=paste0('test_',cellid),amp_batch_id=paste0('test_',cellid),seq_batch_id=paste0('test_',cellid), spike_count=0,row.names=1)
#创建需要的矩阵对象
mat <- scm_new_matrix(matrix, cellmeta, stat_type = "umi")
scdb_add_mat('pbmc', mat) #将矩阵对象添加到database以便后续使用
3、数据质检及过滤
#UMI分布直方图，会在figs目录下生成一个名为pbmc.total_umi_distr.png文件
mcell_plot_umis_per_cell(dbid)
# 过滤不想要的基因和低UMI的细胞
nms <- c(rownames(mat@mat), rownames(mat@ignore_gmat))
ig_genes <- c(grep("^IGJ", nms, v=T),grep("^IGH",nms,v=T),grep("^IGK", nms, v=T),grep("^IGL", nms, v=T))
bad_genes <- unique(c(grep("^MT-", nms, v=T), grep("^MTMR", nms, v=T), grep("^MTND", nms, v=T),"NEAT1","TMSB4X", "TMSB10",ig_genes))
mcell_mat_ignore_genes(new_mat_id='pbmc', mat_id='pbmc', bad_genes, reverse=F)
mcell_mat_ignore_small_cells('pbmc', 'pbmc', 800)
4、选择特征基因
#选择代表性的基因可以加快计算速度
mcell_add_gene_stat(gstat_id='pbmc', mat_id='pbmc', force=T)
mcell_gset_filter_varmean(gset_id='pbmc_feats', gstat_id='pbmc', T_vm=0.08, force_new=T) #（会在database中生成gset.pbmc_feats.Rda）
mcell_gset_filter_cov(gset_id='pbmc_feats', gstat_id='pbmc', T_tot=100, T_top3=2)
mcell_plot_gstats(gstat_id='pbmc', gset_id='pbmc_feats') #（会在database中生成gstat.pbmc.Rda,在figs中生成三个图片pbmc.szcor.png 、pbmc.top3.png、pbmc.varmin.png）
5、构建平衡的细胞图聚类
mcell_add_cgraph_from_mat_bknn(mat_id='pbmc', gset_id='pbmc_feats', graph_id='pbmc_graph', K=100, dsamp=T) #(生成cgraph.pbmc_graph.Rda)
mcell_coclust_from_graph_resamp(coc_id='pbmc_coclust', graph_id='pbmc_graph', min_mc_size=20, p_resamp=0.75, n_resamp=500)
mcell_mc_from_coclust_balanced(coc_id='pbmc_coclust',mat_id='pbmc', mc_id='pbmc_mc',K=30, min_mc_size=30,alpha=2) #(生成coclust.pbmc_coclust.Rda)
6、去除异常值
mcell_plot_outlier_heatmap(mc_id='pbmc_mc', mat_id = 'pbmc', T_lfc=3) #(会在figs中生成pbmc_mc.outlier.png)
mcell_mc_split_filt(new_mc_id='pbmc_mc_f', mc_id='pbmc_mc', mat_id='pbmc', T_lfc=3, plot_mats=F) #（会在database中生成mc.pbmc_mc_f.Rda）
7、选择metacell的mark基因
mcell_gset_from_mc_markers(gset_id='pbmc_markers', mc_id='pbmc_mc_f') #(会在database中生成gset.pbmc_markers.Rda)
8、metacell的二维可视化展示
mc_colorize_default('pbmc_mc_f')#使用默认的颜色给metacell上色，也可提供自定颜色
mcell_mc2d_force_knn(mc2d_id='pbmc_2dproj',mc_id='pbmc_mc_f', graph_id='pbmc_graph')
tgconfig::set_param("mcell_mc2d_cex",1.5, "metacell") #设置metacell在图中的点大小
tgconfig::set_param("mcell_mc2d_height",1100, "metacell") #设置图片的高度
tgconfig::set_param("mcell_mc2d_width",1100, "metacell") #设置图片的宽度
mcell_mc2d_plot('pbmc_2dproj') #(会在figs中生成pbmc_2dproj.2d_graph_proj.png)
以上就可以将数据划分成若干个MC，后续就是对MC的细胞类型做注释。注释之前需要明白一个概念就是lfp，可以表示基因在MC上的富集程度，通过衡量mark基因的lfp值将不同的MC划分成不同的组。
系统性的细胞类型注释
1、将MC层次聚类
mc_hc <- mcell_mc_hclust_confu(mc_id = 'pbmc_mc_f', graph_id = 'pbmc_graph')
mc_sup <- mcell_mc_hierarchy(mc_id = 'pbmc_mc_f', mc_hc = mc_hc, T_gap = 0.04)
#生成pbmc_mc_f.supmc_confu.png
mcell_mc_plot_hierarchy(mc_id = 'pbmc_mc_f', graph_id = 'pbmc_graph', mc_order = mc_hc$order, sup_mc = mc_sup,width = 1500, height = 2400, min_nmc=2, show_mc_ids = T)
2、通过mark基因确定MC的细胞类型
查看聚类具体的某个分支结果
mc_sup[[4]]
$marks
GAPDH CD27-AS1 CXCR6 CARD16 SELL CTLA4 TNFRSF1B ICOS
1.409033 1.425825 1.478982 1.527237 1.571408 1.609004 1.647381 1.729588
TNFRSF4 FOXP3 RGS1 TIGIT CCR8 PIM2 SAT1 CD27
1.908381 1.955422 2.052494 2.065926 2.076183 2.108668 2.131010 2.170227
TNFRSF18 MIR4632 IL2RA BATF
2.185034 2.188458 2.318691 2.542111
$min_marks
CXCR6 SELL CD82 RGS1 SIRPG ICOS IL32 CARD16
0.7418518 0.7487228 0.7501615 0.9337986 1.0267765 1.0475643 1.0644945 1.0738282
CTLA4 CD74 CCR8 GAPDH CD27-AS1 FOXP3 TIGIT SAT1
1.1351219 1.1766789 1.1937338 1.2025314 1.3105893 1.3563242 1.5938838 1.6433828
IL2RA PIM2 CD27 BATF
1.7104313 1.7890192 2.0234116 2.1977205
$marks_gap
PKM IFI6 IL1R2 ENO1 LDHA SRGN SAT1 TNFRSF1B
1.173707 1.179759 1.208043 1.212875 1.219272 1.235725 1.239677 1.257179
APOBEC3C GAPDH ICOS CTSC MIR4632 DDIT4 TNFRSF4 CCR8
1.298520 1.340427 1.342362 1.389729 1.482165 1.484738 1.512940 1.667040
CXCR6 RGS1 TNFRSF18 BATF
1.860873 1.901821 1.971574 2.075978
$marks_gap_anti
GIMAP7 GIMAP4 TXNIP CD52 S100A10
-2.1707354 -1.6104371 -1.5863388 -1.4286298 -1.1785400
FCMR GIMAP1-GIMAP5 SELL LIME1 ALOX5AP
-1.1444732 -1.1084550 -1.0938374 -1.0360991 -0.9999712
AES LDHB GSTK1 RPL13A CCR7
-0.9894325 -0.9848750 -0.9780960 -0.9670300 -0.9525264
GIMAP5 RPS6 UBXN11 RARRES3 LEF1
-0.9523882 -0.9315339 -0.8861049 -0.8690605 -0.8554596
$mcs
[1] 55 53 52 51 54 56
$x_ord
[1] 6.5
$sup_mcs
[1] 40 38 39 58 55 53 52 51 54 56
#mark基因的lfp阈值默认为 > 1.5
* mcs - the ids of metacells in the subtree (those marked in blue)。该分支中MC的id。
* sup_mcs - the ids of metacells in the subtree and in the sibling subtree (those marked in blue and gray)。该父分支(该分支和兄弟分支)中MC的id。
* x_ord - mean rank of the subtree metacells in the confusion matrix。该分支中MC的id在父分支MC的id中的位置rank值的平均值。
* marks - averaging lfp over mcs, showing the top 20 genes。每个mark基因在MC之间的平均值，取前20个基因。
* min_marks - taking the minimal lfp value per gene over mcs, showing the top 20 genes。所有MC中某个mark基因的lfp最小值，显示20个基因。
* marks_gap - showing the genes with maximal difference between the average lfp over mcs and average lfp over metacells in the sibling subtree。该分支中MC的某个mark基因的lfp平均值与兄弟分支中MC的lfp平均值的差值。
* marks_gap_anti - similar to marks_gap but showing genes with minimal difference (e.g. enriched in the sibling subtree。与marks_gap类似，兄弟分支的某个mark基因在该分支里的MC中的lfp平均值与兄弟分支中MC的lfp平均值的差值。
根据mark基因的聚类结果可以将MC大体分成几个分组，然后结合自己关注的mark基因的lfp值来进一步细分分组：
#查看具体mark基因的lfp值分布的步骤如下：
mc <- scdb_mc('pbmc_mc_f')
lfp <- log2(mc@mc_fp)
plt <- function(gene1, gene2, lfp, colors){
plot(lfp[gene1, ], lfp[gene2, ], pch=21, cex=3, bg=colors, xlab=gene1, ylab=gene2)
text(lfp[gene1, ], lfp[gene2, ], colnames(lfp))
}
#查看SLC4A10、KLRB1两个基因的lfp值在MC中的分布
plt(gene1 = 'SLC4A10', gene2 = 'KLRB1', lfp = lfp, colors = mc@colors)
结果如下图：
3、细胞类型注释
根据聚类的结果，将MC和细胞类型整理成注释文件(supmc_key.txt)，制表符分割格式如下：
supid color name
2 #107538 Tfh
4 #59b56f tumor-Treg
8 #a3ffb9 blood-Treg
17 #b3d4f0 naive-TCF7
27 #6d8eaa naive-CXCR6
34 #cee6b9 CD4-cyto
37 #fdb913 CD8-cyto
40 #a76eaf CD8-dysf
51 #d29bc6 CD8-GZMK
supid-聚类热图中的分支节点ID,即热图左上mark基因后的数字
color-细胞类型的颜色
name-细胞类型的名称
如果需要根据基因的lfp值细分，需要另外整理一个注释文件(gene_key.txt)，制表符分割格式如下：
name gene color T_fold
CD8-dysf TUBA1B #a76eaf 2
CD8-MAIT SLC4A10 #ead37a 1
CD8-dysf-KLRC4 KLRC4 #6b3273 3
naive-TCF7 XIST #b3d4f0 1
#根据注释文件重新聚类并生成热图（如果不重名会覆盖之前的热图）
mc_colorize_sup_hierarchy(mc_id = 'pbmc_mc_f',supmc = mc_sup,supmc_key = 'supmc_key.txt', gene_key='gene_key.txt')
heatmapmcell_mc_plot_hierarchy(mc_id = 'pbmc_mc_f',graph_id = 'pbmc_graph',mc_order = mc_hc$order,sup_mc = mc_sup,width=1200, height=2400, min_nmc=2, show_mc_ids= T)
#最后更新一下2D的聚类图，此时生成的图就会注释好细胞类型
mcell_mc2d_plot('pbmc_2dproj')
最终注释效果：
最后，总结来说一下，这个R采用自己的一套聚类方式，引入一个新的MetaCell元胞概念，将细胞分成若干元胞，然后采用聚类和lfp值方式可以将元胞很细致分成不同的细胞类型。对于细胞包的分析也不失为一个不错的选择。这个R包的作者已经将R封装成各种函数，使用起来挺方便的，但也就是因为封装的很厉害导致需要一定的R基础才能修改输出结果。比如说封装的函数包装的基础画图函数plot，且暴露很少的参数修图，所以画出的图不好看。最后附上该包的文章链接：https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1812-2。各位看官们帮忙点个赞吧！！！
