1 先说效果
18个样本，抽平到8500条序列，4344个OTUs，计算beta-NTI共花费时间如下。如果更好的显卡，更大的数据量，节约的时间应该更加可观。
GPU（GTX1050）：1分20秒
iCAMP包 的bNTIn.p()函数4核并行：约16分钟
2 计算beta-NTI提速的努力
- Stegen等最初在2011年发表论文，利用
beta-NTI推断群落构建，但其原始代码过多的使用for循环，效率很低。 - 尝试利用并行的方法加速
beta-NTI的计算过程，取得一定的效果，但picante包的comdistnt函数计算beta-MNTD的效率也同样很低，所以速度还是偏慢（但比之前Stegen的原始代码好很多了，取决于使用的线程数） iCAMP包对beta-MNTD的函数的算法进行了优化，极大提高了计算的速度，同样也支持使用多线程。尽管如此，计算的速度和普通的beta多样性计算相比还是慢了很多（毕竟要进行1000次零模型的模拟）。- 本文尝试利用
pyhon基于cuda的cupy包调用GPU计算beta-NTI（总体的运行依然是在R中运行的） - 下一步可尝试编写自定义的核函数，更高效地调用GPU的计算能力，将时间降到更低。
3 计算beta-NTI的代码
用python写了计算的beta-NTI的函数，但对OTU table的处理、读取遗传发育树、计算OTUs间的遗传距离等依然是在R环境中进行的，然后在R中调用如下的python函数进行计算即可。
3.1 导入需要的python包
import pandas as pd
import numpy as np
import cupy as cp
3.2 基于cupy的计算beta-MNTD的python代码
def bmntd_gpu(otu,phydist):
otu_p = otu/(otu.sum(1).reshape(otu.shape[0],1))
otu2 = (otu_p != 0)+0
comt = cp.zeros(otu2.shape)
min_d = cp.zeros(otu2.shape[1])
for i in range(0,otu2.shape[0]):
id1 = cp.arange(otu2.shape[1])[otu2[i,:] == 1]
id2 = cp.arange(otu2.shape[1])[otu2[i,:] == 0]
min_d[id1] = 0
min_d[id2] = phydist[id1,:][:,id2].min(0)
comt[i,:] = min_d
D = cp.matmul(comt,cp.transpose(otu_p))
D = (D+cp.transpose(D))/2
return D
3.3 基于cupy的计算beta-NTI的python代码
def bnti_gpu(otu,phydist,N):
N=N
otu_cp = cp.array(otu)
phydist_cp = cp.array(phydist)
row = otu_cp.shape[0]
col = phydist_cp.shape[1]
bmntd_rand_cp = cp.zeros([row,row,N])
for i in cp.arange(N):
id = cp.arange(col)
cp.random.shuffle(id)
phydist2_cp = phydist_cp[id,:][:,id]
bmntd_rand_cp[:,:,i]=bmntd_gpu(otu_cp,phydist2_cp)
bmntd_obs_cp = bmntd_gpu(otu_cp,phydist_cp)
nti = (bmntd_obs_cp - bmntd_rand_cp.mean(2))/bmntd_rand_cp.std(2)
return nti.get()
4 与iCAMP包进行比较
library(iCAMP)
library(reticulate)
4.1 iCAMP包计算用时
data("example.data")
comm <- example.data$comm
pd <- example.data$pd
system.time(bNTIn.p(comm, pd, nworker = 4))
用户 系统 流逝
0.72 0.1143.89
4.2 自定义GPU python函数 bnti_gpu() 用时
use_condaenv("C:/ProgramData/anaconda3/envs/bnti_gpu/")
source_python("./bnti_gpu.py")
system.time(bnti_gpu(comm,pd,as.integer(1000)))
用户 系统 流逝
12.36 3.6716.67
4.3 iCAMP随机模拟2000遍
system.time(bNTIn.p(comm, pd, nworker = 4,rand = 2000))
用户 系统 流逝
0.69 0.0687.03
4.4 bnti_gpu() 随机模拟2000遍
system.time(bnti_gpu(comm,pd,as.integer(2000)))
用户 系统 流逝
24.71 6.6631.65
5 实战操作
- 必要的硬件：
具有navidia的显卡 - 推荐安装anaconda，miniconda也可以。
- 导入conda环境的依赖文件
bnti_gpu.yaml
推荐图形界面导入：
或者命令行导入
conda env create -n bnti_gpu -f bnti_gpu.yaml
- 查看自己显卡对应的cuda版本，目前环境安装的是
cuda 9.2，GTX1050以上显卡应该都支持。可以在anaconda中进行相应的升级或降级，如果不报错不建议修改。 -
参考该博文：如何有效查看电脑显卡对应的CUDA版本 http://t.csdn.cn/GapaH
- 启动
Jupyter Lab，加载beta_nti_gpu_r.ipynb,运行即可。
6 测试数据内容与链接
OTU表：otu table.txt
遗传发育树：tree
conda的依赖文件：bnti_gpu.yaml
调用GPU的python函数：bnti_gpu.py
R脚本的jupyter-notebook文件：beta_nti_gpu_r.ipynb
下载链接： 点击这里进行查看
