加速基因體定序分析 Workshop 筆記
這次 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 節點
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 hisat2GPU
連 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
據說 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 還便宜很多)
總之我嘗試了一遍，跑在國網 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 差不多
價格
下圖為國網的價格:
我們算算我這次實測的(因為不同 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
看了看概念大略是 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
我試試 numpy
差距就很明顯出來了
- numpy: 113.043s
- dash: 7.719s
當然用 GPU 更快
- 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 支持我寫文章
