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

news2024/11/26 9:45:14

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.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 实战操作

必要的硬件：具有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

下载链接： 点击这里进行查看