Skip to content

snakemake.md

Latest commit

 

History

History
263 lines (201 loc) · 9.17 KB

snakemake.md

File metadata and controls

263 lines (201 loc) · 9.17 KB
Raw

Introduction

Snakemake is a workflow management system (WMS)

It's based on a 'pull' kind of system. You ask for output files and snakemake will use your Snakefile to figure out how to generate the requested output files (pull). With a bash pipeline you provide the input files and bash will use the bash script to generate the output files (push)

Snakefile tips

expand()

expand('{sample}.{direction}.fastqc.html', sample = ['sample1','sample2'], direction = ['fw', 'rv'])

# This will generate
sample1.fw.fastqc.html
sample1.rv.fastqc.html
sample2.fw.fastqc.html
sample2.rv.fastqc.html

# Useful to compactly state your desired output files in the all rule

input functions

# input function
def get_fastq(wildcards):
    if wildcards.sample == 'sample1':
        return 'sample1/sample1.fw.fastq.gz'
    elif wildcards.sample == 'sample2':
        return 'sample2/sample2.fw.fastq.gz'

# rule using input function
## needs to be stated below the input function in the snakefile
rule fastqc:
    input: get_fastq

# Useful if your rule requires a single input file, but the input file
# is different depending on which upstream rule is used


# input function
def get_pilon_input(wildcards):
    input_dict = {
        'assembly' : '{wildcards.assembly}.fasta'.format(wildcards = wildcards),
        'bam'      : '{wildcards.assembly}.bam'.format(wildcards = wildcards)
    }
    return input_dict

# rule using input function
## needs to be stated below the input function in the snakefile
rule pilon:
    input: unpack(get_pilon_input)
    shell: 'pilon --genome {input.assembly} --bam {input.bam}'

# Useful if your rules requires multipe input files, but that set of input files
# is different depending on which upstream rule is used

python code in your rules

# replace function
DB='panther17.hmm.h3f'
rule hmmscan:
    input: DB
    params: DB.replace('.h3f','')
    shell: 'hmmscan {params.db}'

# lambda expression
rule repeatmasker:
    threads: 8
    params:
        jobs=lambda wildcards, threads : threads // 4
    shell: 'RepeatModeler -p {params.job}'

# Useful if there is a maximum amount of threads you can use per 'job'
# // means floor division. E.g., 9 // 4 = 2

TAB management

Snakemake seems to be fairly sensitive to inconsistent TAB indentations For example if some are indented with literal TAB characters (^I), and others with 4 spaces, things may fail

In Vim, use :retab to make it all consistent

local rules

localrules: all, rule_a, rule_b

# These are rules which will not be submitted to compute nodes on a computer cluster.
# In other words, these are rules executed by the head-node / login-node.

temp()

output: temp('outfile')

# Mark output files with temp() to delete them once all rules
# that use those output files are completed

using {} in shell

shell: 'awk {{ print $0 }}'

shell:
    '''
     function body {{
         IFS= read -r header
         printf '%s\n' "$header"
         "$@"
     }}
    '''

# { and } have special meaning in the shell directive,
# so if they are used in a command you need to add another { and } to escape the special meaning?

Snakemake command line

# Run a snakefile that is not called `Snakefile`
snakemake -S some_snakefile.smk

# Generate a jobgraph and rulegraph
snakemake --dag       | dot -Tpng > jobgraph.png
snakemake --rulegraph | dot -Tnpg > rulegraph.png

# Check `Snakefile` for syntax errors, pipeline connectiveness and which rules and jobs will be executed
snakemake --dry-run / -n
## NOTE: when doing a dry-run the yellow text will report using only 1 thread for a job,
## even if you requested multiple threads in the relevant rule. Not to worry,
## snakemake will execute a rule with the requested number of threads in a normal run

# State reason why each rule was activated
snakemake --reason / -r

# State exact shell commands used for each rule or job
snakemake --printshellcmds / -p

# Set the total number of cores that can be used
snakemake --cores 40

# Use conda environments defined in the rules
## Can be pre-existing conda environment, or an environment YAML file
snakemake --use-conda

# Installs all necessary conda packages, but does NOT run the pipeline
snakemake --conda-create-envs-only

# Activate a conda environment created by snakemake
## Here <env> is the MD5 hash of the environment definition file content ??
source activate .snakemake/conda/<env>

# Define conditions for rerunning a rule. 
## Default is mtime,params,input,software_env,code, which is very aggressive
snakemake --rerun-triggers mtime
snakemake --rerun-triggers mtime code

# Mark a file as up to date, even if its mtime is older than upstream files
snakemake --touch file

# Unlock a working directory
## Snakemake workdirs can get locked if a pipeline crashed in an abnormal way
## For example Ctrl-C, power outage, HPC job that ran out of time
snakemake --unlock

Executing snakemake on a HPC cluster

# Snakemake is able to delegate snakejobs as separate jobs to the queue of a computer cluster
snakemake --cluster 'qsub -V'
snakemake --cluster 'sbatch'
# If you are on a Sun Grid Engine cluster (qsub), make sure you are in
# the environment that has all required tools loaded before executing the snakefile
# `-V` ensures that your environment variables are passed on to the submitted job

# Submit each snakejob with the number of threads specified in the rule
snakemake --cluster 'qsub -V -pe threaded {threads}'

# Limit the number of jobs that can be submitted to the HPC queue at the same time
snakemake --cluster 'qsub -V' --jobs 6

# Specify the location of the STDOUT and STDERR files of each submitted snakejob
snakemake --cluster 'qsub -V -cwd -o logs/{rule}.{jobid}.o -e logs/{rule}.{jobid}.e'
## Here {rule} and {jobid} are special wildcards holding the rule name and jobid of the snakejob, respectively.
## If not specified, the STDOUT and STDERR files are stored in your $HOME.
## -cwd ensures that the relative paths of -o and -e are relative to your working dir instead relative to your $HOME

Experienced issues with snakemake on a HPC cluster

If you have rules that point to specific conda environments, e.g. conda: 'funannotate',
then make sure that there are no conflicting versions of softwares between 
the environment where you are submitting the snakemake pipeline from (e.g. 'proj-ergo', or 'snakemake')
and the conda environment. 

I had one issue where funannotate was calling the wrong version of diamond,
even though the 'funannotate' environment had the correct version installed,
because the environment that I was calling snakemake from, had an earlier version
of diamond, and in the jobs environment that version's path preceeded the correction version's path
in $PATH

Using conda environments on a HPC cluster

If you have a rule with `conda: env.yaml`, snakemake will download 
and install the requested packages under `.snakemake/conda/<env>` 
where `<env>` is an MD5 hash. The rule will be run with the software installed there.

If you put an existing conda environment in the conda: directive, and you get an 
EnvironmentNameNotFound error, try putting the absolute path of the desired
conda environment in the directive

conda: '/scratch2/software/anaconda/envs/signalp6-fast'

UPDATE: the above strategy seems to work in the sense that it now executes the code in the shell: directive.
However, upon ending the execution of the rule, a new error occurrs
subprocess.CalledProcessError: Command 'conda env export --name /scratch2/software/anaconda/envs/signalp6-fast returned non-zero exit status
Because the conda env export --name does not accept, among others '/' characters...

I found that the reason the workflow does not recognize the existing conda environment has to do with the fact that,
while in the snakemake environment on Perun, the conda executable that is used is
/scratch2/software/anaconda/envs/snakemake/bin/conda

while the environment that I wish to activate is managed by a different conda installation:
/scratch2/software/anaconda/condabin/conda

To ensure that the workflow finds the existing conda environment, add this on top of the snakemake file
shell("conda config --add envs_dirs /scratch2/software/anaconda/envs")
This created (it didn't exist before) a .condarc file on my home, containing

envs_dirs:
  - /scratch2/software/anaconda/envs

So I guess just making that file on your own would have been sufficient,
the shell() line in the snakefile is perhaps not necessary


Upgrading a package in a snakemake conda environment:
1. First activate the conda environment
source activate .snakemake/conda/880c7b2708595faf94f84838fc3802df

2. Then run
/scratch2/software/anaconda/envs/snakemake/bin/mamba install -c bioconda diamond=2.1.6

NOTE: I have specifically called the mamba install from the snakemake environment!

Other notes

# Output directories specified in your rules are automatically made my Snakemake. 
# So no need for mkdir in the shell: directive

# After a rule is executed, snakemake will check if all your desired output files 
# stated in output: are actually present in your filesystem
# If at least one of them is not present, snakemake will assume the rule executation failed 
# and remove all files generated by that rule / job