genomics_data_analysis/genome_assembly.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "113df914-b625-46a3-8b82-f522c98f48d8",
   "metadata": {},
   "source": [
    "# Genome assembly and quality control of assembly"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49ff02db-aca0-4676-9da3-330584a10d49",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Introduction \n",
    "\n",
    "After the previous steps of quality control (QC), we have reads still in `.fastq` format, but now we have a summary of their sequencing quality. Furthermore, we have removed regions with poor confidence in base calling (where we cannot be sure if the assigned nucleotides are right) and we removed the adaptor sequences used by the sequencing machine to amplify the DNA. \n",
    "\n",
    "In this series of steps, we will assemble the reads using a tool called `shovill`. Once again, we will mimic how to run the commands in the **Compute Canada (CC)** cluster of analysis. \n",
    "\n",
    "### Required tools and environment set up\n",
    "\n",
    "For these tutorials, tools will be made available using singularity containers, which can be run using the command `singularity run tool_image`. These tools have been made available in the environment already, so there is no need to download them.\n",
    "\n",
    "Tools used in this tutorial:\n",
    "- `Shovill         v1.1.0 (with SPAdes v3.14.1)`\n",
    "- `Quast           v5.0.2`\n",
    "- `Singularity     v3.10.2`\n",
    "- `CheckM          v1.2.2`\n",
    "\n",
    "Tools downloaded for the tutorial are in the a shared folder and in the `tutorials` directory are the results of our commands. \n",
    "\n",
    "**Notes:** \n",
    "- After opening this notebook, make sure that you have the **bash** kernel selected. \n",
    "- Try to run all the commands from your home or base directory, to get there just use `cd`\n",
    "\n",
    "Also, run the source command to make all pre-packed software available in your environment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "496d090f-da9d-4813-a2d7-e0e8d163d1c7",
   "metadata": {
    "tags": [],
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [],
   "source": [
    "# source PATH to use module function\n",
    "source /cvmfs/soft.computecanada.ca/config/profile/bash.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0be6bf2e-4176-49fe-b33c-9195a7a7bda6",
   "metadata": {},
   "source": [
    "Explore the results produced so far using the `tree` command, and get a sense of the structure for the software folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c5d1b2c",
   "metadata": {
    "tags": [],
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [],
   "source": [
    "# get to your HOME directory\n",
    "cd\n",
    "\n",
    "echo -e \"Structure of tutorials folder:\"\n",
    "tree -dL 2 $HOME/tutorials\n",
    "\n",
    "echo -e \"Structure of software directory:\"\n",
    "tree -dL 2 /mnt/cidgoh-object-storage/seagull/jupyter-mdprieto"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "773c3c85",
   "metadata": {},
   "source": [
    "For downstream assembly, we will use the curated reads contained in the `trimmed_reads` subdirectory (20 files containing paired reads for 10 isolates of _P. aeruginosa_). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7787d504-aa7f-47d3-b3c7-00308ecf43c8",
   "metadata": {
    "tags": [],
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [],
   "source": [
    "cd\n",
    "ls -1 $HOME/tutorials/trimmed_reads"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22a79a81-00a5-4cc9-9a72-dab9f7c26fd4",
   "metadata": {},
   "source": [
    "## De-novo assembly with Shovill\n",
    "\n",
    "Shovill is a tool that optimizes the assembler `Spades` to minimize the run time, while maintaining the quality of results. It generates a draft genome using heuristic algorithms and does not require a reference genome to guide the process. See the GitHub repositories of [shovill](https://github.com/tseemann/shovill) and [SPAdes](https://github.com/ablab/spades) for more details. \n",
    "\n",
    "Shovill is not available as a module pre-installed in **CC**, so we must use another strategy. The easiest one is to use a container, we can install a **Docker** container, but Docker containers are not suitable for high performance clusters like **Compute Canada** because they are used mostly when you have administrator privileges. Thus, many HPC allow use of **Singularity** images as an alternative (For more info about what is containerization you can read https://www.melbournebioinformatics.org.au/tutorials/tutorials/docker/docker/).\n",
    "\n",
    "A useful repository of **Singularity** images is located at https://depot.galaxyproject.org/singularity/\n",
    "\n",
    "<font color='darkred'>_**Notes for compute canada:**_ </font>  \n",
    "- Singularity needs to be loaded into the system. On the CIDGOH servers, it is loaded by default. Run the following code to have singularity available in your compute canada session. \n",
    "\n",
    ">    module load singularity\n",
    "\n",
    "- We have already downloaded the **Shovill** singularity container. In CC, you may need to do it, so run the following command pull it from the repository into your local directory. The command tells the system to pull a container from a repository into your local system.\n",
    "\n",
    ">    singularity pull shovill_1.1.sif https://depot.galaxyproject.org/singularity/shovill%3A1.1.0--hdfd78af_1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0295f341-5926-4e79-91c9-d98631079b29",
   "metadata": {
    "tags": [],
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [],
   "source": [
    "# executing shovill\n",
    "\n",
    "singularity exec /mnt/cidgoh-object-storage/images/shovill_1.1.sif shovill --help"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2ba8b22-fddf-4774-b9d0-5c0fb23a33ae",
   "metadata": {},
   "source": [
    "### What does Shovill do?\n",
    "\n",
    "1. Unifies coverage depth (how many times is a region covered by reads on average) for all genomes\n",
    "2. Trims adapters and poor quality reads if necessary\n",
    "3. Assembles using SPAdes\n",
    "4. Polishes genomes (improves quality) and filters low quality contigs\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b145e909-c993-46e8-bf8f-ac42000e5bf2",
   "metadata": {
    "tags": [],
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [],
   "source": [
    "# create output directory\n",
    "CONTIGS_DIR=\"$HOME/tutorials/results/contigs\"\n",
    "\n",
    "# define PATH to trimmed_reads\n",
    "TRIMMED_READS=\"$HOME/tutorials/trimmed_reads\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "858935a0-b97d-4969-aeae-88b3885a5df4",
   "metadata": {},
   "source": [
    "To execute shovill, we run the command from the singularity container we just downloaded. Genome assembly is the most resource intensive process in the pipeline, so it will probably take a while to run (20 - 30 min). As input, we will use the `trimmed_reads` files and only two isolates. The remaining ones are already assembled and can be copied into our environment later.\n",
    "\n",
    "<font color='darkred'>_**Notes for compute canada:**_ </font>  \n",
    "- Allocate sufficiente memory as the size of every genome must be kept in storage while it is assembled\n",
    "- Bioinformatic procedures usually use multiple threads to optimize performance, so their efficiency increases with the number of available cores (including **SPAdes**). \n",
    "- In shovill, the `--ram` option specifies the available ram per thread (core)\n",
    "    - Spades will take input of RAM from shovill as total available mem, better to input limit manually with `--opts \"-m XX\"`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4eacd90-b635-49c7-b907-0729c900b3d5",
   "metadata": {
    "scrolled": true,
    "tags": [],
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [],
   "source": [
    "# for loop to run a command for each sample\n",
    "\n",
    "for READ1 in $(ls $TRIMMED_READS/*R1.fastq.gz | head -n 2)\n",
    "do\n",
    "\n",
    "    \n",
    "    READ2=${READ1/_R1/_R2}                                                              # substitute R1 for R2 in variable\n",
    "    PREFIX_ISOLATE=$(basename $READ1 _R1.fastq.gz)                                      # create variable with isolate name only\n",
    "    echo \"Started processing $PREFIX_ISOLATE\"\n",
    "    \n",
    "    singularity exec /mnt/cidgoh-object-storage/images/shovill_1.1.sif shovill  \\\n",
    "        --R1 $READ1                                                                     `# specify paired read 1` \\\n",
    "        --R2 $READ2                                                                     `# specify paired read 2` \\\n",
    "        --outdir $CONTIGS_DIR                                                           `# define output directory` \\\n",
    "        --force                                                                         `# overwrite results if already available` \\\n",
    "        --ram 140                                                                       `# how much ram memory to use`\n",
    "    \n",
    "    mv \"$CONTIGS_DIR/contigs.fa\" $CONTIGS_DIR/${PREFIX_ISOLATE}_contigs.fa\n",
    "    \n",
    "    echo \"Finished assembly of sample $PREFIX_ISOLATE\"\n",
    "    \n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38b71de2-4acc-47c9-bcf2-9b3892fd47d5",
   "metadata": {},
   "source": [
    "In this tutorial, we optimized the runtime by processing only two samples. The remaining assemblies can be imported in the same directory for future steps using the command below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85edd2c4-8e55-44aa-8c1f-254f8a95784d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cp /mnt/cidgoh-object-storage/seagull/jupyter-mdprieto/precomputed_results/contigs/* $HOME/tutorials/results/contigs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "943e0ddc-c0ed-4de1-b6e7-60550e1a5586",
   "metadata": {},
   "source": [
    "\n",
    "The main output of the **Shovill** pipeline are the files ending in `contigs.fa` which contain assembled reads in fasta format. We can verify that this format contains a header for every contig and then the reads.  \n",
    "\n",
    "More information about the `.fasta` format can be found in this [wikipedia page.](https://en.wikipedia.org/wiki/FASTA_format)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c36552e-64da-4711-b61e-d75cb95e3155",
   "metadata": {
    "tags": [],
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [],
   "source": [
    "head \"$CONTIGS_DIR/ERR10479510_contigs.fa\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9041bbab",
   "metadata": {},
   "source": [
    "As a refresher, this is the current structure of directories for our `tutorials` project folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4514beb",
   "metadata": {
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [],
   "source": [
    "tree -dL 2 $HOME/tutorials"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d38f655-b711-477f-8792-2c1c7d0e30b3",
   "metadata": {},
   "source": [
    "## Quality control of draft genomes\n",
    "\n",
    "**Shovill** produces contigs (overlapping consensus regions of DNA) for every isolate. However, we may have contaminated cultures growing other bacteria besides our organism of interest. Also, given the non-targeted approach of sequencing, the reads from an isolate may have poor quality (low reliability in base calling or poor coverage of certain regions) resulting in inadequate assemblies.\n",
    "\n",
    "Thus, after producing draft genomes, we typically conduct additional checks to verify the characteristics of the resulting files and make sure that we do not have contamination in our samples. \n",
    "\n",
    "### General expectation\n",
    "\n",
    "- The average size of contigs is over 5000 base pairs \n",
    "- No more than 1,000 should be produced in total \n",
    "- The resulting assembly should have a coverage of at least 90% of the reference genome"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f675c8dc-2ca1-4f16-bfb3-d21cdb11e36a",
   "metadata": {},
   "source": [
    "### Quast\n",
    "\n",
    "Quast [(github:quast)](https://github.com/ablab/quast) produces quantitative summaries of the contigs in every assembly. If we specify a reference genome for the organism, it can evaluate misassemblies, unaligned contigs, and overall coverage against the reference genome. \n",
    "\n",
    "**_Some metrics include:_** \n",
    "- Number of contigs and number of contigs > 500bp\n",
    "- **N50** or the length at which the collection of all contigs of at least that length covers half of the assembly \n",
    "- **NG50** is similar to **N50** but measures the coverage of the reference genome\n",
    "- Number of misassemblies including inversions, relocations, and translocations\n",
    "- Number and total length of unaligned contigs (against the reference genome)\n",
    "\n",
    "As the data we are analyzing in the tutorial comes from _P. aeruginosa_ isolates ([PMID:34412676](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8376114/) - [BioProject:PRJEB56397](https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJEB56397)), we will use the [PAO1 reference strain](https://www.ncbi.nlm.nih.gov/data-hub/genome/GCF_000006765.1/) for quality control. \n",
    "\n",
    "Reference genomes contain two files, a fasta sequence file [`.fna, .fa`] and an annotation file [`.gff`]. These reference files are stored in the software directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7cbe0a00-80a9-4853-9165-c6435d69dab0",
   "metadata": {
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [],
   "source": [
    "# load quast into our environment\n",
    "module load StdEnv/2020 gcc/9.3.0 quast/5.0.2\n",
    "\n",
    "echo -e \"Remember the structure of our software directory:\n",
    "\n",
    "$(tree -dL 2 /mnt/cidgoh-object-storage/seagull/jupyter-mdprieto)\n",
    "\n",
    "And the contents of the reference_data subfolder are:\n",
    "\n",
    "$(cd /mnt/cidgoh-object-storage/seagull/jupyter-mdprieto/reference_data/ && ls)\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9574d47-dd93-4f4f-bb8b-abd63697b7fd",
   "metadata": {
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [],
   "source": [
    "# create ENV variable to input and output directory\n",
    "CONTIGS_DIR=\"$HOME/tutorials/results/contigs\"\n",
    "RESULTS_QUAST=\"/$HOME/tutorials/results/assembly_quast\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee0858fb-25ae-4a69-8a9b-9042dafe56fa",
   "metadata": {},
   "source": [
    "The main command `quast.py` produces several reports, with formats such as `.pdf, .html, and .csv`, containing the previously mentioned metrics. \n",
    "\n",
    "They can be opened directly in Jupyter by clicking the file on the explorer or exported to your local computer for further visualization. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd20360a-373c-4e86-94fe-2a353da98649",
   "metadata": {
    "scrolled": true,
    "tags": [],
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [],
   "source": [
    "quast.py $CONTIGS_DIR/*contigs.fa                                                                                       `# pattern for contig files produced by shovill` \\\n",
    "    -r /mnt/cidgoh-object-storage/seagull/jupyter-mdprieto/reference_data/GCF_000006765.1_ASM676v1_genomic.fna.gz       `# reference genome` \\\n",
    "    -g /mnt/cidgoh-object-storage/seagull/jupyter-mdprieto/reference_data/GCF_000006765.1_ASM676v1_genomic.gff.gz                                     `# reference genomic features positions` \\\n",
    "    -o $RESULTS_QUAST                                                                                                  `# output directory` \\\n",
    "    --threads 12"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f801be8-1205-4306-b03a-6f6d7e29c5aa",
   "metadata": {},
   "source": [
    "We can preview metrics such as the N50 of the assemblies and the coverage of the reference genome using the commands below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ca9fb58-07ab-4900-a529-03a3e1500e25",
   "metadata": {
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [],
   "source": [
    "echo -e \"These files are the output of QUAST:\"\n",
    "ls -1 $RESULTS_QUAST\n",
    "\n",
    "# overview of N50 and coverage\n",
    "echo\n",
    "cat $RESULTS_QUAST/report.tsv                      `# print the file` | \\\n",
    "cut -f 1-5                                         `# cut the columns 1 to 5 (separated by tab)` | \\\n",
    "grep -E 'Assembly|N50|fraction'                    `# select lines matching the pattern Assembly OR N50 OR fraction` | \\\n",
    "column -ts $'\\t'                                   `# print table separating at tabs ($'\\t')`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf0a1ad4-5abf-48b9-bea4-3d496ff19bb8",
   "metadata": {
    "tags": []
   },
   "source": [
    "### CheckM\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17902515-d95e-454f-845f-67491e1f046c",
   "metadata": {},
   "source": [
    "**CheckM** infers the quality of the genome assembly based on the presence and uniqueness of sets of gene markers that are specific to species/taxa, and determines the completeness (coverage of reference genome) and the contamination of the input draft genomes.\n",
    "\n",
    "**CheckM** is not available in the CC cluster. So, we use a singularity container with the latest version. For more info, visit https://github.com/Ecogenomics/CheckM. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9e5158e-498e-4b56-93d3-8572efd23fa6",
   "metadata": {
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [],
   "source": [
    "# create output directory and set environment variable \n",
    "RESULTS_CHECKM=\"$HOME/tutorials/results/assembly_checkm\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0107f31e-d8ac-42ba-bf71-8c84b3f6921c",
   "metadata": {},
   "source": [
    "1. The following command creates a dataset with specific genomic markers for a species, taxon or genus \n",
    "\n",
    "`checkm taxon_set <species/genus/taxon> <taxon_name> <marker_file>`\n",
    "\n",
    "**Note:** some of the files used in the tutorials are in the `/mnt` shared drive/folder. Singularity may not be able to recognize that drive by default, so we run `export SINGULARITY_BIND=\"/mnt,/etc\"` to make sure its readable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3bcd0c5d-9d12-4aa1-8283-44bd885667fa",
   "metadata": {
    "scrolled": true,
    "tags": [],
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [],
   "source": [
    "# add drives to singularity\n",
    "export SINGULARITY_BIND=\"/mnt,/etc\"\n",
    "\n",
    "# command\n",
    "singularity exec /mnt/cidgoh-object-storage/images/checkm_1.2.2.sif checkm \\\n",
    "    taxon_set species 'Pseudomonas aeruginosa' /mnt/cidgoh-object-storage/seagull/jupyter-mdprieto/reference_data/pseudomonas.ms"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "460af426-1022-4e98-a6ce-a1c1c91ffa2b",
   "metadata": {},
   "source": [
    "\n",
    "2. Using `checkm analysis` we identify what taxon-specific marker sets are included in every assembly. The process for the samples used in the tutorial (n = 10) should take around 5 min. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a1728b0-af61-4b2c-97ac-9b1e8ccd849a",
   "metadata": {
    "scrolled": true,
    "tags": [],
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [],
   "source": [
    "singularity exec /mnt/cidgoh-object-storage/images/checkm_1.2.2.sif checkm analyze \\\n",
    "    /mnt/cidgoh-object-storage/seagull/jupyter-mdprieto/reference_data/pseudomonas.ms           `#file with checkm marker set for assemblies` \\\n",
    "    $CONTIGS_DIR                                                                                `#dir with assemblies in fasta format` \\\n",
    "    $RESULTS_CHECKM                                                                             `#output directory` \\\n",
    "    -x fa                                                                                       `#extension of assemblies` \\\n",
    "    -t 10 "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e866472e",
   "metadata": {},
   "source": [
    "2. Then, with `checkm qa` we produce a summary of contamination results. This step should not take long (less than 1 min)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90fee905",
   "metadata": {
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [],
   "source": [
    "singularity exec /mnt/cidgoh-object-storage/images/checkm_1.2.2.sif checkm qa                                                                 \\\n",
    "        /mnt/cidgoh-object-storage/seagull/jupyter-mdprieto/reference_data/pseudomonas.ms       `#file with checkm marker set for assemblies` \\\n",
    "        $RESULTS_CHECKM                                                                        `#output directory` \\\n",
    "        --file $RESULTS_CHECKM/checkm_output.tsv                                                                    \\\n",
    "        --tab_table                                                                             `# print tabular output` \\\n",
    "        --threads 10                                                                            `# number of simultaneous threads for process` \\\n",
    "        --out_format 1                                                                          `# format of output 1 = summary, 2 = extended`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bcd13edd-e924-4e26-870e-3bfcf856051f",
   "metadata": {
    "jupyter": {
     "outputs_hidden": true
    },
    "tags": []
   },
   "source": [
    "Now, we can review the output of results for contamination using **CheckM**. \n",
    "\n",
    "- *Completeness* is a measure of the coverage of gene marker sets spected for a species in a given contig. \n",
    "- *Contamination* shows the presence of multi-copy marker genes in the genome assembly. \n",
    "- *Strain heterogeneity* is determined by the number of multy-copy gene marker sets that have an amino acid identity >=  90%.\n",
    "\n",
    "A high heterogeneity suggests that a majority of the contamination comes from closely related organisms. A smaller value may come from phylogenetically distinct sources"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ac32a13-2121-4bcd-a12b-0265958f3e76",
   "metadata": {
    "vscode": {
     "languageId": "shellscript"
    }
   },
   "outputs": [],
   "source": [
    "column -ts $'\\t' $RESULTS_CHECKM/checkm_output.tsv  "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Bash",
   "language": "bash",
   "name": "bash"
  },
  "language_info": {
   "codemirror_mode": "shell",
   "file_extension": ".sh",
   "mimetype": "text/x-sh",
   "name": "bash"
  },
  "toc-autonumbering": true,
  "toc-showmarkdowntxt": false,
  "vscode": {
   "interpreter": {
    "hash": "f6693f812c7ab8e381cab3b9fa9a750f08ad5ffac90c3671106fa5f1758428d8"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}