Exercise 1A 5 – BIOI2 – Integrative BIOInformatics platforme

Getting started with Snakemake

About this course | Before the session | About Snakemake | Course material | Exercises

Exercise 1A - create your first snakefile

objective > setup > o1 > o2 > o3 > o4 > o5 > o6 > recap

Objective 5

Create a new snakefile named ex1_o5.smk in which we add an extra input to the fastqc rule.

Where to start?

Tired of explicitly writing all input and output file names?

=> Use Snakemake’s expand() function to manage all your input RNA-seq files at once.

So, how should we go about it?

Create a list: create a Python list at the beginning of the snakefile containing all the base names of your input files (don’t include the .fastq.gz suffix)
NB: in Python, a list of strings can be defined like this:
list_name = ["item1","item2",...,"itemN"]
Integrate expand(): replace the list of file names in the input and output directives of your rules by the expand() function.

For example, the input directive of your fastqc rule could look like this:

Explanation: This function expects a string as input defining the file paths, in which you replace the variable parts by “placeholders” called wildcards. In this case, there is only one, that we named sample and that should be written between braces: {sample}.

As further inputs to this function, we also have to specify by what elements we would like to replace each wildcard with. In this case, we have to give the function our list of base names (sample=list_name)

more about expand()

Solution

Your code for ex1_o5.smk should look like this (we added the SRR3099587_chr18 sample to the previous 2):

SAMPLES=["SRR3099585_chr18","SRR3099586_chr18","SRR3099587_chr18"]

rule all:
  input:
    expand("FastQC/{sample}_fastqc.html", sample=SAMPLES),    
    expand("FastQC/{sample}_fastqc.zip", sample=SAMPLES), 
    "multiqc_report.html",
    "multiqc_data",

rule fastqc:
  input:
    expand("Data/{sample}.fastq.gz", sample=SAMPLES)
  output:
    expand("FastQC/{sample}_fastqc.zip", sample=SAMPLES),
    expand("FastQC/{sample}_fastqc.html", sample=SAMPLES)
  shell: "fastqc --outdir FastQC {input}"

rule multiqc:
  input:
    expand("FastQC/{sample}_fastqc.zip", sample = SAMPLES)
  output:
    "multiqc_report.html",
    directory("multiqc_data")
  shell: "multiqc {input}"

Test the script

Next, let’s run your pipeline again:

You should see something similar to the following output on your screen:

Observe the output

Snakemake detects that the output files of SRR3099587_chr18 are missing (“reason: Missing output files: FastQC/SRR3099587_chr18_fastqc.html, FastQC/SRR3099587_chr18_fastqc.zip“) and also that “Set of input files has changed since last execution” so it re-ran the fastqc rule.

Have a look at your output folder, you should now have 6 files in there:

Have you noticed that, up until now, Snakemake re-executed FastQC on all inputs every time the fastqc rule was run, indifferent of the fact that the output files were already generated or not?

No? Have a closer look at the output log that you just got… In fact, if you look closely, fastqc is always run only once but on a variable number of inputs (depending on the input of our rule):

Snakemake re-runs fastqc on all inputs because we gave as input to the fastqc rule a list of files rather than just a single file. However, this is sub-optimal, especially on a computer cluster on which you have access to a large amount of available resources to run your jobs on. We’ll see how to change this in the next objective.