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

如果您只需要快速获取计算结果，不想折腾环境配置的话，可以后台联系我，代为计算beta-NTI。目前笔者使用GPU计算betaNTI的速度极快，核对数据完整性后一般可于十分钟内返回结果，提高您的科研效率。

  beta-NTI（beta nearest taxon index）是微生物生态学经常计算的指标。由于beta-NTI的计算需要反复构建零模型1000次，往往比较耗时，可通过并行的方式减少计算时间。目前R包（iCAMP）利用CPU并行的方式进行计算，并优化了beta-mntd的计算方法，相比stegen等的方法而言大幅减少了计算时间。GPU相比CPU会有更好的并行效果，因而我可以尝试利用GPU并行的方式进一步加快计算速度。

1 先说效果

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

GPU（GTX1050）：1分20秒;
GPU（RTX4090）：16秒;
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 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.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 create -n bnti_gpu python=3.8.16
conda activate bnti_gpu
conda install numpy=1.19.5
conda install cudatoolkit=9.2
conda install cupy=8.3.0

4. 如果第3步安装成功可忽略该部分。 也可导入conda环境的依赖文件bnti_gpu.yaml
   推荐图形界面导入：
   

或者命令行导入

conda env create -n bnti_gpu -f bnti_gpu.yaml

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

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

8. 启动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
R脚本的Rmarkdown文件：beta_nti_gpu_r.Rmd
R脚本的HTML文件：beta_nti_gpu_r.html
 
下载链接： 点击这里进行查看
 
如果没有英伟达显卡或者不愿费劲自行安装conda等相关包，可查看上述链接的代码后，在[autoDL]注册账号，使用分享的公共镜像点我查看公共镜像创建GPU实例，镜像中的运行环境已配置好，可以开机直接使用。

觉得博主写地不错，可以买箱苹果，支持博主。