snakemake/index.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Intro to workflows for data analysis, using snakemake\n",
    "### Authors: Meghan Correa and Ming Tang\n",
    "### Bioinformatics Coffee Hour - May 26, 2020\n",
    "\n",
    "Modified from https://github.com/ctb/2019-snakemake-ucdavis\n",
    "\n",
    "No license; the below content is under CC0. (Do with it what you will, and I hope it's useful!)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction and thoughts\n",
    "\n",
    "Hi! I'm Meghan Correa, I a member of the Software Operations team at FAS informatics. I support data analysis for the Bauer core among other things.\n",
    "\n",
    "This 25min teaser meant to show you how to get started using snakemake.\n",
    "\n",
    "Note, the only way you'll really learn to do all of this is by applying it to your own research and spending time on it :)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Why use a workflow management tool?\n",
    "\n",
    "* Dependency management\n",
    "* Reentrancy - start back up where you left off\n",
    "* Reusable\n",
    "* Documented\n",
    "* Portable\n",
    "\n",
    "There are many to choose from:\n",
    "https://github.com/pditommaso/awesome-pipeline\n",
    "\n",
    "FAS Informatics choose to use snakemake for the post sequencing pipeline and other workflows mainly because it is written in python and is popular in the bioinformatics community."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tasks we are going to do\n",
    "\n",
    "Inside the `data` folder, there are 4 fastq.gz files. We will use `Snakemake` do `fastqc` on each fastq files and compile multiple html files into a single one using `multiQC`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Software we're going to use\n",
    "\n",
    "We're going to be using [conda](https://conda.io/en/latest/) and [snakemake](https://snakemake.readthedocs.io/en/stable/), a well as packages from [bioconda](https://bioconda.github.io). If you wanted to run all of this on your own computer, you'll need to follow the bioconda install instructions.\n",
    "\n",
    "We'll be implementing a short read quality check pipeline, using [fastqc](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and [multiqc](https://multiqc.info/). No worries if you don't know what any of this means, it's not super critical to the snakemake side of things :)\n",
    "\n",
    "You can see the full set of installed software requirements for python in `environment.yml` file and installs from apt repository in `apt.txt` in the binder folder of the repository.\n",
    "\n",
    "**Note we have installed all the required tools in the binder container**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Download some data\n",
    "\n",
    "Execute the below cell to fetch some data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "curl -L https://github.com/ctb/2019-snakemake-ucdavis/raw/9db09bc0b6a3469f8a0d4996d4b2995bf36e5d27/data/0Hour_001_1.fq.gz > data/0Hour_001_1.fq.gz\n",
    "curl -L https://github.com/ctb/2019-snakemake-ucdavis/raw/9db09bc0b6a3469f8a0d4996d4b2995bf36e5d27/data/6Hour_001_1.fq.gz > data/6Hour_001_1.fq.gz\n",
    "curl -L https://github.com/ctb/2019-snakemake-ucdavis/raw/9db09bc0b6a3469f8a0d4996d4b2995bf36e5d27/data/0Hour_001_2.fq.gz > data/0Hour_001_2.fq.gz\n",
    "curl -L https://github.com/ctb/2019-snakemake-ucdavis/raw/9db09bc0b6a3469f8a0d4996d4b2995bf36e5d27/data/6Hour_001_2.fq.gz > data/6Hour_001_2.fq.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running snakemake!\n",
    "\n",
    "### Getting started - your first Snakefile\n",
    "\n",
    "Open the `Snakemake` file with a text editor and paste the below into it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rule fastqc_a_file:\n",
    "  shell:\n",
    "    \"fastqc data/0Hour_001_1.fq.gz\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(I suggest copy/pasting this into the Snakefile.)\n",
    "\n",
    "remember to click `File` --> `Save File`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "snakemake"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and you should see:\n",
    "```\n",
    "Building DAG of jobs...\n",
    "Using shell: /bin/bash\n",
    "Provided cores: 1\n",
    "Rules claiming more threads will be scaled down.\n",
    "...\n",
    "Approx 95% complete for 0Hour_001_1.fq.gz\n",
    "Analysis complete for 0Hour_001_1.fq.gz\n",
    "[Wed Feb 27 13:09:51 2019]\n",
    "Finished job 0.\n",
    "1 of 1 steps (100%) done\n",
    "Complete log: /home/jovyan/.snakemake/log/2019-02-27T130941.260352.snakemake.log\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and there will be two new files,\n",
    "\n",
    "Points to make:\n",
    "* the snakemake configuration file is by default called `Snakefile`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Updating the Snakefile to track inputs and outputs\n",
    "\n",
    "At the moment this is basically just a shell script with extra syntax... what's the point?\n",
    "\n",
    "Well, shell scripts will rerun the command every time you run the file, even if there's no reason to do so because the file hasn't changed but with snakemake, you can annotate the rule with input and output files!\n",
    "\n",
    "**Digression:** This is particularly important for large or long workflows, where you're dealing with 10s to 100s of files that may take hours to days to process! It can be hard to figure out which files to rerun, but (spoiler alert) snakemake can really help you do this!\n",
    "\n",
    "Change your snakefile to look like this:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rule fastqc_a_file:\n",
    "  input:\n",
    "    \"data/0Hour_001_1.fq.gz\"\n",
    "  output:\n",
    "    \"data/0Hour_001_1_fastqc.html\",\n",
    "    \"data/0Hour_001_1_fastqc.zip\"\n",
    "  shell:\n",
    "    \"fastqc data/0Hour_001_1.fq.gz\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "here, we've annotated the rule with the required\n",
    "**input** file, as well as the expected **output** files.\n",
    "\n",
    "Question: how do we know what the output files are?\n",
    "\n",
    "Now run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "snakemake"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and you should see:\n",
    "```\n",
    "Building DAG of jobs...\n",
    "Nothing to be done.\n",
    "Complete log: /home/jovyan/.snakemake/log/2019-02-27T132031.813143.snakemake.log\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What happened??\n",
    "\n",
    "snakemake looked at the file, saw that the output files existed, and figured out that it didn't need to do anything!\n",
    "\n",
    "### Forcibly re-running things\n",
    "\n",
    "You can tell snakemake to run all the rules no matter what with `-F`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "snakemake -F"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also remove an output file and it will automatically re-run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rm data/*.html\n",
    "snakemake"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "note that you don't need to remove *all* the output files to rerun a command - just remove *one* of them.\n",
    "\n",
    "You can *also* update the timestamp on an *input* file, and snakemake will figure out that the output file is older than the input file, and rerun things."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "touch data/*.fq.gz\n",
    "snakemake"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another usefule snakemake option is the --dryrun option, this will show you what snakemake would do but won't actually run any of the rules.  I highly recommend running this to check that was will be run is what you expect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "snakemake --dryrun"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multiple rules\n",
    "\n",
    "Let's add a rule to run fastqc on a second file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rule fastqc_a_file:\n",
    "  input:\n",
    "    \"data/0Hour_001_1.fq.gz\"\n",
    "  output:\n",
    "    \"data/0Hour_001_1_fastqc.html\",\n",
    "    \"data/0Hour_001_1_fastqc.zip\"\n",
    "  shell:\n",
    "    \"fastqc data/0Hour_001_1.fq.gz\"\n",
    "\n",
    "rule fastqc_a_file2:\n",
    "  input:\n",
    "    \"data/6Hour_001_1.fq.gz\"\n",
    "  output:\n",
    "    \"data/6Hour_001_1_fastqc.html\",\n",
    "    \"data/6Hour_001_1_fastqc.zip\"\n",
    "  shell:\n",
    "    \"fastqc data/6Hour_001_1.fq.gz\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "snakemake"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, if you run this, the Right Thing won't happen: snakemake will do nothing.  Why?\n",
    "\n",
    "Well, snakemake only runs the *first* rule in a Snakefile, by default. You can give a rule name on the command line, if you like, **or** you can tell snakemake what output file(s) you want. "
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "snakemake data/0Hour_001_1_fastqc.html data/6Hour_001_1_fastqc.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But the conventional way to do this in snakemake is to add an \"all\" rule as the first rule, let's try it. \n",
    "\n",
    "\n",
    "### A first refactoring: adding a better default rule\n",
    "\n",
    "Let's start refactoring (cleaning up) this Snakefile.\n",
    "\n",
    "First, let's add a rule at the top:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rule all:\n",
    "  input:\n",
    "    \"data/0Hour_001_1_fastqc.html\",\n",
    "    \"data/6Hour_001_1_fastqc.html\"\n",
    "\n",
    "rule fastqc_a_file:\n",
    "  input:\n",
    "    \"data/0Hour_001_1.fq.gz\"\n",
    "  output:\n",
    "    \"data/0Hour_001_1_fastqc.html\",\n",
    "    \"data/0Hour_001_1_fastqc.zip\"\n",
    "  shell:\n",
    "    \"fastqc data/0Hour_001_1.fq.gz\"\n",
    "\n",
    "rule fastqc_a_file2:\n",
    "  input:\n",
    "    \"data/6Hour_001_1.fq.gz\"\n",
    "  output:\n",
    "    \"data/6Hour_001_1_fastqc.html\",\n",
    "    \"data/6Hour_001_1_fastqc.zip\"\n",
    "  shell:\n",
    "    \"fastqc data/6Hour_001_1.fq.gz\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "snakemake"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a blank rule that gathers together all of the various files you want produced, and says \"hey, snakemake, I depend on all of these files for my input - make them for me!\"\n",
    "\n",
    "Let's run snakemake again and now you should see the second fastqc command run, with the appropriate output files!\n",
    "\n",
    "Note that snakemake only runs the second rule, because it looks at the output files and sees that the first file you wanted, `0Hour_001_1_fastqc.html` already exists!\n",
    "\n",
    "This is a blank rule that gathers together all of the various files you want produced, and says \"hey, snakemake, I depend on all of these files for my input - make them for me!\"\n",
    "\n",
    "### A second refactoring: doing a bit of templating\n",
    "\n",
    "There's a lot of repetition in each of these rules. Let's collapse it down a little bit by replacing the filename in the fastqc command with a magic variable, `{input}`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rule all:\n",
    "  input:\n",
    "    \"data/0Hour_001_1_fastqc.html\",\n",
    "    \"data/6Hour_001_1_fastqc.html\"\n",
    "\n",
    "rule fastqc_a_file:\n",
    "  input:\n",
    "    \"data/0Hour_001_1.fq.gz\"\n",
    "  output:\n",
    "    \"data/0Hour_001_1_fastqc.html\",\n",
    "    \"data/0Hour_001_1_fastqc.zip\"\n",
    "  shell:\n",
    "    \"fastqc {input}\"\n",
    "\n",
    "rule fastqc_a_file2:\n",
    "  input:\n",
    "    \"data/6Hour_001_1.fq.gz\"\n",
    "  output:\n",
    "    \"data/6Hour_001_1_fastqc.html\",\n",
    "    \"data/6Hour_001_1_fastqc.zip\"\n",
    "  shell:\n",
    "    \"fastqc {input}\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This all works as before, but now the rule is a bit more generic and will work with any input file. Sort of.\n",
    "\n",
    "### Refactoring 3: templating output files, too\n",
    "\n",
    "The output filenames ALSO depend on the input file names in some way - specifically, fastqc replace part of the filename with `_fastqc.html` and `_fastqc.zip` to make its two output files.\n",
    "\n",
    "Let's rewrite the rule using some snakemake pattern matching:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rule all:\n",
    "  input:\n",
    "    \"data/0Hour_001_1_fastqc.html\",\n",
    "    \"data/6Hour_001_1_fastqc.html\"\n",
    "\n",
    "rule fastqc_a_file:\n",
    "  input:\n",
    "    \"{filename}.fq.gz\"\n",
    "  output:\n",
    "    \"{filename}_fastqc.html\",\n",
    "    \"{filename}_fastqc.zip\"\n",
    "  shell:\n",
    "    \"fastqc {input}\"\n",
    "\n",
    "rule fastqc_a_file2:\n",
    "  input:\n",
    "    \"{filename}.fq.gz\"\n",
    "  output:\n",
    "    \"{filename}_fastqc.html\",\n",
    "    \"{filename}_fastqc.zip\"\n",
    "  shell:\n",
    "    \"fastqc {input}\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we've done here is tell snakemake that anytime we say we *want* a file that ends with `_fastqc.html`, it should look for a file that ends in `.fq.gz` and then run `fastqc` on it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you look at the rule above, they are now identical rules! \n",
    "\n",
    "Let's remove one, to get a trimmer, leaner, and above all *functional* snakefile:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rule all:\n",
    "  input:\n",
    "    \"data/0Hour_001_1_fastqc.html\",\n",
    "    \"data/6Hour_001_1_fastqc.html\"\n",
    "\n",
    "rule fastqc_a_file:\n",
    "  input:\n",
    "    \"{filename}.fq.gz\"\n",
    "  output:\n",
    "    \"{filename}_fastqc.html\",\n",
    "    \"{filename}_fastqc.zip\"\n",
    "  shell:\n",
    "    \"fastqc {input}\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the variable name in input and output does not have to be \"filename\", it can be anything as long as you're consistant.  So let's run snakemake again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rm data/*.html\n",
    "snakemake"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Building out the workflow\n",
    "\n",
    "So, we've gotten fastqc sorted out. What's next?\n",
    "\n",
    "Let's add in a new rule - multiqc, to summarize our fastqc results.\n",
    "\n",
    "multiqc takes a directory name under which there are one or more fastqc reports, and then summarizes them.\n",
    "\n",
    "Running it on the command line,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "multiqc data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "you can see that it creates two outputs, `multiqc_report.html` and the directory `multiqc_data/` which contains a bunch of files. Let's create a snakemake rule for this; add:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rule all:\n",
    "  input:\n",
    "    \"data/0Hour_001_1_fastqc.html\",\n",
    "    \"data/6Hour_001_1_fastqc.html\"\n",
    "\n",
    "rule fastqc_a_file:\n",
    "  input:\n",
    "    \"{filename}.fq.gz\"\n",
    "  output:\n",
    "    \"{filename}_fastqc.html\",\n",
    "    \"{filename}_fastqc.zip\"\n",
    "  shell:\n",
    "    \"fastqc {input}\"\n",
    "\n",
    "rule run_multiqc:\n",
    "  input:\n",
    "    \"data/0Hour_001_1_fastqc.html\",\n",
    "    \"data/6Hour_001_1_fastqc.html\",\n",
    "  output:\n",
    "    \"multiqc_report.html\",\n",
    "    directory(\"multiqc_data\")\n",
    "  shell:\n",
    "    \"multiqc data/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "snakemake"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This ...doesn't really do what we want, for a few reasons.\n",
    "\n",
    "First of all, the output of run_multiqc is not specified in the all rule so snakemake doesn't look for a rule to create this. \n",
    "\n",
    "Second of all, `multiqc_report.html` already exists, so snakemake won't run the rule. \n",
    "\n",
    "Let's fix the first two things first:\n",
    "\n",
    "* add `multiqc_report.html` to the inputs for the first all.\n",
    "* then remove `multiqc_report.html` and re-run snakemake.\n",
    "\n",
    "Your snakefile should look like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rule all:\n",
    "  input:\n",
    "    \"data/0Hour_001_1_fastqc.html\",\n",
    "    \"data/6Hour_001_1_fastqc.html\",\n",
    "    \"multiqc_report.html\"\n",
    "\n",
    "rule fastqc_a_file:\n",
    "  input:\n",
    "    \"{filename}.fq.gz\"\n",
    "  output:\n",
    "    \"{filename}_fastqc.html\",\n",
    "    \"{filename}_fastqc.zip\"\n",
    "  shell:\n",
    "    \"fastqc {input}\"\n",
    "\n",
    "rule run_multiqc:\n",
    "  input:\n",
    "    \"data/0Hour_001_1_fastqc.html\",\n",
    "    \"data/6Hour_001_1_fastqc.html\",\n",
    "  output:\n",
    "    \"multiqc_report.html\",\n",
    "    directory(\"multiqc_data\")\n",
    "  shell:\n",
    "    \"multiqc data/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rm multiqc_report.html\n",
    "snakemake"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yay, that seems to work!\n",
    "\n",
    "Points to make:\n",
    "\n",
    "* other than the first rule, rules can be in any order\n",
    "* the rule name doesn't really matter, it's mostly for debugging. It just needs to be \"boring\" (text, underscores, etc. only)\n",
    "\n",
    "\n",
    "But providing input files explicitly to the multiqc rule is not great because those are the same files we have in the all rule and to add new files we'd now have to add them to two places. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can fix the first issue by using **variables**. \n",
    "\n",
    "To use variables, let's make a Python list at the very top, containing all of our expected output files from fastqc:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fastqc_output = [\"data/0Hour_001_1_fastqc.html\", \"data/6Hour_001_1_fastqc.html\",\n",
    "  \"data/0Hour_001_2_fastqc.html\", \"data/6Hour_001_2_fastqc.html\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and modify the `all` and `multiqc` rules to contain this list. The final snakefile looks like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fastqc_output = [\"data/0Hour_001_1_fastqc.html\", \"data/6Hour_001_1_fastqc.html\",\n",
    "  \"data/0Hour_001_2_fastqc.html\", \"data/6Hour_001_2_fastqc.html\"]\n",
    "\n",
    "rule all:\n",
    "  input:\n",
    "    fastqc_output,\n",
    "    \"multiqc_report.html\"\n",
    "\n",
    "rule fastqc_a_file:\n",
    "  input:\n",
    "    \"{filename}.fq.gz\"\n",
    "  output:\n",
    "    \"{filename}_fastqc.html\",\n",
    "    \"{filename}_fastqc.zip\"\n",
    "  shell:\n",
    "    \"fastqc {input}\"\n",
    "\n",
    "rule run_multiqc:\n",
    "  input:\n",
    "    fastqc_output\n",
    "  output:\n",
    "    \"multiqc_report.html\",\n",
    "    directory(\"multiqc_data\")\n",
    "  shell:\n",
    "    \"multiqc data/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "snakemake"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Refactoring this to make it slightly more concise --\n",
    "\n",
    "We've got one more redundancy in this file - the `fastqc_output` is listed in the `all` rule, but you don't need it there! Why?\n",
    "\n",
    "Well, `multiqc_report.html` is already in the all rule, and the multiqc rule depends on `fastqc_output`, so `fastqc_output` already needs to be created to satisfy the all rule, so... specifying it in the all rule is redundant! And you can remove it!\n",
    "\n",
    "(It's not a big deal and I usually leave it in. But I wanted to talk about dependencies!)\n",
    "\n",
    "The Snakefile now looks like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fastqc_output = [\"data/0Hour_001_1_fastqc.html\", \"data/6Hour_001_1_fastqc.html\",\n",
    "  \"data/0Hour_001_2_fastqc.html\", \"data/6Hour_001_2_fastqc.html\"]\n",
    "\n",
    "rule all:\n",
    "  input:\n",
    "    \"multiqc_report.html\"\n",
    "\n",
    "rule fastqc_a_file:\n",
    "  input:\n",
    "    \"{filename}.fq.gz\"\n",
    "  output:\n",
    "    \"{filename}_fastqc.html\",\n",
    "    \"{filename}_fastqc.zip\"\n",
    "  shell:\n",
    "    \"fastqc {input}\"\n",
    "\n",
    "rule run_multiqc:\n",
    "  input:\n",
    "    fastqc_output\n",
    "  output:\n",
    "    \"multiqc_report.html\",\n",
    "    directory(\"multiqc_data\")\n",
    "  shell:\n",
    "    \"multiqc data/\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and we can rerun it from scratch by doing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rm data/*.html multiqc_report.html\n",
    "snakemake"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Recap\n",
    "\n",
    "So, we've put all this work into making this snakefile with its input rules and its output rules... and there are a lot of advantages to our current approach already! Let's list a few of them --\n",
    "\n",
    "* we've completely automated our analysis!\n",
    "* we can easily add new data files into fastqc and multiqc!\n",
    "* we can rerun things easily, and (even better) by default only things that *need* to be run will be run.\n",
    "* the snakefile is actually pretty reusable - we could drop this into a new project, and, with little effort, run all of these things on new data!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## More advanced snakemake"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running things in parallel\n",
    "\n",
    "You can use `snakemake --cores 4` to run four jobs in parallel. If you are running on a cluster then --cores becomes the number of jobs that can be submitted at a time.  \n",
    "\n",
    "### Specifying software required for a rule\n",
    "\n",
    "You can specify software on a per-rule basis! This is really helpful when you have incompatible software requirements for different rules, or want to run on a cluster, or just want to pin your snakemake workflow to a specific version.\n",
    "\n",
    "For example, if you create a file `env_fastqc.yml` with the following content,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "channels:\n",
    "  - bioconda\n",
    "  - defaults\n",
    "  - conda-forge\n",
    "dependencies:\n",
    "  - fastqc==0.11.8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and then change the fastqc rule to look like this:"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "rule fastqc_a_file:\n",
    "  input:\n",
    "    \"{filename}.fq.gz\"\n",
    "  output:\n",
    "    \"{filename}_fastqc.html\",\n",
    "    \"{filename}_fastqc.zip\"\n",
    "  conda:\n",
    "    \"env_fastqc.yml\"\n",
    "  shell:\n",
    "    \"fastqc {input}\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "you can now run snakemake like so,"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "snakemake --use-conda"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and for that rule, snakemake will install just that software, with the specified version.\n",
    "\n",
    "This aids in reproducibility, in addition to the practical advantages of isolating software installs from each other.\n",
    "\n",
    "(You can also do this with docker and singularity containers, too!)\n",
    "\n",
    "### Running on a cluster\n",
    "\n",
    "You can specify a cluster submit command:\n",
    "\n",
    "snakemake --cluster sbatch -j 32\n",
    "\n",
    "You can also specify per rule cluster config in a file (example below):\n",
    "\n",
    "snakemake --cluster sbatch --cluster-config cluster-config.json -j 32 \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "{\n",
    "   \"default\":{\n",
    "      \"time\":\"12:00:00\",\n",
    "      \"nodes\":1,\n",
    "      \"cores\":24,\n",
    "      \"mem\":100,\n",
    "      \"partition\":\"shared\",\n",
    "      \"job-name\":\"{rule}\",\n",
    "      \"output\":\"log/{rule}-%j.out\",\n",
    "      \"error\":\"log/{rule}-%j.err\"\n",
    "   },\n",
    "   \"fastqc\":{\n",
    "      \"time\":\"4-12:00:00\",\n",
    "      \"mem\":32000\n",
    "   },\n",
    "   \"multiqc\":{\n",
    "      \"mem\":64000\n",
    "   }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or you can use snakemake profiles and bundle cluster config files in a directory:\n",
    "\n",
    "snakemake --profile slurm -j 32\n",
    "\n",
    "https://snakemake.readthedocs.io/en/v5.18.1/executing/cli.html#profiles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding in some Python...\n",
    "\n",
    "You can add in some Python to load in the input files, like so:"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "import glob, sys\n",
    "fastqc_input = glob.glob('data/?Hour_00?_?.fq.gz')\n",
    "\n",
    "fastqc_output = []\n",
    "for filename in fastqc_input:\n",
    "  new_filename = filename.split('.')[0] + '_fastqc.html'\n",
    "  fastqc_output.append(new_filename)\n",
    "\n",
    "rule all:\n",
    "  input:\n",
    "    \"multiqc_report.html\"\n",
    "\n",
    "rule clean:\n",
    "  shell:\n",
    "    \"rm -f {fastqc_output} multiqc_report.html\"\n",
    "\n",
    "rule fastqc_a_file:\n",
    "  input:\n",
    "    \"{arglebarf}.fq.gz\"\n",
    "  output:\n",
    "    \"{arglebarf}_fastqc.html\",\n",
    "    \"{arglebarf}_fastqc.zip\"\n",
    "  conda:\n",
    "    \"env_fastqc.yml\"\n",
    "  shell:\n",
    "    \"fastqc {input}\"\n",
    "\n",
    "\n",
    "rule run_multiqc:\n",
    "  input:\n",
    "    fastqc_output\n",
    "  output:\n",
    "    \"multiqc_report.html\",\n",
    "    directory(\"multiqc_data\")\n",
    "  shell:\n",
    "    \"multiqc data/\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Go forth and workflow!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Bash",
   "language": "bash",
   "name": "bash"
  },
  "language_info": {
   "codemirror_mode": "shell",
   "file_extension": ".sh",
   "mimetype": "text/x-sh",
   "name": "bash"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}