加速基因體定序分析 Workshop 筆記

linnil1
12 min readFeb 23, 2022

--

封面圖

這次 workshop 只有三小時,所以其實很趕,中後段講很快

這次課程分成三個部分,分別為 Parabricks, Rapids, GTC

課程的網址 https://hackmd.io/@Phoenix1016/SkIdTEgnY

國網

我們在台灣杉三號跑資料

跟生醫相關部分的使用說明
https://man.twcc.ai/@Ldk_QYrOR2yo3m8Cb1549A/rkegDKslF

所有使用說明在這裡(老實說 放在 hackmd 的說明 比他們網站的 UX 好太多了)
https://man.twcc.ai/@TWCC-III-manual/H1bEXeGcu
(國網、台灣杉的註冊地方叫 iservice ,第一次會以為走錯ㄌ)

總之國網對生醫不錯,所以除了有”足夠”的儲存空間(每人10T) 以外,還有 一堆大記憶體節點(可以理解成server) 甚至還有 GPU 節點

https://man.twcc.ai/@Ldk_QYrOR2yo3m8Cb1549A/rkegDKslF 截圖

Job Scheduling system

國網使用 slurm 系統 , 其實蠻像原本的 pbs 的

  • submit: sbatch {xx.sh}
    Submitted batch job 894324
  • cancel: scancle {ID}
  • list: squeue
 JOBID PARTITION NAME USER ST TIME NODES NODELIST(REASON)
894324 ngs2gpu NGS20140 linnil1t PD 0:00 1 (QOSMinMemory)
894325 ngs2gpu NGS20140 linnil1t R 0:05 1 gpn3001

然後 shell script 的 header

#!/usr/bin/sh
#SBATCH -A ACD111008 (計畫名稱)
#SBATCH -J linnil1_test (job name)
#SBATCH -p ngs2gpu (queue name)(listed in usage)
#SBATCH -c 12 (cpu)
#SBATCH — mem=180g (memory)
#SBATCH — gres=gpu:2 (GPU)
#SBATCH -o NGS20140605B-2.log (output)
#SBATCH -e NGS20140605B-2.err.log (error output)
#SBATCH — mail-user=linnil1.886@gmail.com (email)
#SBATCH — mail-type=BEGIN,END (email when job started or ended)
echo 123

Container

雖然台灣杉三號沒有 docker (其實有方法拉),但其實有 singularity 可以使用,所以不用太擔心程式版本、或者裝不起來,跑不起來的問題。

module load libs/singularity/3.7.1
singularity pull docker://quay.io/biocontainers/hisat2:2.2.1 — h1b792b2_3
./hisat2_2.2.1 — h1b792b2_3.sif hisat2

GPU

連 login node 都是 V100,好猛

NVIDIA Clara Parabricks

https://www.nvidia.com/en-us/clara/genomics/

Parabricks 不算是一個全新概念軟體,就只是加速舊有的演算法,裡面重寫像 bwa, haplotypecaller 之類的 code 使得該軟體可以在 GPU 上加速。阿,還有一點是參數應該都幫你預設設好了,所以等下執行時只要給他 sample 跟 reference 的位置就好

目前 Parabrciks 已經把 GATK 相關的工具都加速了,未來可能連 GWAS 也會改寫

以下列出所有能跑 GPU 的程式,其實還蠻多的,大概就是 GATK 的那些 tools

https://docs.nvidia.com/clara/parabricks/3.7.0/_images/pb_tools_v37.png

據說 Parabricks 的 result 跟原本的軟體幾乎一樣,只有 randomness 的差別,然後最近還發布 7.5 hr 就能找到 disease 的 benchmark https://blogs.nvidia.com/blog/2022/01/12/world-record-genome-sequencing-parabricks

實際跑一次

在國網裡面已經安裝,只要使用 module load (ml) 就可以了

module purge
module load biology/parabricks/3.7.0

實際使用方式可以參考這個使用說明
https://docs.nvidia.com/clara/parabricks/3.7.0/Tutorials/GettingTheSampleData.html

follow 這個教學,可以從 fastq 做到 vcf 檔,sample 是一個 2.6G 的 WES 的 pair-end 資料,總共跑的時間

  • Read Mapping: 7mins
  • Variant calling: 10mins
  • QCbyBam: 4.5mins

當天他們下的結論

  • 速度快了很多
  • 但會使用 GPU (gpu instance 比較貴)
  • 總共算起來會省下 1/3

實際 sample:

我拿個 Illumina 30x pair-end read (33Gx2),來實測速度。但是這個課程帳號只有 100G 的 quota,跟剛剛說明的不太一樣,寄信後才發現原來是要自己加國網的 Quota ,計費方式就是 quota * 月 ,而不是 GB*s (超過 1500G 每多 1G/月 收 0.2NT 比 aws 的 1.29NT 還便宜很多)

Quota 介面

總之我嘗試了一遍,跑在國網 V100 * 2 花了 114 分鐘,跟講義的 72mins (A100 * 2) 也差太大,其中

  • Read mapping: 79mins
  • BQSR + Haplotype caller: 34mins
  • VCFqc: 11mins

而雖然他們當天說需要兩個以上的 GPU 才能跑,但我設定了一個 GPU 仍然跑得起來,但是需要 > 90G 的 RAM,所以我用 ngs2gpu 這個有兩個 GPU queue 然後強制指定 GPU 使用一個 --num-gpu=1

  • Read mapping: 125mins
  • Haplotype caller: 44mins
  • VCFqc: 30mins (沒辦法指定一個 GPU)

這樣看起來多一個 GPU 可以快 30%,當然可以用四顆 GPU 的 ngs4gpu 跑。不過我好像會 not permitted ,所以直接看講義數據是 能再快 44% (72mins -> 40mins with 4*A100)

另外我發現,其實 Parabricks 已經幫你把 pipeline 包好了 只要下 pbrun germline 就能從 fastq 直接到 vcf,實測是可以動的,但是如果在 BQSR 加上 knownsites 的參數的時候,就會使用超過 180G 的 RAM 然後 crash 掉,相信未來應該會修好(這版是 3.7.0)

CPU 版的話,我使用 GATK 實測 (用 singularity) 在 32threads, 64Gram 的環境

  • bwa: 3:37
  • MarkDuplicates: 2:37
  • BQSR: 7:28
  • HaplotypeCaller: 20:00

總共 33 小時 42 分,跟講義上的 33hrs 差不多

價格

下圖為國網的價格:

國網價格。 GU = CPU*hour or GPU*hour

我們算算我這次實測的(因為不同 GPU 可能影響蠻大的,所以有實測,但是 CPU 部分就沒有)

  • 32x CPU 版 是 33.7 hours(我實際跑的) = 0.08 * 32 * 33.7 = NT$ 86
  • 2x V100 GPU 是 114mins(我實際跑的) = 10 * 2 * 1.9 = NT$ 38

大概在學術價格下可以省一半

至於使用公有雲計算會差多少,我就相信就計算成本而言會便宜 1/3 價格吧

aws 版本: https://aws.amazon.com/marketplace/pp/prodview-apbngojlskcyq ($4.2/hr(4 GPU) ~= 30/GU = 3 * 科技部計畫的計算價格,GPU 還只有 T4)

逆風一下,還有一個工具叫 Sention 他也是把原本的 GATK(java) 用 c 改寫,也是一樣在 CPU 上跑,但速度上快了很多,我們的 30x 資料只花了 8hr 就跑完了,可以想像是剛剛 GATK 版的 1/4 價格,但這裡一樣我沒考慮 LICENSE 的價格。

Rapids

https://rapids.ai/start.html

看了看概念大略是 scale up,也就是可以讓資料跨 CPU(Dask) 甚至跑在 GPU 上(Rapids),甚至跑在多個 GPU 上 (Rapids+Dask)。讓大資料不是問題 ! (note: 原始的 numpy 只能跑 single thread)

再把常用的 (numpy, pandas, scikit-learn) 包裝成統一格式,當然用法跟原本的很像,我們才容易無痛轉移。注: 剩下的 pytorch 跟 plotly 本來就有支援 GPU

當然,常用的演算法如 dataframe, UMAP, t_SNE, k-mean 都有加速版 https://github.com/rapidsai/cuml#supported-algorithms

Rapids 是 開源的,所以可以自行安裝在想要的地方,所有 issue 也可以第一時間自己修
(parabricks 並沒有開源 只有 90 天試用期)

加速方法其實跟 torch 再算 gradient 的概念很像,就是每一行 code 都沒有真實執行,而是在建 graph, 只有 .compute() 後才一次跑整個 graph,就能在其中加速了,這就是 Lazy execution

https://tutorial.dask.org/_images/grid_search_schedule.gif

我試試 numpy

Testing script: numpy vs dash

差距就很明顯出來了

  • numpy: 113.043s
  • dash: 7.719s

當然用 GPU 更快

Testing cupy speed
  • cupy: 0.598s

有用 GPU 加速方法是使用 cuda ,有 [cuDF](https://github.com/rapidsai/cudf), [cuGraph](https://github.com/rapidsai/cugraph), [cuML](https://github.com/rapidsai/cuml) 之類的 c++ library 可以讓 python 使用,都是在 Rapids 這個 project 下的一份子

實際跑一次

現場示範這個 sample
https://github.com/clara-parabricks/rapids-single-cell-examples/blob/master/notebooks/hlca_lung_gpu_analysis-visualization.ipynb

  • Size: 65662 x 26485
  • Loading : 11s
  • Preprocessing: 5.8s
  • Regress out + filter outlier: 41s
  • Size after filtered: 65462 × 5000
  • Cluster by UMAP(without PCA): 30.5s
  • Cluster by UMAP(with PCA 50 component): 2s

GPU 的 memory 其實沒有吃很多,大概 6.6G 而已,我想說不定連 1080 就可以跑了

速度的詳細比較 作者都寫在 repo 的 README 了 https://github.com/clara-parabricks/rapids-single-cell-examples#acceleration

連錢都幫你算好了 !

GTC

GTC 是 NVIDIA 舉辦的 conference,就是再說 NVIDIA 除了原本的做晶片以外,還往各種方向去延伸,說明 GPU 好棒棒,可以加速很多東西

課程的話就非常簡短的介紹幾個 NVIDIA 跟 生醫的合作關係,還有一些會在 GTC 給 TALK

反正就是
1. NVIDIA 加速 -> 研究可以多測試幾種參數 -> 論文的 quality 上升
2. NVIDIA 加速 -> 打敗其他種方式 -> 穫得 獎項 或者是 合作

在今年的 3/22
報名網址 https://www.nvidia.com/gtc/

結論

我覺得 不管 cost 值不值
測試的時候 當然是越快越好 節省研發成本跟人力精神
至於實際的話,說不定 CPU 運算還是比較便宜的,但那也是要跑在大量 samples 後才會荷包有感

更何況他是給出能節省 1/3 的運算成本(Note 這應該不包含 LICENSE),所以我覺得是應該要學習看看的

寫完了,按下 clap 支持我寫文章

--

--

linnil1

目前做生物資訊與演算法,過去做過 Machine Vision(Deep learning),維護伺服器(k8s, docker),部分IOT(rpi, arduino)