Posts /

R的那些包（2）——NMF

Back to posts /
Twitter Facebook Google+
11 May 2017

本文为NMF包说明文档。主要参考NMF-vignette和NMF参考手册 {NMF}特性如下：

  • 11个内置NMF算法；
  • 4中设置种子(初始化)的算法；
  • 所有算法和种子算法的接口良好的集合在了一起；
  • 提供了一个普适测试框架，用以比较和魔改不同的NMF算法；
  • 可以嵌入自己的算法和设置中资的方法；
  • 提供了绘图及可视化的解释方法；
  • 针对并行运算进行了优化；
  • 植入C++标准算法，对内存cost进行了优化；
  • Optional layer for bioinformatics using BioConductor (Gentleman et al. 2004a);

1 概论

1.1 安装载入NMF

# 在本文的例子中需要用到一个数据集，如果要使用这个数据集，需首先安装Biobase再安装NMF(或在安装完Biobase后重新安装NMF)，安装方法如下：

# 要使用Golub数据集，需要先安装并载入Biobase
source("https://bioconductor.org/biocLite.R")
biocLite("Biobase")

# 如果出现提示更新包的版本，按自己需要选择
Update all/some/none? [a/s/n]:

# 安装NMF
install.packages("NMF")

# 载入NMF
library(NMF)

1.2 NMF包中的算法

# nmfAlgorithm
> nmfAlgorithm()
 [1] "brunet"    "KL"        "lee"       "Frobenius" "offset"   
 [6] "nsNMF"     "ls-nmf"    "pe-nmf"    "siNMF"     "snmf/r"   
[11] "snmf/l"   

# 测试获取其中某种算法参数  
> nmfAlgorithm("lee")
<object of class: NMFStrategyIterative>
 name: lee [NMF]
 objective: 'euclidean'
 model: NMFstd
 <Iterative schema>
  onInit: none
  Update: function (i, v, x, rescale = TRUE, copy = FALSE, eps = 10^-9,
            weight = NULL, ...)
  Stop: 'connectivity'
  onReturn: none

1.3 设置种子

Key Description
ica Independent Component Analysis (ICA) (from the {fastICA} (Marchini et al. 2013)). （只有正值的部分被用来初始化向量）
nnsvd Nonnegative Double Singular Value Decomposition. The basic algorithm contains no randomization and is based on two SVD processes, one approximating the data matrix, the other approximating positive sections of the resulting partial SVD factors utilizing an algebraic property of unit rank matrices. It is well suited to initialize NMF algorithms with sparse factors. Simple practical variants of the algorithm allows to generate dense factors. Reference: (Boutsidis et al. 2008)
none 使用none，用户可以手动调节初始值
random 将矩阵初始化为均匀分布的值，值域范围在[0,max(V)]

1.4 核心接口nmf函数

2 使用实例

这一节里使用Golub数据集展示使用NMF包进行非负矩阵分解的实例。 Golub数据集在许多NMF研究的论文中都被使用过，在{NMF}中Golub数据的存储在一个自定义数据结构ExpressionSet中。

2.1 使用内置模型

data(esGolub) # 载入数据
esGolub<-esGolub[1:200,]  # 取其中前两百个基因作为测试样本（走个流程节省时间）

> res<-nmf(esGolub,3,method="ns",seed="random",theta=0.7, nrun=5)
> res # 查看res
> algorithm(res) # 查看res模型的算法
> fit(res)  # 可以用函数fit查看模型的参数
> View(fitted(res))


> W <- basis(res)
> H <- coef(res) # 得到分解矩阵W和H

# 使用summary可以查看res的参数，如残差、计算时间等
# 同时还可以指定target和class获取更多信息。
# 需注意这里的Cell是esGolub数据集中的一个元素
> summary(res)
> summary(res,target=esGolub)
> summary(res,class=esGolub$Cell)

# 有趣的是，NMF对象就像matrix，data.frame一样，可以随意取行/列的子集
# 这样得到的结果也是一个NMFfit对象
> res[1:10,]
> res[,1:10]
> res[2:33,2:33]

问题：这里metagene-specific features没有看懂

2.2 使用自己的模型

# 直接调用关键字
> nmf(esGolub,3,seed="nndsvd")

# 使用数字种子（random）
> nmf(esGolub,3,seed=123456)

# 使用自己的初始化方法（直接传入nmfModel类的初始W/H矩阵）
> init <- nmfModel(3, esGolub, W=0.5, H=0.3)
> res <- nmf(esGolub, 3, seed=init)
# 这里查看W和H与2.1中一样使用`basis`和`coef`

2.3 并行计算

# 申请4个运算单元进行并行运算，如出现问题将报错
nmf(esGolub,3,nrun=5,.opt="vP4")
# 可以使用`-p`或.pdackkend=NA来强制顺序运算
nmf(esGolub,3,nrun=5,.opt="v-p",seed=123)
nmf(esGolub,3,nrun=5,.opt="v",.pbackend=NA,seed=123)

NMF也可以进行集群运算（HPC cluster），但是笔者条件有限，没有进行测试，有需要的朋友可以参阅[1]中2.5.4节。

2.4 估计向量化的秩

# 从2到6每一个r值都会跑十遍
> estim.r<-nmf(esGolub,2:6,nrun=10)

# 观察不同r值的评测指标
> plot(estim.r)

# 观察不同rank值的一致性矩阵，可以加入其他的注解信息
> consensusmap(estim.r, annCol=esGolub,labCol=NA,labRow=NA)

# 生成与esGolub同维度的随机矩阵做矩阵分解，并与真实数据NMF结果进行比较
> V.random <- randomize(esGolub)
> estim.r.random<-nmf(V.random,2:6,nrun=10)
> plot(estim.r, estim.r.random)

2.5

3 扩展包

Reference

https://www.bioconductor.org/packages/release/bioc/html/Biobase.html


Twitter Facebook Google+