使用SvABA进行结构变异分析

Posted on by KevinZhou

官方github: https://github.com/walaj/svaba
conda 安装后直接运行即可

构建批量运行脚本：

#!/bin/bash
# ============================================================
# SvABA somatic 全基因组流程 —— 脚本生成命令
# 运行本文件会生成可执行的 svaba.sh
# sample_pairs.txt 格式(tab 分隔): 第1列=tumor($a[0])  第2列=normal($a[1])
# 在有 ../align 的工作目录下运行生成出的 svaba.sh
#
# 三步:
#   1) svaba run        全基因组 call(不加 -k,即不限制染色体)
#   2) svaba postprocess 合并/排序/去重/建索引
#   3) svaba tovcf       转成 VCFv4.5
# 输出前缀 -a 用 tumor 名 $a[0],各样本结果落在以 $a[0] 命名的目录
# ============================================================

REF=/home/zhoukaiwen/database/GRCh38/genecode_GRCh38.p14.genome.fa
BLACKLIST=/home/zhoukaiwen/database/gridss/ENCFF356LFX.bed
THREADS=8

perl -ne 'chomp;next if /^$/;@a=split/\t/;
$ref="/home/zhoukaiwen/database/GRCh38/genecode_GRCh38.p14.genome.fa";
$bl="/home/zhoukaiwen/database/gridss/ENCFF356LFX.bed";
$t=8;
print "mkdir -p $a[0] && cd $a[0] && svaba run -t ../../align/$a[0]_bqsr.bam -n ../../align/$a[1]_bqsr.bam -G $ref -a $a[0] -p $t --blacklist $bl > $a[0].run.log 2>&1 && svaba postprocess -i $a[0] -t $t -m 32G >> $a[0].run.log 2>&1 && svaba tovcf -i $a[0].bps.sorted.dedup.txt.gz -b ../../align/$a[0]_bqsr.bam -a $a[0] >> $a[0].run.log 2>&1 && echo $a[0] svaba ok\n";
' sample_pairs.txt > svaba.sh

生成文件参考：

mkdir -p FETB06-BLPT-E && cd FETB06-BLPT-E && svaba run -t ../../align/FETB06-BLPT-E_bqsr.bam -n ../../align/FETB06-N_bqsr.bam -G /home/zhoukaiwen/database/GRCh38/genecode_GRCh38.p14.genome.fa -a FETB06-BLPT-E -p 8 --blacklist /home/zhoukaiwen/database/gridss/ENCFF356LFX.bed > FETB06-BLPT-E.run.log 2>&1 && svaba postprocess -i FETB06-BLPT-E -t 8 -m 32G >> FETB06-BLPT-E.run.log 2>&1 && svaba tovcf -i FETB06-BLPT-E.bps.sorted.dedup.txt.gz -b ../../align/FETB06-BLPT-E_bqsr.bam -a FETB06-BLPT-E >> FETB06-BLPT-E.run.log 2>&1 && echo FETB06-BLPT-E svaba ok
mkdir -p FETB06-BLPT-M && cd FETB06-BLPT-M && svaba run -t ../../align/FETB06-BLPT-M_bqsr.bam -n ../../align/FETB06-N_bqsr.bam -G /home/zhoukaiwen/database/GRCh38/genecode_GRCh38.p14.genome.fa -a FETB06-BLPT-M -p 8 --blacklist /home/zhoukaiwen/database/gridss/ENCFF356LFX.bed > FETB06-BLPT-M.run.log 2>&1 && svaba postprocess -i FETB06-BLPT-M -t 8 -m 32G >> FETB06-BLPT-M.run.log 2>&1 && svaba tovcf -i FETB06-BLPT-M.bps.sorted.dedup.txt.gz -b ../../align/FETB06-BLPT-M_bqsr.bam -a FETB06-BLPT-M >> FETB06-BLPT-M.run.log 2>&1 && echo FETB06-BLPT-M svaba ok
mkdir -p FETB06-FA-OL-E && cd FETB06-FA-OL-E && svaba run -t ../../align/FETB06-FA-OL-E_bqsr.bam -n ../../align/FETB06-N_bqsr.bam -G /home/zhoukaiwen/database/GRCh38/genecode_GRCh38.p14.genome.fa -a FETB06-FA-OL-E -p 8 --blacklist /home/zhoukaiwen/database/gridss/ENCFF356LFX.bed > FETB06-FA-OL-E.run.log 2>&1 && svaba postprocess -i FETB06-FA-OL-E -t 8 -m 32G >> FETB06-FA-OL-E.run.log 2>&1 && svaba tovcf -i FETB06-FA-OL-E.bps.sorted.dedup.txt.gz -b ../../align/FETB06-FA-OL-E_bqsr.bam -a FETB06-FA-OL-E >> FETB06-FA-OL-E.run.log 2>&1 && echo FETB06-FA-OL-E svaba ok

用于提交sbatch的脚本：

#!/bin/bash
#SBATCH -J svaba
#SBATCH -N 1
#SBATCH -n 1
#SBATCH -c 8
#SBATCH --mem=32G
#SBATCH --array=1-59%30
#SBATCH -o slurm_log/svaba_%A_%a.out
#SBATCH -e slurm_log/svaba_%A_%a.err

# 按行取出 svaba.sh 中的命令(每行一个样本)
CMD=$(sed -n "${SLURM_ARRAY_TASK_ID}p" svaba.sh)

echo "================================================="
echo "JobID      : $SLURM_JOB_ID"
echo "Array Task : $SLURM_ARRAY_TASK_ID"
echo "Command    : $CMD"
echo "Start Time : $(date)"
echo "================================================="

eval "$CMD"

echo "Finished : $(date)"