GATK's Pipeline for calling variants with wgs data

#!/bin/bash

PICARD=/home/lq/storage/call_variants/softwares/picard.jar

gatk=/home/lq/storage/call_variants/softwares/gatk_noDownSample/GenomeAnalysisTK.jar

db=/home/lq/software/data/known_database

reference_file=/home/lq/software/data/reference_sequence/hs37d5.fasta

tmp_dir=/home/lq/Storage/erdai/tmp

result_path=/home/lq/Storage/erdai

##################################################################################################

############################################# QC #######################################

##################################################################################################

echo "####################################### QC #######################################"

in_dir=/home/lq/Storage

out_dir=${result_path}/quality_control

fastqc -t 32 ${in_dir}/HAP1_264_DHG15556-S_1.fq.gz ${in_dir}/HAP1_264_DHG15556-S_2.fq.gz --outdir=${out_dir}

### If necessary:

##################################################################################################

############################################# Cutadapt #######################################

##################################################################################################

echo "####################################### Cutadapt #######################################"

out_dir=${result_path}/1_cutadapt_output

cutadapt -q 5 -o ${out_dir}/1_cutadapt.fastq.gz ${in_dir}/HAP1_264_DHG15556-S_1.fq.gz

cutadapt -q 5 -o ${out_dir}/2_cutadapt.fastq.gz ${in_dir}/HAP1_264_DHG15556-S_2.fq.gz

##################################################################################################

############################################# BWA #######################################

##################################################################################################

echo "####################################### BWA #######################################"

in_dir=${out_dir}

out_dir=${result_path}/2_bwa_output

if [ ! -e ${out_dir}/cutadapt_bwa.sam ]; then

bwa mem -M -P -t 32 -R '@RG\tID:normal\tSM:normal\tLB:normalLib\tPU:runname\tCN:GenePlus\tPL:illumina' $reference_file ${in_dir}/HAP1_264_DHG15556-S_1.fq.gz ${in_dir}/HAP1_264_DHG15556-S_2.fq.gz > ${out_dir}/cutadapt_bwa.sam

####################################################################################################

##################################### Samtools: Compress SAM to BAM #############################

####################################################################################################

echo "########################### Use SAMTOOLS to compress SAM to BAM ###########################"

in_dir=${out_dir}

out_dir=${result_path}/3_samtools_output

if [ ! -e ${out_dir}/cutadapt_bwa.bam ]; then

time samtools view -@ 32 -bS ${in_dir}/cutadapt_bwa.sam > ${out_dir}/cutadapt_bwa.bam

###################################################################################################

########################################### Samtools: Sort #####################################

###################################################################################################

echo "################################## Use SAMTOOLS to sort ###################################"

if [ ! -e ${out_dir}/cutadapt_bwa_sort.bam ]; then

time samtools sort -@ 32 -m 1536M ${out_dir}/cutadapt_bwa.bam -o ${out_dir}/cutadapt_bwa_sort.bam

###################################################################################################

########################################### Samtools: index ####################################

###################################################################################################

echo "################################## Use SAMTOOLS to index ##################################"

if [ ! -e ${out_dir}/cutadapt_bwa_sort.bam.bai ]; then

time samtools index ${out_dir}/cutadapt_bwa_sort.bam ${out_dir}/cutadapt_bwa_sort.bam.bai

###################################################################################################

##################################### Picard: MarkDuplicates ####################################

###################################################################################################

echo "############################# Use Picard to Markduplicates ################################"

in_dir=${out_dir}

out_dir=${result_path}/4_picard_output

if [ ! -e ${out_dir}/cutadapt_bwa_sort_markdup.bam ]; then

time java -d64 -server -XX:+UseParallelGC -XX:ParallelGCThreads=2 -Xms8g -Xmx16g \

-Djava.io.tmpdir=${tmp_dir} -jar $PICARD MarkDuplicates \

I=${in_dir}/cutadapt_bwa_sort.bam O=${out_dir}/cutadapt_bwa_sort_markdup.bam \

METRICS_FILE=${out_dir}/cutadapt_bwa_sort_markdup.bam.metrics ASO=coordinate \

OPTICAL_DUPLICATE_PIXEL_DISTANCE=100 VALIDATION_STRINGENCY=SILENT

###################################################################################################

########################################### Samtools: index ####################################

###################################################################################################

echo "################################## Use SAMTOOLS to index ##################################"

if [ ! -e ${out_dir}/cutadapt_bwa_sort_markdup.bam.bai ]; then

time samtools index ${out_dir}/cutadapt_bwa_sort_markdup.bam ${out_dir}/cutadapt_bwa_sort_markdup.bam.bai

####################################################################################################

########################################### Local Realignment ####################################

####################################################################################################

echo "#################################### Local Realignment ####################################"

in_dir=${out_dir}

out_dir=${result_path}/5_realign_output

if [ ! -e ${out_dir}/cutadapt_bwa_sort_markdup_realign.bam ]; then

time java -d64 -server -XX:+UseParallelGC -XX:ParallelGCThreads=2 -Xms32g -Xmx120g \

-Djava.io.tmpdir=${tmp_dir} -jar $gatk -T RealignerTargetCreator \

-I ${in_dir}/cutadapt_bwa_sort_markdup.bam \

-o ${out_dir}/coRealign.intervals \

-known Mills_and_1000G_gold_standard.indels.b37.vcf \

-R ${reference_file} \

-nt 32

time java -d64 -server -XX:+UseParallelGC -XX:ParallelGCThreads=2 -Xms32g -Xmx120g \

-Djava.io.tmpdir=${tmp_dir} -jar $gatk -T IndelRealigner \

-I ${in_dir}/cutadapt_bwa_sort_markdup.bam \

-targetIntervals ${out_dir}/coRealign.intervals \

-known Mills_and_1000G_gold_standard.indels.b37.vcf.gz \

-R ${reference_file} \

-o ${out_dir}/cutadapt_bwa_sort_markdup_realign.bam \

-nt 32

####################################################################################################

############################################### BQSR #############################################

####################################################################################################

echo "######################################## BQSR #############################################"

in_dir=${out_dir}

out_dir=${result_path}/6_bqsr_output

if [ ! -e ${out_dir}/cutadapt_bwa_sort_markdup_realign_recal.bam ]; then

time java -d64 -server -XX:+UseParallelGC -XX:ParallelGCThreads=2 -Xms30g -Xmx30g \

-Djava.io.tmpdir=${tmp_dir} -jar $gatk -T BaseRecalibrator \

-I ${in_dir}/cutadapt_bwa_sort_markdup_realign.bam \

-o ${out_dir}/cutadapt_bwa_sort_markdup_realign.grp \

-R ${reference_file} \

-knownSites Mills_and_1000G_gold_standard.indels.b37.vcf \

-knownSites dbsnp_138.b37.vcf \

-knownSites b37_cosmic_v73_061615.vcf \

-nct 32

time java -d64 -server -XX:+UseParallelGC -XX:ParallelGCThreads=2 -XX:+UseParallelOldGC -Xms30g -Xmx30g \

-Djava.io.tmpdir=${tmp_dir} -jar $gatk -T PrintReads \

-I ${in_dir}/cutadapt_bwa_sort_markdup_realign.bam \

-BQSR ${out_dir}/cutadapt_bwa_sort_markdup_realign.grp \

-o ${out_dir}/cutadapt_bwa_sort_markdup_realign_recal.bam \

-R ${reference_file} \

-nct 32

####################################################################################################

############################################ GATK: HaplotypeCaller ####################################

###################################################################################

echo "################################## Use GATK HaplotypeCaller ##################################"

in_dir=${out_dir}

out_dir=${result_path}/7_haplotypeCaller_output

time java -d64 -server -XX:+UseParallelGC -XX:ParallelGCThreads=2 -Djava.io.tmpdir=${tmp_dir} \

-jar $gatk -T HaplotypeCaller \

-nct 32 \

-R ${reference_file} \

-ploidy 1 \

-I ${in_dir}/cutadapt_bwa_sort_markdup_realign_recal.bam \

--dbsnp dbsnp_138.b37.vcf \

-o ${out_dir}/cutadapt_bwa_sort_markdup_realign_recal.raw.snps.indels.vcf

####################################################################################################

#######################################VQSR: Build the SNP recalibration model############################

####################################################################################################

echo "########################VQSR: Build the SNP recalibration model######################################"

in_dir=${out_dir}

out_dir=${result_path}/8_vqsr_output

java -jar GenomeAnalysisTK.jar \

-T VariantRecalibrator \

-R ${reference_file} \

-nt 20 \

-input ${in_dir}/cutadapt_bwa_sort_markdup_realign_recal.raw.snps.indels.vcf \

-resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap_3.3.b37.vcf \

-resource:omni,known=false,training=true,truth=true,prior=12.0 1000G_omni2.5.b37.vcf \

-resource:1000G,known=false,training=true,truth=false,prior=10.0 1000G_phase1.snps.high_confidence.b37.vcf \

-resource:dbsnp,known=true,training=false,truth=false,prior=2.0 dbsnp_138.b37.vcf \

-an QD \

-an MQ \

-an MQRankSum \

-an ReadPosRankSum \

-an FS \

-an SOR \

-an DP \

-mode SNP \

-tranche 100.0 -tranche 99.9 -tranche 99.0 -tranche 90.0 \

-recalFile ${out_dir}/recalibrate_SNP.recal \

-tranchesFile ${out_dir}/recalibrate_SNP.tranches \

-rscriptFile ${out_dir}/recalibrate_SNP_plots.R

####################################################################################################

######################### Apply the desired level of recalibration to SNPs in the call set#######################

####################################################################################################

echo "################### Apply the desired level of recalibration to SNPs in the call set########################"

java -jar GenomeAnalysisTK.jar \

-T ApplyRecalibration \

-ef \

-nt 20 \

-R ${reference_file} \

-input ${in_dir}/cutadapt_bwa_sort_markdup_realign_recal.raw.snps.indels.vcf \

-mode SNP \

--ts_filter_level 99.0 \

-recalFile ${out_dir}/recalibrate_SNP.recal \

-tranchesFile ${out_dir}/recalibrate_SNP.tranches \

-o ${out_dir}/exfilter_recalibrated_snps_raw_indels.vcf

####################################################################################################

#######################################VQSR: Build the INDEL recalibration model############################

####################################################################################################

echo "########################VQSR: Build the INDEL recalibration model######################################"

java -jar GenomeAnalysisTK.jar \

-T VariantRecalibrator \

-nt 20 \

-R ${reference_file} \

-input ${out_dir}/exfilter_recalibrated_snps_raw_indels.vcf \

-resource:mills,known=false,training=true,truth=true,prior=12.0 Mills_and_1000G_gold_standard.indels.b37.vcf \

-resource:dbsnp,known=true,training=false,truth=false,prior=2.0 dbsnp_138.b37.vcf \

-an QD \

-an DP \

-an FS \

-an SOR \

-an MQRankSum \

-an ReadPosRankSum \

-mode INDEL \

-tranche 100.0 -tranche 99.9 -tranche 99.0 -tranche 90.0 \

--maxGaussians 4 \

-recalFile ${out_dir}/recalibrate_INDEL.recal \

-tranchesFile ${out_dir}/recalibrate_INDEL.tranches \

-rscriptFile ${out_dir}/recalibrate_INDEL_plots.R

####################################################################################################

######################### Apply the desired level of recalibration to INDELs in the call set#######################

####################################################################################################

echo "################### Apply the desired level of recalibration to INDELs in the call set########################"

java -jar GenomeAnalysisTK.jar \

-T ApplyRecalibration \

-ef \

-nt 20 \

-R ${reference_file} \

-input ${out_dir}/exfilter_recalibrated_snps_raw_indels.vcf \

-mode INDEL \

--ts_filter_level 99.0 \

-recalFile ${out_dir}/recalibrate_INDEL.recal \

-tranchesFile ${out_dir}/recalibrate_INDEL.tranches \

-o ${out_dir}/exfilter_recalibrated_variants.vcf

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 213,254评论 6赞 492
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 90,875评论 3赞 387
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 158,682评论 0赞 348
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 56,896评论 1赞 285
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 66,015评论 6赞 385
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 50,152评论 1赞 291
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 39,208评论 3赞 412
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 37,962评论 0赞 268
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 44,388评论 1赞 304
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 36,700评论 2赞 327
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 38,867评论 1赞 341
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 34,551评论 4赞 335
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 40,186评论 3赞 317
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 30,901评论 0赞 21
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 32,142评论 1赞 267
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 46,689评论 2赞 362
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 43,757评论 2赞 351

GATK's Pipeline for calling variants with wgs data

推荐阅读更多精彩内容