R语言NMF挖掘基因模块

前言

NMF（非负矩阵分解）目前在生物信息领域主要用于肿瘤分型，目前网上大多数的教程也是针对于对样本聚类的，但其实NMF这个方法还可以对基因进行聚类，已经有许多文章使用了这个方法来识别肿瘤表达的元程序，其实本质也就是挖掘基因模块。

具体原理其实是对样本的聚类分型其实也是基于基因特征来的，因此当样本进行聚类时也会随之产生与各样本簇最相关的一些基因模块，当我们着重关注这些基因模块时，NMF就可以用于挖掘基因模块。

安装

NMF这个R包的安装本身是非常简单的，一句代码即可

install.packages("NMF")   # 安装包的命令
library(NMF) # 加NMF包 

坑在这里，上述代码默认安装最新版的NMF（发文时是0.26版本），但是在某些系统当中，最新版的NMF可能不能够识别全部的CPU，具体可以在library之后看一下反馈

如果红框中的CPU核心数目与实际一致就没什么问题，如果不一致的话（最新版的NMF有可能仅显示2），就需要对NMF的版本进行降级安装（安装0.25版本）

library("devtools")
install_version("NMF", version = "0.25")

NMF的计算由于本身并不是特别复杂，但是运行次数较多的话（也就是后面提到的nrun参数），充分利用CPU的并行运算可以大大提高程序效率，大幅缩短程序运行时间

除此之外，还需检查NMF的依赖情况，首次安装时shared memory支持可能会出现问题

需要检查方框内是否显示OK字样，如果缺少依赖包的话方框内会直接显示缺少的包名，此时直接对应安装缺少的包即可，shared memory共享内存可以在程序运行时大幅节省内存占用，尽可能降低内存不足发生的概率，但是根据官方说明，共享内存不支持Windows ，所以如果你是Windows系统的话，这步可以忽视掉

使用 

#首先对模块数量进行评估，选择最佳的模块数目，即聚类数
#tmp3为要做分解的表达矩阵，行为基因，列为样本
#rank为一个数值向量，代表评估的聚类数目范围，这里我评估了聚成5-10类的情况
#nrun为运行次数，运行次数越多结果鲁棒性也越强，但代价就是程序运行时间变长，对资源的占用也会增加，此时采用多核并行运算可以大幅提高效率，但要注意内存是否溢出，并行核心越多占用的内存也会成倍增长，要结合自己机器的配置来考量
#.option为一些附加参数 +代表启用本参数 -代表关闭本参数 常用参数 v代表verbose打印过程输出(默认关闭) m为共享内存(默认启用) p为并行运算(默认启用，后可跟数值，代表启用多少个核心) 举例vm代表启用verbose启用共享内存 v-m代表启用verbose关闭共享内存 但要注意+在开头可以省略 但中间不可省略 例如v-mp代表启用verbose关闭共享内存和并行运算 v-m+p10则代表启用verbose关闭共享内存开启10核的并行运算
#stop参数代表在运算过程中出现问题是直接停止程序还是跳过，默认是false，即出现问题跳过而不终止程序，这里的跳过指的是rank，例如rank=5时出现问题，跳过继续跑rank=6，这里为了方便排查问题，我的习惯是将stop直接设为true，一旦程序出现问题，立即终止
estimate <- nmf(
    tmp3,
    rank = 5:10,
    nrun = 10,
    .options = "v-m",
    stop = TRUE
  )
#绘制评估结果图表
ggsave(
    paste0(tmp1, "-estimate.pdf"),
    plot(estimate),
    scale = 1,
    width = 12,
    height = 8
  )
#计算最佳rank 官方给出的最佳rank选择方法是找到cophenetic变化最大的前点，这在评估图表上可以观察到，但是不如直接计算来的更方便，尤其是在批量循环中，具体原理是计算后者减去前者的cophenetic值并取绝对值，即是变化值，然后找到这个变化值最大的区间并取前点
#需要注意的是，我这里rank是从5开始的，所以后面是+4，如果你的rank区间起始点不是5的话需要进行相应的调整，否则算出来就不对了，例如rank从2开始的话，后面应该+1
rank = estimate[["measures"]][["cophenetic"]] %>% diff() %>% abs() %>% which.max() + 4
#根据上述计算的最佳rank值进行无监督聚类
estimate2 <- nmf(
    tmp3,
    rank = rank,
    nrun = 10,
    .options = "v-m",
    stop = TRUE
  )
#提取每个基因模块贡献度前50的基因（各模块间的基因可能会有重复）
sig.order <- extractFeatures(estimate2, 50L)
#按照这些基因提取表达谱
exp_re <-
    tmp3[sig.order %>% unlist(),] %>% na.omit() #sig.order有时候会有缺失值
#提取样本的聚类标签
group <- predict(estimate2)
#绘制共识矩阵（跟一致性聚类相似）
pdf(paste0(tmp1, "-consensusmap.pdf"))
  consensusmap(
    estimate2,
    labRow = NA,
    labCol = NA,
    annCol = data.frame("cluster" = group[colnames(exp_re)])
  )
dev.off()
#绘制表达热图
bk <- seq(0, 18, by = 0.01)
tmp4 <- NULL
for (i in 1:rank) {
    tmp4 <- c(tmp4, rep(i, 50))
  }
anno <-
    data.frame(Program = tmp4 %>% as.factor(),
               row.names = rownames(exp_re))
pheatmap(
    exp_re,
    cluster_cols = T,
    cluster_rows = F,
    annotation_row = anno,
    annotation_names_row = F,
    color = colorRampPalette(rev(brewer.pal(
      n = 3, name = "RdYlBu"
    )))(length(bk)),
    breaks = bk,
    show_rownames = F,
    show_colnames = F,
    clustering_method = "ward.D",
    border_color = NA,
    gaps_row = seq(50, dim(exp_re)[1], 50),
    filename = paste0(tmp1, "-heatmap.pdf")
  )
#输出表达矩阵
write.table(
    exp_re,
    file = paste0(tmp1, "-exp.txt") ,
    quote = F,
    sep = "\t"
  )
#输出各模块的基因列表
geneset <- c()
for (i in 1:length(sig.order)) {
    geneset <- cbind(geneset, rownames(tmp3)[sig.order[[i]]])
  }
colnames(geneset) <- paste0(tmp1, "-", 1:rank)
write.table(
    geneset,
    file = paste0(tmp1, "-geneset.txt") ,
    quote = F,
    sep = "\t",
    row.names = F
  )

评估结果图

共识矩阵

表达热图

可能会遇到的坑

1. 程序报错 显示出现NA或者Inf

这是数据问题导致的，如果给的矩阵有某行或者有某些行均为0的话，在做矩阵分解时可能会出现问题，因此在数据方面需要进行质量控制，例如删除全为0的行。

2. 程序报错 显示big.matrix类型超出内存分配上限

目前这个问题疑似是big memory包的bug，NMF的共享内存机制就是基于这个包来实现的，当数据比较大的时候可能就会出现这个问题，我遇到的八千多细胞就已经出现这个问题了，解决办法就是关掉共享内存，在.option参数中设定“-m”即可，但这样的话就会导致高内存占用，但是没办法，总归比报错跑不起来要好。

3. 程序报错 显示 'list' object cannot be coerced to type 'double'

搜遍全网也没找到这个问题出现的原因，更多的可能性还是NMF这个包本身的bug，解决办法是降低nrun数值，即减小运行数目。

4. 绘制consensusmap时报错 显示C堆栈达到上限

当时运行的平台是centos7，目前尚不清楚这个问题在其他平台会不会出现，解决办法是使用R的终端来跑，在打开R session前首先在终端中输入

ulimit -s unlimited #解除C堆栈的限制

然后在终端中输入R 回车 进入R session 把之前的RData load进去然后绘图输出即可

这种方法无法在R的IDE中生效，例如Rstudio，因此只能在终端界面完成这个绘图（可以在绘图输出以后再转回rstudio继续，但这步要放到终端中单独完成）。