Homec4science
Diffusion eSCT pipeline interoperability (master) Browse

R2915/example_pipelines/MelArray/pipeline/subworkflows/HaploCaller.snake
88e9313e0804master

Last Change
Raw File

HaploCaller.snake
View Options

##############################################
## MelArray - Haplotype caller sub-pipeline ##
##############################################
#
# @author: Phil Cheng, USZ (for original pipeline)
# @author: Balazs Laurenczy, ETHZ (for Snakemake pipeline)
# @date: 2017-2018
########
# INIT #
########
# if this Snakefile runs as a subworkflow
if "subw_outputs_dict" in locals() and "IDs" in locals():
# specify the final required outputs by expanding the final output with some existing IDs
subw_outputs_dict['HaploCaller'] = expand(paths['recal_indel_vcf'], md5 = IDs_MD5)
#########
# RULES #
#########
# haplotype caller step
rule haplotype_caller:
input:
bqsr_bam = paths['bqsr_bam'],
bqsr_bai = paths['bqsr_bai']
output:
hc_vcf = paths['hc_vcf'],
log:
hc_vcf_out = paths['hc_vcf_out']
shell:
# wrap the whole command in one string
(
# append the Singularity container call (if needed)
SINGULARITY_CMD +
# create the command for GATK with the amount of memory specified in the config
'gatk --java-options "-Xmx{memory}G" '.format(**soft['gatk']) +
# specify the GATK JAR path, the GATK command, the reference genome and some more references
'HaplotypeCaller -R {ref_genome} --dbsnp {dbsnp} '.format(**rule_paths) +
# specify the inputs and outputs
'-I {bqsr_bam} -O {hc_vcf} '.format(**rule_paths) +
# dump all standard and error output to the log file
'&> {hc_vcf_out}'.format(**paths)
# replace all instances of '{id}' with '{wildcards.id}'
).replace('{id}', '{wildcards.id}')
# create the *.g.vcf file list based on the IDs
rule create_gvcf_list:
input:
# use the output of all the haplotype caller steps (reduction / convergence step)
hc_vcf = expand(paths['hc_vcf'], id = IDs)
output:
# use the MD5-hash of the sorted list of IDs as name for the file, to make it unique. This
# way, every time a new ID is added, this list gets re-calculated
vcf_list = paths['vcf_list']
shell:
# wrap the whole command in one string
(
# append the Singularity container call (if needed)
SINGULARITY_CMD +
# get the list of all VCF files names
'find {hc_vcf_root} -name "*{hc_vcf_pattern}" '.format(**rule_paths) +
# filter for .g.vcf files and dump all the file names in the list file
'> ' + paths['vcf_list']
# replace all instances of '{md5}' with '{wildcards.md5}'
).replace('{md5}', '{wildcards.md5}')
# genotype gvcfs step
rule genotype_gvcf:
input:
vcf_list = paths['vcf_list']
output:
hc_geno_vcf = paths['hc_geno_vcf'],
log:
hc_geno_vcf_out = paths['hc_geno_vcf_out']
shell:
# wrap the whole command in one string
(
# append the Singularity container call (if needed)
SINGULARITY_CMD +
# create the command for GATK with the amount of memory specified in the config
'gatk --java-options "-Xmx{memory}G" '.format(**soft['gatk']) +
# specify the GATK JAR path, the GATK command, the reference genome and some more references
'GenotypeGVCFs -R {ref_genome} --dbsnp {dbsnp} '.format(**rule_paths) +
# specify the inputs and outputs
'-V {vcf_list} -O {hc_geno_vcf} '.format(**rule_paths) +
# dump all standard and error output to the log file
'&> {hc_geno_vcf_out}'.format(**paths)
# replace all instances of '{md5}' with '{wildcards.md5}'
).replace('{md5}', '{wildcards.md5}')
# SNV variant recalibrator step
rule var_snv_recal:
input:
hc_geno_vcf = paths['hc_geno_vcf']
output:
var_snv_recal = paths['var_snv_recal'],
var_snv_recal_tranches = paths['var_snv_recal_tranches']
log:
var_snv_recal_out = paths['var_snv_recal_out']
shell:
# wrap the whole command in one string
(
# append the Singularity container call (if needed)
SINGULARITY_CMD +
# create the command for GATK with the amount of memory specified in the config
'gatk --java-options "-Xmx{memory}G '.format(**soft['gatk']) +
# specify the GATK JAR path, the GATK command and the reference genome
'VariantRecalibrator -R {ref_genome} '.format(**rule_paths) +
# specify the input
'-V {hc_geno_vcf} '.format(**rule_paths) +
# specify the outputs
'--recal-file {var_snv_recal} --tranches-file {var_snv_recal_tranches} '.format(**rule_paths) +
# specify all resources / references
# The locations of a space and a colon reversed between GATK 3 and 4. . First try below.
'-resource {hapmap}: '.format(**soft['var_snv_recal']) + # parameter string
'{hapmap} '.format(**rule_paths) + # path of the reference file
'-resource {omni}: '.format(**soft['var_snv_recal']) + # parameter string
'{omni} '.format(**rule_paths) + # path of the reference file
'-resource {1000G}: '.format(**soft['var_snv_recal']) + # parameter string
'{1000G} '.format(**rule_paths) + # path of the reference file
'-resource {dbsnp}: '.format(**soft['var_snv_recal']) + # parameter string
'{dbsnp} '.format(**rule_paths) + # path of the reference file
# add some parameters from the configuration that are specific for this step
'--mode {mode} '.format(**soft['var_snv_recal']) +
# transform the list of annotations in a '-an XXXX -an XXXX -an XXXX' list
'-an {} '.format(' -an '.join(soft['var_snv_recal']['an_list'])) +
# dump all standard and error output to the log file
'&> {var_snv_recal_out}'.format(**paths)
# replace all instances of '{md5}' with '{wildcards.md5}'
).replace('{md5}', '{wildcards.md5}')
# apply SNV variant recalibration step
rule apply_snv_recal:
input:
var_snv_recal = paths['var_snv_recal'],
var_snv_recal_tranches = paths['var_snv_recal_tranches'],
hc_geno_vcf = paths['hc_geno_vcf']
output:
recal_snv_vcf = paths['recal_snv_vcf']
log:
recal_snv_vcf_out = paths['recal_snv_vcf_out']
shell:
# wrap the whole command in one string
(
# append the Singularity container call (if needed)
SINGULARITY_CMD +
# create the command for GATK with the amount of memory specified in the config
'gatk --java-options "-Xmx{memory}G '.format(**soft['gatk']) +
# specify the GATK JAR path, the GATK command and the reference genome
'ApplyVQSR -R {ref_genome} '.format(**rule_paths) +
# specify the inputs
'-V {hc_geno_vcf} --tranches-file {var_snv_recal_tranches} '.format(**rule_paths) +
'--recal-file {var_snv_recal} '.format(**rule_paths) +
# specify the output
'-O {recal_snv_vcf} '.format(**rule_paths) +
# add some parameters from the configuration that are specific for this step
'--ts_filter_level {ts_filter_level} -mode {mode} '.format(**soft['apply_snv_recal']) +
# dump all standard and error output to the log file
'&> {recal_snv_vcf_out}'.format(**paths)
# replace all instances of '{md5}' with '{wildcards.md5}'
).replace('{md5}', '{wildcards.md5}')
# variant indel recalibrator step
rule var_indel_recal:
input:
recal_snv_vcf = paths['recal_snv_vcf']
output:
var_indel_recal = paths['var_indel_recal'],
var_indel_recal_tranches = paths['var_indel_recal_tranches']
log:
var_indel_recal_out = paths['var_indel_recal_out']
shell:
# wrap the whole command in one string
(
# append the Singularity container call (if needed)
SINGULARITY_CMD +
# create the command for GATK with the amount of memory specified in the config
'gatk --java-options "-Xmx{memory}G '.format(**soft['gatk']) +
# specify the GATK JAR path, the GATK command and the reference genome
'VariantRecalibrator -R {ref_genome} '.format(**rule_paths) +
# specify the input
'-V {recal_snv_vcf} '.format(**rule_paths) +
# specify the outputs
'--recal-file {var_indel_recal} -tranches-file {var_indel_recal_tranches} '.format(**rule_paths) +
# specify all resources / references
'-resource {mills}: '.format(**soft['var_indel_recal']) + # parameter string
'{mills} '.format(**rule_paths) + # path of the reference file
'-resource {dbsnp}: '.format(**soft['var_indel_recal']) + # parameter string
'{dbsnp} '.format(**rule_paths) + # path of the reference file
# add some parameters from the configuration that are specific for this step
'--max-gaussians {max-gaussians} -mode {mode} '.format(**soft['var_indel_recal']) +
# transform the list of annotations in a '-an XXXX -an XXXX -an XXXX' list
'-an {} '.format(' -an '.join(soft['var_indel_recal']['an_list'])) +
# dump all standard and error output to the log file
'&> {var_indel_recal_out}'.format(**paths)
# replace all instances of '{md5}' with '{wildcards.md5}'
).replace('{md5}', '{wildcards.md5}')
# apply indel variant recalibration step
rule apply_indel_recal:
input:
var_indel_recal = paths['var_indel_recal'],
var_indel_recal_tranches = paths['var_indel_recal_tranches'],
hc_geno_vcf = paths['hc_geno_vcf']
output:
recal_indel_vcf = paths['recal_indel_vcf']
log:
recal_indel_vcf_out = paths['recal_indel_vcf_out']
shell:
# wrap the whole command in one string
(
# append the Singularity container call (if needed)
SINGULARITY_CMD +
# create the command for GATK with the amount of memory specified in the config
'gatk --java-options "-Xmx{memory}G '.format(**soft['gatk']) +
# specify the GATK JAR path, the GATK command and the reference genome
'ApplyVQSR -R {ref_genome} '.format(**rule_paths) +
# specify the inputs
'-V {hc_geno_vcf} --tranches-file {var_indel_recal_tranches} '.format(**rule_paths) +
'--recal-file {var_indel_recal} '.format(**rule_paths) +
# specify the output
'-O {recal_indel_vcf} '.format(**rule_paths) +
# add some parameters from the configuration that are specific for this step
'--ts_filter_level {ts_filter_level} -mode {mode} '.format(**soft['apply_indel_recal']) +
# dump all standard and error output to the log file
'&> {recal_indel_vcf_out}'.format(**paths)
# replace all instances of '{md5}' with '{wildcards.md5}'
).replace('{md5}', '{wildcards.md5}')
c4science · Help