Snakemake—simple data processing for researchers

How to learn Snakemake from scratch using official materials?

1. Go through the official tutorial and best practices
2. You are probably ready to implement your own pipeline
3. If you need to run Snakemake on HPC or a cloud, go through the executor tutorial and the documentation about the cluster/cloud execution
4. Check out the rest of the documentation once you get a better understanding of Snakemake
5. Keep an eye on Snakemake's changelog for new features

There are quite a few completely opposite properties about these two scenarios, which are listed in the table.

For the reproducible analysis scenario, the computing environment is known and limited only when we try to reproduce it on our own machines. In my experience, having a clear and simple instructions is easier for other researchers to work with the analysis than say a fully automated and dockerized setup.

The task dependency graph has the property as a directed acyclic graph (DAG). For Snakemake, such graph is constructed backwards from the final outputs, also known as the pull-based approach.

Reading the pipeline backwards can be a bit counterintuitive at first. My tip is to always start with the final outputs (usually at the end of the file) and trace the workflow backwards.

I implemented two RNA-seq pipelines in Snakemake that are lossly based on GDC's mRNA analysis pipeline:

• https://github.com/ding-lab/cptac_rna_expression that starts from GDC aligned BAMs
• https://github.com/ding-lab/HTAN_bulkRNA_expression that starts from FASTQs

These pipelines were developed a while ago, when some of my computing environments were quite old and limited. So I took a relative conservative approach by using minimal features and very few dependencies (just the Python standard library). These pipelines can certainly be improved to better conform to the community's best practices. With the ubiquitous conda and Docker support and new Snakemake features, some possible improvements include:

• Replace {WORKFLOW_ROOT} with workflow.source_path()
• Use pandas to process data frames
• Modularization, though my pipelines are very short (less than 300 lines)

If you have wrote a Makefile before, snakemake rules share the same spirit as a Makefile's rule. And they run in a similar way.

Snakemake syntax explained

{sample}

Known as wildcards that matches the output file pattern with the existing files.

input:, output:

File patterns for inputs and outputs. Use key-value pairs key='...' when there are multiple patterns. They are just Python strings and thus can be referred elsewhere by rules.myrule.output.key.

shell:

The command that will run in a bash shell to generate the desired outputs.

Check out the documentation for the rest parts of the rule.

How snakemake creates a task

Due to the way snakemake executes, it's easier to resolve the output first by realizing all the wildcard variables. Then expand the rest of the rule (inputs, parameters, and etc.)

For example, to generate readcount/my_sample_a.tsv using this task, {sample} expands to my_sample_a in the rest parts of the rule. So Snakemake will look for the corresponding inputs and run the following shell commands:

featureCounts -T 16 -a /path/to/genes.gtf \
	-o readcount/my_sample_a.tsv \
	inputs/my_sample_a.bam \
	2> logs/readcount/my_sample_a.log
		

Process multiple samples using snakemake all

To make the pipeline more extensible, use a run-all rule and a config file to process all samples of a batch. As shown here, we obtain the list of samples to process from config file by reading in the list from the file at config['sample_list']. In this case, SAMPLES is a Python list of ["my_sample_a", "my_sample_b", ..., "my_sample_d"]. In the all rule, we use a Snakemake function expand() to generates all the output targets by "filling in" the wildcards of the output pattern using every value of SAMPLES.

Other commonly used features

When the simple string pattern matching is not enough, Snakemake also allows many parts of the rule to be determined by Python functions, including input, params, and resources. These functions are known as input functions in Snakemake's documentation. They make the rule more flexible. As shown here, different samples may have their input files at different locations. The same task will ask for more memory if it has been failed before.

unpack() is a helper function to automatically unpack the return value of the input function. So a function can return multiple key-value pairs.

file_map_df is a pandas DataFrame object loaded from the tabular file specified by the config (config["file_map"]).

I want to mention the dry run mode is very useful too (snakemake -n).

How to set up the conda environment or the Docker image?

I actually had better experience to simply put together all the dependencies in one big fat Docker image or conda environment. This setup can live outside Snakemake, so it doesn't take care of pulling the docker image or creating the conda environment. This comes with some advantages:

• It's easier to debug and learn when one is new to Snakemake
• If one creates a pipeline from scratch, it's likely you will first try out the steps in a conda environment anyway. Then one can re-use the environment
• Some computing environments have very non-standard (broken) ways to use conda and Docker, which breaks Snakemake

However, this approach is against the best practice because it reduce the reusability of the rules and sometimes different rules have incompatible dependencies. One should consider switching to a rule-based environments/images later on.

Advantages of using batch folders, config, file map, and output manifest

Batch folders make the pipeline source code and the outputs live in different locations. Good for version tracking the pipeline and reusing the pipeline. We can also create a template folder for new batches. Config file is a natural solution to specify the batch-specific parameters.

Input file map is useful when the upstream files are tracked externally (say, in a database or some centralized system). It also retains the file metadata (e.g. checksums, UUID, experiment details, and other cross-references). For example, our lab maintains a file catalog to track all CPTAC data on GDC. It also tracks the files we have downloaded to the different storage systems in our lab. Thus, our example pipeline looks up this catalog to find the input BAM files.

Advantages of using checksums everywhere

There are quite a few possible reasons lead to a input change or corruption:

• With increasing data size, it's more likely to get rotten bits from failing storage
• Silent file transfer errors
• Updates from the upstream

In general, it's a good idea to keep checksums all the time. Here is a snippet to generate MD5 checksums of everything under a folder in parallel and verify the content using the generated checksums:

fd --type file --exclude 'README*' --exclude '*CHECKSUMS' \
	| sort \
	| parallel --keep-order -j4 --bar 'md5sum {}' \
	> MD5_CHECKSUMS
md5sum -c MD5_CHECKSUMS