在超算上用sentieon快速完成NGS的变异检测

sentieon

很久之前就听说sentieon在跑calling variants的速度非常快，能甩GATK 不知道多少条街，但是一直缺少一次机会去进行测试。这里感谢sentieon软件公司的张春风提供的测试名额，让我不需要在使用超算测试的时候不用那么漫长的等待就能够测试软件。

这篇教程分为两个部分，第一步部分是sentieon的使用环境，第二部分记录我如何登陆超算服务器

使用sentieon进行变异检测

变异检测我是按照最简单的分析流程来，也就是比对，标记重复，检测变异。下图中的BQSR部分是人类专属，我一个搞植物的就不凑热闹了

流程图

基因组的比对

在基因组比对之前，你需要先建立索引。建立索引不能在你登录节点上运行（会被服务器自动kill掉），你也需要进行任务投递。但是bwa建立索引是单线程，投递任务则会占用一个结点的所有24核，计算核时却是按照24核来计算，所以建议在本地完成。

这里为了演示投递任务的使用方法，所以依旧在超算上投递任务。新建一个文件，明明为make_bwa_index.sh，内容如下

#!/bin/bash# load bwa enviornment
module load bwa/0.7.15-gcc-4.8.5-bwakitWORKDIR=/BIGDATA1/ac_sippe_hzhou_1/SNP/ref
cd $WORKDIR# build index
bwa index ref.fa

然后用chmod修改权限，用yhrun进行任务投递

chmod 755 make_bwa_index.sh
yhrun -n 1 -p free make_bwa_index.sh   &

之后就可以通过yhq查看任务运行情况

建立索引

下一步是用bwa-mem进行基因组比对。新建一个align.sh，复制如下内容

#!/bin/bashmodule load bwa/0.7.15-gcc-4.8.5-bwakit
module load sentieon/201711THREADS=24
WORKDIR=/BIGDATA1/ac_sippe_hzhou_1/SNP
REF=/BIGDATA1/ac_sippe_hzhou_1/SNP/ref/ref.facd $WORKDIRecho "running alignment and sorting"bwa mem -t $THREADS -M -R "@RG\tID:A1\tSM:A1\tPL:illumina" $REF A1_1.fq.gz A1_2.fq.gz 2> log.txt |sentieon util sort --sam2bam --sortblock_thread_count $THREADS -o A1.bam -

增加可执行权限，然后进行运行

chmod 755 align.sh
yhrun -n 1 -p free align.sh   &

当我投递完任务之后，我就会想到一个事情，这个任务到底有没有在运行呢？但是直接在登录节点输入top是看不到运行的进程的，因为你的任务是在计算节点运行，那么我们应该怎么做呢？

经过我的一番摸索，发现不不太难，首先用yhq找到NODELIST(REASON)这一栏下面的节点名，然后用ssh 节点名就可以登录到对应的节点。比如说我投递的bwa目前就是在满负荷运行中。

运行情况

尽管程序在计算节点里计算，但是运行过程中的信息还是会出现在屏幕上，所以理论上，我应该用yhrun -n 1 -p free align.sh 2> align.log来运行程序比较合适。最后的输出信息如下

algo: util-sort
output file size: 14164441164
reads: 238542749
merge_execute: 1 calls 1.438 user 1.102 sys 237.617 real
bam_mem_sort: 103 calls 160.872 user 9.389 sys 175.551 real
bam_write: 103 calls 480.751 user 62.691 sys 602.597 real
execute: 1 calls 1.438 user 1.103 sys 4818.172 real
merge: 1192 calls 3459.932 user 114.072 sys 4960.035 real
parse_chunk: 1639 calls 443.749 user 37.959 sys 521.809 real
read_chunk: 1640 calls 108.445 user 59.480 sys 4476.118 real
sort_block: 1 calls 109.759 user 59.688 sys 4580.548 real
overall: 1783939072 mem 5704.932 user 372.131 sys 4818.194 real

由于每个计算节点只有24个CPU，因此用bwa比对和sentieon的排序速度并不是特别的快，如果你本地的服务器资源很多的话，可以将这一步先跑完，直接上传BAM文件。

标记重复

标记重复这一步可以用sentieon driver提供的Dedup完成，我写了一个markdup.sh脚本用于提交任务

#!/bin/bashmodule load sentieon/201711THREADS=24
WORKDIR=/BIGDATA1/ac_sippe_hzhou_1/SNP
REF=/BIGDATA1/ac_sippe_hzhou_1/SNP/ref/ref.facd $WORKDIRsentieon driver -t $THREADS -i A1.bam --algo LocusCollector --fun score_info A1_score.txtsentieon driver -t $THREADS -i A1.bam --algo Dedup \--rmdup --score_info  A1_score.txt --metrics A1_dedup_metrics.txt A1_deduped.bam

修改运行权限，递交任务。

chmod 755 markdup.sh
time yhrun -n 1 -p free markdup.sh  2> markdup.log  &
real    5m27.257s
user    0m0.036s
sys     0m0.016s

由于使用sentieon的核时费用是原来的三倍，所以我还测试了sambamba，因为这货的速度也非常的快。我也用了24个CPU进行对比，发现时间需要40分钟。所以还是用sentieon处理数据比较省时间。

注：节点为独占式使用方法，每个节点为24核心。每个节点运行一小时的核时是24，考虑到计费标准是X毛/核时，也就是一个结点使用一小时的费用是X*24/10元。

变异检测

最后，就是最重要的一步，也就是原本要一天，现在可能吃顿午饭功夫就能完成的变异检测流程，这一步用到的是sentieon driver中的--algo Haplotyper, 它允许一次提供多个输入文件，我写了一个snp_calling.sh脚本用于提交任务

#!/bin/bashmodule load sentieon/201711THREADS=24
WORKDIR=/BIGDATA1/ac_sippe_hzhou_1/SNP
REF=/BIGDATA1/ac_sippe_hzhou_1/SNP/ref/ref.facd $WORKDIRsentieon driver -t $THREADS --algo Haplotyper -i A1_deduped.bam -i A2_deduped.bam --algo  Haplotyper output.vcf.gz

还是增加运行权限，然后提交任务。

chmod 755 snp_calling.sh
time yhrun -n 1 -p free snp_calling.sh  2> Haplotyper.log  &

这里，我就用到了2个样本，还不能体现出sentieon的性能，最好的测试应该是人类30X重测序，然后给个1000个群体样本，那就是真的可怕了。

使用感想

相对于实际的运行时间，本次分析的真正限速步骤应该是数据上传了。3M/s的水管速度，已经2M/s的下载速度真的很让我抓狂，尤其是在我测试的时候，还看到一个5G评测视频，下载速度可以达到600M/s, 上传可以达到80M/s，真的是没有比较就没有伤害。

当然这次主要是测试Sentieon，这个软件在运行速度上无可挑剔，用起来也很顺手，集成了比对软件(bwa-mem), 排序和去重。

当然了，在命令行参数写的有些小瑕疵。

sentieon driver --help --algo Haplotyper和 sentieon driver --algo Haplotyper --help 居然不一样，前者才能够给出使用说明。
Haplotyper的帮助文档没有提示-i参数能够添加多个样本，我需要先去搜索到sentieon的文档才弄清楚多个样本的用法。

当然任何软件在使用的时候都应该先看帮助文档，如果我先去搜索到https://support.sentieon.com/manual/DNAseq_usage/dnaseq/的话，那么看到他们的流程图和然后按照提供的案例跑一遍的话，感觉压力就会小很多。

当然还有一个问题，普通用户是否会愿意花钱提高速度？

这篇文章因为在测试，所以用的是普通的shell脚本，如果我是批量样本运行的话，我应该写一个snakemake脚本自动控制任务投递才是最优的选择。

下面是介绍我登陆超算的过程，可以不看哦

虚拟专用网络验证

为了保证服务器的安全，我们不能直接通过SSH的方式直接连接服务器，服务器也无法用wget从外网下载数据。

这一步主要参考「虚拟专用网络客户端使用手册」。

由于我是互联网用户，所以需要从如下的下载地址中下载虚拟专用网络软件(这一步需要账号和密码)

https://vpn3.nscc-gz.cn:4433
https://vpn1.nscc-gz.cn:4433

之后打开Hillstone Secure Connect，按照互联网登录的方式填写账号和密码

会话窗口

登录之后，在Windows的右下角有一个小图标出现

登录状态

如果不需要使用的话，右击它可以选择退出

使用Xshell进行登陆。

有很多方式可以登录服务器，这里只介绍Windows的Xshell，因为教育版免费，而且用起来很爽。

第一步是新建一个会话

新建一个会话

第二步是填入IP地址和端口

设置

第三步是修改用户验证方式，这里要导入服务方提供的用户密钥

导入密钥

第三步：在登录的时候输入用户名即可

用户密码

数据上传

如果Xshell安装了Xftp, 那么可以点击工具栏的文件夹图标即可。

启动上传

关于上传速度，上传的速度会慢慢的飚起来

起步速度

最快速度应该可以达到3M/S

我的最快速度

数据从服务器到本地的传输大概在了2M/s左右。

软件使用

天河二号使用module进行软件管理（正好也是我的服务器的管理方式），目前有370个软件（如果包括不同软件版本，则有700多个）

$ module avai > avail
$ tail -n +2 avail | awk '{print $1"\n" $2"\n"$3}'| cut -d '/' -f 1 | sort | uniq | wc -l
370

以我们测试的sentieon为例，加载方式为

module load sentieon/201711

对于基因组组装而言，我目前就看到了Canu 1.8版本，个人觉得够用了。

作业提交

天河二号超算任务投递相关的常用命令如下:

yhi: 查看用户可用节点和节点状态，主要关注输出结果中的STATE, 如果有idle, 则表示有空余资源可以使用
yhq: 查看用户的作业运行情况，主要关注TIME和NODE，前者是运行时间，后者是占用节点数。作者状态用R(运行中)、 PD(paidui,排队)、 S(作业挂起)、 CG(作业运行完毕，等待资源释放) 表示。
yhcancer: 任务取消
yhrun: 任务提交, 核心命令