PRIMUS续集：如何使用BIC确定 k 值 by 单细胞天地

接上篇

贝叶斯信息准则

在上篇提到，PRIMUS 的k值很关键，但还没有指出如何确定k值。作者提及Bayesian information criterion (BIC) 是一个统计学概念。在统计学中，贝叶斯信息准则（BIC）是在有限的模型集合中选择模型的准则；通常首选BIC较低的型号。它部分基于似然函数，与Akaike信息准则（AIC）密切相关。

BIC=k*ln(n)-2 ln(L)

注：L是在该模型下的最大似然，n是数据数量，k是模型的变量个数。

那么如何实现PRIMUS BIC分析？

BIC 代码实现

我们回到前面的代码：

读入数据

rm(list= ls())
library(SingleCellExperiment)
library(scuttle)
library(data.table)
library(dplyr)
rm(list = ls())
## 数据准备
counts = fread('./dataset4/myData_pancreatic_5batches.txt.gz') %>% as.data.frame()
rownames(counts) <- counts$V1; counts <- counts[,2:ncol(counts)]
dim(counts)
head(counts)[1:5,1:2]
meta = fread('./dataset4/mySample_pancreatic_5batches.txt.gz',) %>% as.data.frame()
rownames(meta) <- meta$V1; meta <- meta[,2:ncol(meta)]
head(meta)[1:5,1:3]
colnames(meta)
dim(meta)

抽样（考虑后续循环次数较多，减少细胞量）

table(meta$celltype, meta$batch) 
#>                  1    2    3    4    5
#>   acinar        958  219  185    6    0
#>   alpha        2326  812  886  190  886
#>   beta         2525  448  270  111  472
#>   delta         601  193  114    9   49
#>   ductal       1077  245  386   96    0
#>   endothelial   252   21   16    0    0
#>   epsilon        18    3    7    0    0
#>   gamma         255  101  197   18   85
#>   macrophage     55    0    0    0    0
#>   mast           25    0    7    0    0
#>   mesenchymal     0   80    0   27    0
#>   MHC class II    0    0    5    0    0
#>   schwann        13    0    0    0    0
#>   stellate      457    0   54    0    0
#>   t_cell          7    0    0    0    0
allCells=rownames(meta)
allbatch = levels(factor(meta$batch))
choose_Cells = unlist(lapply(allbatch, function(x){
  cgCells = allCells[meta$batch == x]
  cg=sample(cgCells,100,replace = F)  # 减少至100 个
  cg
}))

Ct <- counts[,choose_Cells]
mt <- meta[choose_Cells,]

table(mt$celltype, mt$batch)             
#>                1  2  3  4  5
#>   acinar      13 11  9  1  0
#>   alpha       25 40 36 36 53
#>   beta        27 23 16 24 38
#>   delta        6 11  9  3  3
#>   ductal      14 11 18 22  0
#>   endothelial  3  0  1  0  0
#>   gamma        3  3  9  5  6
#>   mesenchymal  0  1  0  9  0
#>   stellate     9  0  2  0  0

构建SingleCellExperiment对象

simData = logNormCounts(SingleCellExperiment(assays = list(counts = Ct )))

simData
library(prism) # 使用作者推荐的方式prism::prism.gain 获取 scaling factor
Ct <- as.matrix(Ct)
dim(Ct)
mt$sizeFactor  <- prism::prism.gain(X = Ct,method = 'poi',alpha = 1 )$v

PRIMUS 运算

library(PRIMUS)

labels2weights <- function(L, levels = sort(unique(L)))
  return (outer(L, levels, `==`) + 0L)

D = t(labels2weights(mt$batch))

library(tidyverse)
### 使用循环，设置k值1-20， 每个循环5次
fits <- lapply(seq.int(5),function(x)
    {
    set.seed(x*1234)
    lapply(seq.int(20),function(x){
    runPrimus(Y = Ct, D = D, g = mt$sizeFactor,
            k = x, max.iter = 100)
})}
)

BIC 运算

### 求BIC值
BIC <- lapply(fits,function(x){
  PRIMUS::calcBIC(Y = Ct,g = mt$sizeFactor,fits = x)
  })
### 数据处理，为可视化做准备
BIC <- as.data.frame(BIC)
colnames(BIC) <- seq.int(5)
BIC <- cbind(k=rownames(BIC),BIC)
df <- reshape::melt(BIC)

可视化

library(ggplot2)
df$k<- factor(df$k,levels = seq.int(20))
ggplot(df,aes(x=k, y= value))+
  geom_boxplot()+theme_bw()+
  geom_vline(aes(xintercept=14),colour='#AA0000')

结束语

至此，完成了PRIMUS的整合的全部流程及k值的选择过程。在本例中，由于fits 值单程计算费时，故在k =1~20 之后并未再进一步运算，后续进一步扩大k值范围，可能能够找到最低点。考虑到该份数据来自于真实数据，各个对应的细胞类型是否能再次细分暂未可知，后续根据BIC值计算的最佳k值>> 实际细胞类别数与此又是否有关联。不求甚解，到此为止。