08.Module_8_serotyping.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Module 8: Serotyping of *Streptococcus pneumoniae* and *Streptococcus agalactiae* \n",
    "\n",
    "## Overview\n",
    "\n",
    "# This module will be developed in two parts:\n",
    "\n",
    "1. Part 1: Predicting serotypes of *S. pneumoniae*\n",
    "\n",
    "2. Part 2: Predicting serotypes of *S. agalactiae*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Install condacolab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install -q condacolab\n",
    "import condacolab\n",
    "condacolab.install()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Install software"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Note**: In this module, we will use SeroBA for the prediction of *S. pneumoniae* serotypes. However, the original repository for this tool is no longer supported. Therefore, we will use the [Bentley-group](https://github.com/sanger-bentley-group/seroba) repository, where the use of the tool with Docker is recommended.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What is Docker?\n",
    "\n",
    "Docker is an open platform for developing, shipping, and running applications. It provides the ability to package and run an application in a loosely isolated environment called a container. Containers are lightweight and contain everything needed to run the application, so you do not need to rely on what is currently installed on the host; i.e., it involves bundling an application together with all of the necessary configuration files, libraries, and dependencies to ensure the software can run in a reproducible fashion across a diversity of computing environments. You can easily share containers while you work, and be sure that everyone you share with gets the same container that works in the same way.\n",
    "\n",
    "![docker](images/docker.png)\n",
    "\n",
    "*Taken from: https://docs.docker.com/get-started/docker-overview/*\n",
    "\n",
    "If you want to delve deeper into how Docker and its containers work, you can visit: https://docs.docker.com/get-started/introduction/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Install udocker\n",
    "%%shell\n",
    "!pip install udocker\n",
    "!udocker --allow-root install"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get SeroBA\n",
    "!udocker --allow-root pull sangerbentleygroup/seroba"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Chek if SeroBA is installed\n",
    "!udocker --allow-root run sangerbentleygroup/seroba seroba version"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Note**: In this module, we will use srst2 for the prediction of *S. agalactiae* serotypes. However, this tool does not work properly with Python3. Therefore, we will install Python2 and run the tool using the python=2.7 version with the commands `!conda run -n py2_env`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a Python 2.7 environment with conda. The environment is called py2_env\n",
    "!conda create -n py2_env python=2.7\n",
    "#!apt-get install python2\n",
    "!conda run -n py2_env python --version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Install srst2\n",
    "!conda run -n py2_env conda install -c bioconda srst2 --yes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check if srst2 is installed\n",
    "!conda run -n py2_env srst2 --help"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://zenodo.org/records/14231070/files/Module_8.tar.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extract the .tar.gz file "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!tar xvf Module_8.tar.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1: Predicting *S. pneumoniae* serotypes\n",
    "\n",
    "To date, there are >100 known serotypes described for *S. pneumoniae* based on differing biochemical and antigenic properties of the capsule. There are a number of in-silico methods to detect the cps locus, which can then be used to predict serotypes from WGS data. Accurate identification of pneumococcal serotypes is important for tracking the distribution and evolution of serotypes following the introduction of effective vaccines.\n",
    "\n",
    "[SeroBA](https://github.com/sanger-pathogens/seroba?tab=readme-ov-file#installation) was developed and it makes efficient use of computational resources in addition to accurately detecting the cps locus at low coverage, and it predict serotypes from WGS data using a database adapted from PneumoCaT. SeroBA can predict serotypes, by identifying the cps locus, directly from raw whole genome sequencing read data with 98% concordance using a k-mer based method, can process 10,000 samples in just over 1 day using a standard server and can call serotypes at a coverage as low as 10x. SeroBA is implemented in Python3 and is freely available under an open source GPLv3.\n",
    "\n",
    "*Further reading*: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6113868/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Explore usage of serobar by running ``seroba -h``"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run seroBA\n",
    "!udocker --allow-root run sangerbentleygroup/seroba seroba --help"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SeroBA requires only three inputs: \n",
    "\n",
    "1. Database with kmc (utility designed for counting k-mers) and ariba (Antimicrobial Resistance Identification By Assembly). \n",
    "\n",
    "For this module, we will use the database `seroba_k71_14082017` located in the module folder.\n",
    "\n",
    "If you want to create the database yourself on your computer, you can follow the steps provided [here](https://github.com/sanger-pathogens/seroba?tab=readme-ov-file#setting-up-the-database).\n",
    "\n",
    "2. Forward and reverse sequence files in fastq\n",
    "\n",
    "We will use the paired files ERR331173_1.fastq.gz and ERR331173_2.fastq.gz, which correspond to run accession ERR331173 from the [PRJEB3084](https://www.ebi.ac.uk/ena/browser/view/PRJEB3084) project.\n",
    "\n",
    "Some important data about the sample:\n",
    "\n",
    "- Country of origin: Perú\n",
    "- Organism: *Streptococcus pneumoniae*\n",
    "- Instrument Platform: ILLUMINA\n",
    "- Instrument Model: Illumina MiSeq\n",
    "- Read Count: 1203646\n",
    "- Base Count: 240729200\n",
    "- Center Name: Wellcome Sanger Institute; SC\n",
    "- Library Layaout: PAIRED\n",
    "- Library strategy: WGS\n",
    "\n",
    "\n",
    "3. Output prefix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### To predict the serotype of a single strain, we will use the command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run seroBA\n",
    "!udocker --allow-root run -v /content/Module_8/seroba_k71_14082017:/seroba_k71_14082017 -v /content/Module_8:/fastq_files sangerbentleygroup/seroba seroba runSerotyping /seroba_k71_14082017 /fastq_files/ERR1795461_1.fastq.gz /fastq_files/ERR1795461_2.fastq.gz /fastq_files/output_test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An explanation of this command is as follows:\n",
    "\n",
    "`udocker` is a tool that allows running Docker containers without needing administrator (root) privileges. It is designed for environments where Docker cannot be installed directly, such as shared systems or cloud computing platforms where users do not have root permissions.\n",
    "\n",
    "`--allow-root` This option allows udocker to run as a superuser (root). Normally, udocker is designed to function without elevated privileges, but this option allows it to be run by the root user.\n",
    "\n",
    "`-v` This option is used to mount volumes between the host and the container. This makes files in the host's folder accessible inside the container, facilitating data exchange between the host and the container.\n",
    "\n",
    "`run` is a function to start a container.\n",
    "\n",
    "`sangerbentleygroup/seroba` is the Docker image [sangerbentleygroup represents the repository and seroba represents the container image].\n",
    "\n",
    "`seroba` is the tool/software\n",
    "\n",
    "`runSerotyping` specifies that program will perform serotyping \n",
    "\n",
    "`seroba_k71_14082017` specifies where the seroba directory\n",
    "\n",
    "`ERR331173_1.fastq.gz` and `ERR331173_2.fastq.gz` are the forward and reverse fastq files\n",
    "\n",
    "`ERR331173_output` specifies the output prefix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the output folder,  you will find a **pred.tsv** including your predicted serotype."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### To predict the serotype of multiple strains\n",
    "\n",
    ">**Note**: In this module, we will not run the multiple serotyping due to the lack of resources in Colab. However, here is an example of how to do it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First, we will create a folder for each pair of compressed FASTQ files and the strain identification name using the command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Do not execute\n",
    "# Move files\n",
    "# for x in *1.fastq.gz; do mkdir ${x%%_1.fastq.gz} ; mv $x ${x%%_1.fastq.gz}; mv ${x%%1.fastq.gz}2.fastq.gz ${x%%_1.fastq.gz}; done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An explanation of this command is as follows:\n",
    "\n",
    "`for x en 1.fastq.gz;` This starts a for loop that iterates over all files in the current directory that end with \"1.fastq.gz\".\n",
    "\n",
    "`do` This starts the code block that will be executed for each file.\n",
    "\n",
    "`mkdir ${x%%_1.fastq.gz}` This creates a new directory with the same name as the file, but with the \"_1.fastq.gz\" suffix removed.\n",
    "\n",
    "`mv $x ${x%%_1.fastq.gz}` This moves the original file into the new directory created in the previous step.\n",
    "\n",
    "`mv ${x%%1.fastq.gz}2.fastq.gz ${x%%_1.fastq.gz}` This moves a second file that has the same prefix as the first file, but with a \"2.fastq.gz\" suffix, into the new directory created in the first step.\n",
    "\n",
    "`done` This ends the for loop.\n",
    "\n",
    "Overall, this script is designed to take paired-end sequencing data that is stored in two separate files with names that end in \"_1.fastq.gz\" and \"_2.fastq.gz\", and organize it into directories based on the prefix of the file names. \n",
    "\n",
    "*Further reading on loop commands*: https://www.gnu.org/software/bash/manual/bash.html#Looping-Constructs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### we will then run seroBA using the command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Do not execute\n",
    "# Run seroBA for all samples\n",
    "#for x in * ; do seroba runSerotyping seroba_k71_14082017 $x/${x}_1.fastq.gz $x/${x}_2.fastq.gz $x\"_output\"; done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An explanation of this command is as follows:\n",
    "\n",
    "`for x in *;` This starts a for loop that iterates over all files in the current directory.\n",
    "\n",
    "`do` This starts the code block that will be executed for each file.\n",
    "\n",
    "`seroba runSerotyping seroba_k71_14082017` This executes the command \"runSerotyping\" within the Docker container, using the \"seroba_k71_14082017\" database. \n",
    "\n",
    "`$x/${x}_1.fastq.gz $x/${x}_2.fastq.gz` These are the input files for the command, located in a directory named after the file (with the `#` character removed) and with \"_1.fastq.gz\" or \"_2.fastq.gz\" appended to the name.\n",
    "\n",
    "`$x\"_output\"` This is the output directory for the command, also named after the file (with \"_output\" appended to the name).\n",
    "\n",
    "`;` This separates the Docker command from the end of the for loop.\n",
    "\n",
    "`done` This ends the for loop."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### We will then compile the results from the runs above using the command:\n",
    "\n",
    "This command will combine the seroba outputs in one tsv file. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Do not execute\n",
    "# Summary of seroBA results\n",
    "#!seroba summary ./"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: Predicting *S. agalactiae* serotypes \n",
    "\n",
    "*Streptococcus agalactiae* (Group B Streptococcus, or GBS) are currently divided into ten serotypes based on type-specific capsular antigens and are designated as Ia, Ib, II, III, IV, V, VI, VII, VIII, and IXs. \n",
    "\n",
    "Group B Streptococcus Serotyping by Genome Sequencing repository contains a curated reference file which can be used for serotyping *Streptococcus agalactiae* in silico with whole genome sequencing data. The reference file (GBS-SBG.fasta) is designed to be usable for both short-read mapping and assembly-based strategies.\n",
    "\n",
    "*Further reading*: https://github.com/swainechen/GBS-SBG "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The program [SRST2](https://github.com/katholt/srst2) (Short Read Sequence Typing for Bacterial Pathogens) is designed to take Illumina sequence data, an MLST database, and/or a gene sequence database (e.g., resistance genes, virulence genes, etc.) and report the presence of STs and/or reference genes. This program performs rapid and accurate detection of genes and alleles from short-read WGS sequencing. SRST2 can type reads using any sequence database and can calculate combinatorial sequence types defined in MLST-type databases.\n",
    "\n",
    "*Additional Reading*: https://genomemedicine.biomedcentral.com/articles/10.1186/s13073-014-0090-6\n",
    "\n",
    "### SRST2 requires:\n",
    "\n",
    "1. Paired reads:\n",
    "\n",
    "In this section, we will analyze data from the study [Near-term pregnant women in the Dominican Republic experience high rates of Group B Streptococcus rectovaginal colonization with virulent strains](https://doi.org/10.1371/journal.pgph.0002281) conducted by Laycock KM, Acosta F, Valera S, Villegas A, Mejia E, Mateo C, et al. (2023). We will use the paired FASTQ files SRR23874884_1.fastq.gz and SRR23874884_2.fastq.gz, which correspond to run accession SRR23874884 from the [PRJNA945321](https://www.ebi.ac.uk/ena/browser/view/PRJNA945321) project.\n",
    "\n",
    "Some important data about the sample:\n",
    "\n",
    "- Country of origin: República Dominicana\n",
    "- Organism: *Streptococcus agalactiae*\n",
    "- Instrument Platform: ILLUMINA\n",
    "- Instrument Model: Illumina NovaSeq 6000\n",
    "- Read Count: 1916182\n",
    "- Base Count: 453475786869647750\n",
    "- Center Name: Wellcome Sanger Institute; SC\n",
    "- Library Layaout: PAIRED\n",
    "- Library strategy: WGS\n",
    "\n",
    "2. A database with the sequences in fasta file:\n",
    "\n",
    "GBS-SBG.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### To predict the serotype of a single strain, we will use the command: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, we need to move to the directory where the files are located\n",
    "%cd /content/Module_8/s_agalactiae/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run srst2\n",
    "!conda run -n py2_env srst2 --input_pe SRR23874884_1.fastq.gz SRR23874884_2.fastq.gz --output SRR23874884 --log --gene_db GBS-SBG.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An explanation of this command is as follows:\n",
    "\n",
    "`conda` is the command-line tool for managing environments and packages in Conda, which is a package and environment management system.\n",
    "\n",
    "`run` This Conda subcommand is used to execute commands within a specific Conda environment without needing to activate the environment first.\n",
    "\n",
    "`-n py2_env`: The `-n` option specifies the name of the Conda environment in which the commands should be executed. In this case, `py2_env` is the name of the environment we have created.\n",
    "\n",
    "`srst2` is the tool\n",
    "\n",
    "`--input_pe` specifies the input file are paired end reads which are SRR23874884_1.fastq.gz and SRR23874884_2.fastq.gz \n",
    "\n",
    "`--output`specifies the output file SRR23874884\n",
    "\n",
    "`--log` switch on logging to file, rather than standard output\n",
    "\n",
    "`--gene_db` specifies the database GBS-SBG.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the command `ls -lh` to check the contents in the folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# List files\n",
    "!ls -lh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output file from the above run is `SRR23874884__genes__GBS-SBG__results.txt`\n",
    "\n",
    "\n",
    "So, `cat` `SRR23874884__genes__GBS-SBG__results.txt` to view the contents of this file "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show results\n",
    "!cat SRR23874884__genes__GBS-SBG__results.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### To execute SRST2 on multiple strains, run the command:\n",
    "\n",
    ">**Note**: In this module, we will not run the multiple serotyping due to the lack of resources in Colab. However, here is an example of how to do it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Do not execute\n",
    "# Run srst2 for all samples\n",
    "#!conda run -n py2_env srst2 --input_pe *.fastq.gz --output s.agalactiae --log --gene_db GBS-SBG.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`--input_pe .fastq.gz` specifies the input file are multiple compressed fastq.gz files. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## BONUS!\n",
    "\n",
    "If you are working with BASH in your computer or in a HPC and you have too many files you can optimize commands, loops are very useful for large datasets.\n",
    "\n",
    "Here's a way to do it. \n",
    "\n",
    "Create a new bash script using nano named `serotype.sh`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/bin/bash\n",
    "#Nombre de archivo: Serotype.sh\n",
    "#Este script para serotipificar S.pneumoniae a partir de fastq.gz recortados\n",
    "\n",
    "function docker_run() { docker run --rm=True -u $(id -u):$(id -g) -v $(pwd):/data \"$@\" ;}\n",
    "wordir=/path/to/your/directory/\n",
    "cd $wordir\n",
    "\n",
    "for i in $(ls *_1.trimmed.fastq.gz); do\n",
    "\n",
    "NAME=$(basename $i _1.trimmed.fastq.gz)\n",
    "echo \"$NAME\"\n",
    "j=\"${NAME}_1.trimmed.fastq.gz\"\n",
    "echo \"$j\"\n",
    "k=\"${NAME}_2.trimmed.fastq.gz\"\n",
    "echo \"$k\"\n",
    "\n",
    "docker_run staphb/seroba seroba runSerotyping seroba_k71_14082017 ./$j ./$k ${NAME}_serotype_output;\n",
    "\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We move all the outputs to a new folder `mv *_serotype_output serotype_results`, then compile all the data within the new directory with the command `seroba summary ./`, you should obtain a tsv file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a new bash script using nano named `serotype_2.sh`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/bin/bash\n",
    "#Este script para serotipificar múltiples lecturas genómicas de S. agalactaeae\n",
    "\n",
    "function docker_run() { docker run --rm=True -u $(id -u):$(id -g) -v $(pwd):/data \"$@\" ;}\n",
    "wordir=/path/to/your/directory/\n",
    "cd $wordir\n",
    "mkdir -p serotyping_output\n",
    "\n",
    "for i in $(ls *_1.fastq.gz); do\n",
    "\n",
    "NAME=$(basename $i _1.fastq.gz)\n",
    "echo \"$NAME\"\n",
    "j=\"${NAME}_1.fastq.gz\"\n",
    "echo \"$j\"\n",
    "k=\"${NAME}_2.fastq.gz\"\n",
    "echo \"$k\"\n",
    "\n",
    "docker_run staphb/srst2 srst2 --input_pe $j $k --output ./serotyping_out/${NAME}_output --log --gene_db analysis/clean_data/GBS-SBG.fasta;\n",
    "\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Adapted from:*\n",
    "\n",
    "- Advanced Bioinformatics Course developed for the GPS and JUNO projects - Wellcome Sanger Insitute\n",
    "\n",
    "*Modified by Luisa Sacristán (Universidad de los Andes-CABANA)*"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}