sentieon——更具效率的变异检测软件

目前在全基因组范围内做variants calling的软件还是GATK，Samtools为主，我之前的工作基本上属于这两个都用一点，昨天需要跑一批数据的时候发现GATK更新了，很多命令和参数都有了改变，觉得再去看manual就很麻烦，反正都是看不如就学习一下最近非常火的sentieon。
sentieon是一款商业软件（据说价格很贵），相比于目前主流的软件，它的主要特点就是 快、准。
相比于GATK，sentieon用的和GATK一样的数学模型，不一样的是对算法进行了改进，用sentieon的结果和GATK一样，而速度据说能比GATK快20-50倍。
因为我们实验室购买了这个软件，也都全部安装配置到服务器上可以直接调用，下面就贴一下我使用过程中的脚本，十分感谢杨庆勇老师在群里的无私分享。
因为样本比较多，所以我是通过批量提交作业的方法进行的。

#batch_run.lsf

 #BSUB -J sentieon
    #BSUB -n 1
    #BSUB -R span[hosts=1]
    #BSUB -o %J.out
    #BSUB -e %J.err

for sample in ./SRR*_2.fastq.gz
        do
        sh sentieon.sh ${sample} & 
        sleep 1
done
wait

#sentieon.sh

#!/bin/sh
#样本名称
sample=$1
i=$(basename $sample |sed 's/_2.fastq.gz//')
# 加载所需软件
module load sentieon/201808.07
export SENTIEON_LICENSE=mn01:9000
module load SAMtools/1.9
release_dir=/public/home/software/opt/bio/software/Sentieon/201808.07

# fastq文件路径
fq1=${i}_1.fastq.gz
fq2=${i}_2.fastq.gz
# 参考基因组文件
fasta=xxx.fasta
# 输出文件路径
workdir=`pwd`/sentieon_result
# 需要使用的核心数
nt=10
# 比对信息
group_prefix="read_group_name"
platform="ILLUMINA"
mq=30

[ ! -d $workdir ] && mkdir -p $workdir
cd $workdir
# 输出文件
rawCram=$i.cram
sortedCram=$i.q$mq.sorted.cram
depCram=$i.deduped.cram
realnCram=$i.realn.cram
outvcf=$i.sentieon.vcf
exec > $workdir/$i.callVCF.log 2>&1 # call vcf的日志文件
# ******************************************
# 1. 利用 BWA-MEM 进行比对并排序
# ******************************************
bsub -K -J Sentieon -n 10 -o %J.out -e %J.err -R span[hosts=1] "( $release_dir/bin/sentieon bwa mem -M -R \"@RG\tID:${i}\tSM:${i}\tPL:$platform\" \
-t $nt -K 10000000 $fasta $fq1 $fq2 || echo -n 'error' ) | samtools sort -@ $nt  --output-fmt CRAM \
--reference $fasta -o $rawCram - && samtools index -@ $nt $rawCram 
samtools view -hCS -T $fasta -q $mq -o $sortedCram $rawCram && \
samtools index -@ $nt $sortedCram
samtools flagstat $rawCram > $i.stat.raw.txt && \
samtools flagstat $sortedCram > $i.stat.q$mq.txt &

# ******************************************
# 2. Calculate data metrics
# ******************************************
$release_dir/bin/sentieon driver -r $fasta -t $nt -i $sortedCram --algo MeanQualityByCycle ${i}_mq_metrics.txt \
--algo QualDistribution ${i}_qd_metrics.txt --algo GCBias --summary ${i}_gc_summary.txt ${i}_gc_metrics.txt \
--algo AlignmentStat --adapter_seq '' ${i}_aln_metrics.txt --algo InsertSizeMetricAlgo ${i}_is_metrics.txt
$release_dir/bin/sentieon plot metrics -o ${i}_metrics-report.pdf gc=${i}_gc_metrics.txt \
qd=${i}_qd_metrics.txt mq=${i}_mq_metrics.txt isize=${i}_is_metrics.txt
$release_dir/bin/sentieon driver -r $fasta -t $nt -i $sortedCram --algo LocusCollector --fun score_info ${i}_score.txt

# ******************************************
# 3. 去除 Duplicate Reads
# ******************************************
$release_dir/bin/sentieon driver -r $fasta -t $nt -i $sortedCram --algo Dedup --rmdup --cram_write_options version=3.0 \
--score_info ${i}_score.txt --metrics ${i}_dedup_metrics.txt $depCram && rm -f $sortedCram

# ******************************************
# 4. Indel 重排序 (可选)
# 如果只需要最终的比对结果文件，到这里就可以了，这条命令下面的命令都可以注释掉
# ******************************************
$release_dir/bin/sentieon driver -r $fasta -t $nt -i $depCram --algo Realigner --cram_write_options version=3.0 \
$realnCram && rm -f $depCram

# ******************************************
# 5. Variant calling
# ******************************************
$release_dir/bin/sentieon driver -t $nt -r $fasta -i $realnCram --algo Genotyper $outvcf"

需要注意的是在这个地方：

bsub -K -J Sentieon -n 10 -o %J.out -e %J.err -R span[hosts=1] “( $release_dir/bin/sentieon bwa mem -M -R \”@RG\tID:${i}\tSM:${i}\tPL:$platform\”
用了”\”转义了头文件参数的两个引号，否则在执行的时候会因为提交作业命令后的引号（第一个）而引起不能识别-R后面的参数。