利用GPU并行计算beta-NTI，大幅减少群落构建计算时间

1 先说效果

  18个样本，抽平到8500条序列，4344个OTUs，计算beta-NTI共花费时间如下。如果更好的显卡，更大的数据量，节约的时间应该更加可观。

GPU（GTX1050）：1分20秒
iCAMP包 的bNTIn.p() 函数4核并行：约16分钟

2 计算beta-NTI提速的努力

1. Stegen等最初在2011年发表论文，利用beta-NTI推断群落构建，但其原始代码过多的使用for循环，效率很低。
2. 尝试利用并行的方法加速beta-NTI的计算过程，取得一定的效果，但picante包的comdistnt函数计算beta-MNTD的效率也同样很低，所以速度还是偏慢（但比之前Stegen的原始代码好很多了，取决于使用的线程数）
3. iCAMP包对beta-MNTD的函数的算法进行了优化，极大提高了计算的速度，同样也支持使用多线程。尽管如此，计算的速度和普通的beta多样性计算相比还是慢了很多（毕竟要进行1000次零模型的模拟）。
4. 本文尝试利用pyhon基于cuda的cupy包调用GPU计算beta-NTI（总体的运行依然是在R中运行的）
5. 下一步可尝试编写自定义的核函数，更高效地调用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)+0comt = 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] = 0min_d[id2] = phydist[id1,:][:,id2].min(0)comt[i,:] = min_dD = cp.matmul(comt,cp.transpose(otu_p))D = (D+cp.transpose(D))/2return D

3.3 基于cupy的计算beta-NTI的python代码

def bnti_gpu(otu,phydist,N):N=Notu_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.11 43.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.67 16.67

4.3 iCAMP随机模拟2000遍

system.time(bNTIn.p(comm, pd, nworker = 4,rand = 2000))

用户 系统 流逝
0.69 0.06 87.03

4.4 bnti_gpu() 随机模拟2000遍

system.time(bnti_gpu(comm,pd,as.integer(2000)))

用户 系统 流逝
24.71 6.66 31.65

5 实战操作

1. 必要的硬件：具有navidia的显卡
2. 推荐安装anaconda，miniconda也可以。
3. 导入conda环境的依赖文件bnti_gpu.yaml
   推荐图形界面导入：
   

或者命令行导入

conda env create -n bnti_gpu -f bnti_gpu.yaml

4. 查看自己显卡对应的cuda版本，目前环境安装的是cuda 9.2，GTX1050以上显卡应该都支持。可以在anaconda中进行相应的升级或降级，如果不报错不建议修改。

5. 参考该博文：如何有效查看电脑显卡对应的CUDA版本 http://t.csdn.cn/GapaH
   

6. 启动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
 
下载链接： 点击这里进行查看