13.Modulo_13_filogenetica.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Módulo 13: Filogenética\n",
    "\n",
    "## Descripción general\n",
    "\n",
    "La filogenética es el estudio de las relaciones evolutivas entre entidades biológicas, a menudo especies, individuos o genes (que pueden denominarse taxones). Los principales elementos de la filogenética se resumen en la siguiente figura.\n",
    "\n",
    "![Intro](images/phylo1.png)\n",
    "\n",
    "*Tomado de: https://www.ebi.ac.uk/training/online/courses/introduction-to-phylogenetics/what-is-phylogenetics/*\n",
    "\n",
    "\n",
    "Los árboles filogenéticos basados en datos del genoma completo nos informan de las relaciones entre aislados bacterianos a una escala muy fina. Cuando combinamos esa información de alta resolución sobre las relaciones evolutivas de los aislados con datos geográficos, podemos comprender mejor la distribución actual del patógeno e inferir los procesos epidemiológicos que han actuado sobre la bacteria a lo largo del tiempo. El ejemplo más sencillo sería que una filogenia mostrara que un patógeno está limitado geográficamente (por ejemplo, los aislados de la misma región siempre se agrupan juntos). Esto podría indicar que el patógeno no se está extendiendo rápidamente. Mientras que un patógeno con una filogenia que muestra que los aislados de regiones distantes están probablemente relacionados con los aislados de regiones cercanas, la interpretación es que es probable que el patógeno se propague a través de las fronteras regionales. La referenciación geográfica de los datos genómicos también puede combinarse con información temporal para estudiar el movimiento de los patógenos en el espacio y el tiempo. Esto es más útil cuando se hace en tiempo real y, por tanto, puede ser útil para la detección y el seguimiento de brotes.\n",
    "\n",
    "En este módulo está dividido en tres partes donde abordaremos: \n",
    "1. Construcción de árboles filogenéticos\n",
    "2. Identificación de regiones de recombinación\n",
    "3. Agrupación mediante popPUNK."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instalar condacolab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install -q condacolab\n",
    "import condacolab\n",
    "condacolab.install()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!conda config --add channels defaults\n",
    "!conda config --add channels bioconda\n",
    "!conda config --add channels conda-forge"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instalar programas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instalar FastTree\n",
    "!conda install bioconda::fasttree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instalar snp-sites\n",
    "!conda install snp-sites"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Istalar Gubbins\n",
    "!conda install gubbins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instalar PopPunk\n",
    "!conda install poppunk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Descargar datos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://zenodo.org/records/13883141/files/Module_13.tar.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extraer el archivo .tar.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!tar xvf Module_13.tar.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parte 1: Creando un árbol filogenético usando FastTree\n",
    "\n",
    "### Paso 1: Realice un alineamiento  de sólo SNP utilizando snp-sites\n",
    "\n",
    "Crear una filogenia a partir de secuencias genómicas completas puede ser un proceso muy lento y de gran carga computacional. Podemos acelerarlo utilizando sólo los sitios variables (SNPs). Sin embargo, debemos ser conscientes de que incluir sólo los sitios variables puede afectar a las estimaciones de la tasa evolutiva realizadas por el software filogenético, por lo que debemos tener en cuenta los sitios que eliminamos en nuestro análisis.\n",
    "\n",
    "Para ello utilizaremos snp-sites. Puede ver las opciones de snp-sites con el comando:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd Module_13"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ejecutar snp-sites\n",
    "!snp-sites -h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En primer lugar, elimine todos los sitios invariantes y cree una alineación de secuencias múltiples de sólo SNP. Utilizaremos los resultados de las ejecuciones rápidas descritas en la página anterior. Ejecute el comando:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!snp-sites -o clean.full.SNPs.aln clean.full.aln"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La explicación de este comando es la siguiente:\n",
    "\n",
    "`snp-sites`: es la herramienta/programa\n",
    "\n",
    "`-o clean.full.SNPs.aln`: especifica el archivo de salida\n",
    "\n",
    "`clean.full.aln`: especifica el archivo de entrada - que es una salida de snippy\n",
    "\n",
    "Podemos ver cuántos sitios invariantes se eliminaron (y qué proporción de A, T, G, C tenían) utilizando:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!snp-sites -C clean.full.aln"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Obtendrá la siguiente salida:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Paso 2: Creación de un árbol filogenético a partir de los SNPs utilizando FastTree\n",
    "\n",
    "Puede ver las opciones para fasttree de la siguiente manera:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!FastTree -h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generaremos un árbol filogenético de máxima verosimilitud utilizando este comando:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!FastTree -nt -gtr clean.full.SNPs.aln > clean.full.SNPs.aln.tree"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La explicación de este comando es la siguiente:\n",
    "\n",
    "`FastTree`: es la herramienta/programa\n",
    "\n",
    "`-nt`: especifica que el alineamiento de entrada es de nucleótidos\n",
    "\n",
    "`gtr`: especifica el modelo evolutivo\n",
    "\n",
    "`clean.full.SNPs.aln`: alineamiento de entrada\n",
    "\n",
    "`clean.full.SNPs.aln.tree`: especifica el nombre del árbol de salida\n",
    "\n",
    "Exploremos la salida de nuestro comando anterior utilizando el comando:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls -lh clean.full.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nuestro árbol de máxima verosimilitud se etiqueta clean.full.SNPs.aln.tree. Podemos visualizarlo con figtree o iTOL."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Paso 3: Visualización de un arbol\n",
    "\n",
    "Descargue el archivo clean.full.SNPs.aln.tree en su ordenador en el ícono de \"Archivo\" en la parte izquierda de Colab y diríjase a https://itol.embl.de/ \n",
    "\n",
    "En la parte superior izquierda, seleccione la opción \"Upload\" como se muestra a continuación:\n",
    "\n",
    "![iTOL](images/itol3.png)\n",
    "\n",
    "Posteriormente, seleccione el botón de \"Browse...\", seleccione su archivo .tree y seleccione \"Upload\":\n",
    "\n",
    "![iTOL](images/itol4.png)\n",
    " \n",
    "Finalmente, obtendrá el siguiente arbol:\n",
    "\n",
    "![iTOL](images/itol.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "O utilizando el paquete ggtree en R:\n",
    "\n",
    "```\n",
    "library(ggtree)\n",
    "\n",
    "mltree <- midpoint.root(read.tree(\"clean.full.SNPs.aln.treef\")) \n",
    "\n",
    "ggtree(mltree) + # plot basic tree\n",
    "\n",
    "geom_tiplab(size=3) + # add tip labels\n",
    "\n",
    "geom_treescale() + # add scale bar\n",
    "\n",
    "xlim(0, 0.0005) # set limits so the plot fits nicely on the screen\n",
    "\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Paso 3: Interpretar un árbol filogenético\n",
    "\n",
    "Debemos recordar que el árbol que hemos creado procede de un alineamiento de SNP únicamente.  Recuerde que cuanto más larga es la rama de una cepa, más mutaciones tiene. La cepa de la muestra ERR2667737 y la cepa de referencia son parecidas entre sí, y presentan algunas mutaciones que también se encuentran en la cepa de la muestra ERR2667694. Las cepas de las muestras ERR2667707 y ERR2667708 comparten mutaciones y presentean mutaciones diferentes al resto de cepas."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parte 2: Identificación de regiones de recombinación usando Gubbins\n",
    "\n",
    "Muchas bacterias participan en altas tasas de recombinación homóloga. Esto significa que donan y reciben segmentos de ADN entre sí. En el contexto de un árbol filogenético, en el que comparamos regiones similares y disímiles para determinar el parentesco de los aislados, esto puede resultar problemático y dar lugar a longitudes de rama que reflejen recombinaciones en lugar de divergencia. Esto es muy importante en el caso del *Streptococcus pneumoniae*, que es naturalmente competente, lo que significa que puede captar fácilmente ADN.\n",
    "\n",
    "[Gubbins](https://github.com/nickjcroucher/gubbins/blob/master/docs/gubbins_manual.md) (Genealogies Unbiased By recomBinations In Nucleotide Sequences) es un algoritmo que identifica iterativamente los loci que contienen elevadas densidades de sustituciones de bases, al tiempo que construye una filogenia basada en las mutaciones puntuales putativas fuera de estas regiones. Las simulaciones demuestran que el algoritmo genera reconstrucciones muy precisas bajo modelos realistas de diversificación a corto plazo de las secuencias mediante mutaciones puntuales y recombinación, y puede ejecutarse en alineaciones de muchos cientos de secuencias de genomas bacterianos. Por tanto, no es adecuado para estudiar la recombinación en la diversidad de toda la especie, algo que puede hacerse gen a gen con programas como fastGEAR. En cambio, funciona con muestras de diversidad limitada, que comparten un ancestro común reciente: una cepa o linaje.\n",
    "\n",
    "El archivo de entrada necesario para Gubbins es un alineamiento FASTA de todo el genoma. Cada secuencia debe tener un identificador único y deben evitarse los caracteres especiales. Las secuencias sólo deben utilizar los caracteres ACGT (bases de ADN),N (base desconocida) o - (hueco en la alineación). Si se va a incluir un árbol de partida, éste deberá tener un formato Newick. La alineación se genera más fácilmente mediante el mapeo de secuencias contra una secuencia de referencia. Esto puede hacerse con el conocido programa de mapeo Snippy.\n",
    "\n",
    "Puede ver los comandos Gubbins de la siguiente manera:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!run_gubbins.py -h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ejecutaremos la herramienta gubbins en un alineamiento completo del genoma y no en un alineamiento SNPs. Ahora vamos a ejecutar este comando:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!run_gubbins.py --mar -p output clean.full.aln"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La explicación de este comando es la siguiente:\n",
    "\n",
    "`run_gubbins.py`: es la herramienta/programa\n",
    "\n",
    "`--mar`: Esta opción indica a Gubbins que utilice un algoritmo de Monte Carlo con cadena de Markov (MCMC) para inferir la historia evolutiva de los genomas.\n",
    "\n",
    "`-p output`: Esta opción especifica el nombre del directorio de salida donde se guardarán los resultados del análisis de Gubbins. Como /data dir es el punto de montaje de PWD de su host debido a -v, verá los archivos en el PWD de su host\n",
    "\n",
    "`clean.full.aln`: fichero de entrada\n",
    "\n",
    "Este comando puede tardar unos minutos en ejecutarse.\n",
    "\n",
    "Veamos lo que gubbins ha hecho utilizando el comando (`ls -l output.*`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Una explicación de los archivos de salida:\n",
    "\n",
    "`output.branch_base_reconstruction.embl`: Reconstrucción de sustitución de bases en formato EMBL\n",
    "\n",
    "`output.recombination_predictions.embl`: Predicciones de recombinación en formato EMBL.\n",
    "\n",
    "`salida.predicciones_recombinación.gff`: Predicciones de recombinación en formato GFF\n",
    "\n",
    "`salida.sitios_polimorfos_filtrados.phylip`: Alineación en formato Phylip de los sitios polimórficos filtrados utilizados para generar la filogenia en la iteración final\n",
    "\n",
    "`output.final_tree.tre`: este archivo contiene la filogenia final en formato Newick; las longitudes de rama están en mutaciones puntuales\n",
    "\n",
    "`output.node_labelled.final_tree.tre`: árbol filogenético final en formato Newick pero con etiquetas internas en los nodos; las longitudes de las ramas están en mutaciones puntuales\n",
    "\n",
    "`output.log`: archivo de registro que especifica el software utilizado en cada paso del análisis, con las citas correspondientes\n",
    "\n",
    "`output.per_branch_statistics.csv`: informe por rama de las sustituciones de bases dentro y fuera de los eventos de recombinación\n",
    "\n",
    "`output.summary_of_snp_distribution.vcf`: archivo VCF que resume la distribución de las mutaciones puntuales."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Puedes explorar estos archivos de salida usando el comando `head`. Por ejemplo, \"output.recombination_predictions.gff\" es un archivo GFF que contiene un registro de cada bloque de recombinación identificado, cuántos SNPs contiene y qué muestras están afectadas.\n",
    "\n",
    "EJEMPLO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**output.final_tree.tre** es una filogenia sin regiones de recombinación. Lea en el árbol filtrado gubbins y grafíquelo con ggtree o visualice en Figtree o iTOL como se explicó previamente.\n",
    "\n",
    "Usted tendrá la salida:\n",
    "\n",
    "![iTOL](images/itol2.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "O utilizando el paquete ggtree en R:\n",
    "```\n",
    "!gubbins.tree <- midpoint.root(read.tree(\"gubbins.final_tree.tre\")) \n",
    "\n",
    "ggtree(gubbins.tree) + # plot basic tree\n",
    "\n",
    "geom_tiplab(size=3) + # add tip labels\n",
    "\n",
    "geom_treescale() + # add a scale bar\n",
    "\n",
    "xlim(0, 300) # change the sizing of the plot so it fits nicely on the screen \n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "También podemos visualizar los bloques de recombinación utilizando una herramienta web llamada Phandango. En su navegador, navegue a:  https://jameshadfield.github.io/phandango/#/\n",
    "\n",
    "Necesitarás los siguientes archivos (arrástrelos y suéltelos)::\n",
    "\n",
    "1. output.final_tree.tre\n",
    "2. output.recombination_predictions.gff\n",
    "3. reference.gff (output from Prokka of Reference_sequence_GPSC1.fa)\n",
    "\n",
    "Phandango debería mostrar automáticamente los bloques de recombinación en rojo (ancestral) y azul (específico de una muestra).\n",
    "\n",
    "![rec](images/phylocanvas.png)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parte 3: Agrupación mediante popPUNK\n",
    "\n",
    "[PopPUNK](https://poppunk.readthedocs.io/en/latest/index.html) es una herramienta para agrupar genomas. Nos referimos a los clusters como clusters de k-mer de longitud variable, o VLKCs. Desde el punto de vista biológico, estos grupos suelen representar cepas distintas. Los subgrupos de cepas se denominan linajes.\n",
    "\n",
    "La siguiente figura muestra un resumen de cómo ejecutar popPUNK\n",
    "\n",
    "![punk](images/poppunk_flowchart.png)\n",
    "\n",
    "*Tomado de: https://poppunk.readthedocs.io/en/latest/overview.html*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Base de datos\n",
    "\n",
    "Como estamos trabajando con *Streptococcus pneumoniae*, usaremos la base de [datos de referencia GPS](https://gps-project.cog.sanger.ac.uk/GPS_v9.tar.gz) y la [designación GPS](https://gps-project.cog.sanger.ac.uk/GPS_v9_external_clusters.csv) que utilizaremos para agrupar nuestro genoma. También puede acceder a los genomas de referencia de otras especies bacterianas desde esta base de datos. Si una especie no está incluida en esta base de datos, se recomienda que cree su propia base de datos.\n",
    "\n",
    "La base de datos del genoma de referencia GPS de *Streptococcus pneumoniae* se guarda en su directorio como (GPS_v9) y la designación GPS como (GPS_v9_external_clusters.csv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Descomprimir la base de datos\n",
    "!tar xvf GPS_v9.tar.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Archivo de texto con los detalles de sus muestras\n",
    "\n",
    "Necesita un archivo que enumere los nombres de sus muestras y las rutas a sus datos de secuencia. En este caso es el archivo: poppunk_input.tsv\n",
    "\n",
    "Puede revisarlo con el siguiente comando:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat poppunk_input.tsv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Este archivo de texto contiene los nombres de las muestras y sus datos de secuencia. No tiene encabezamiento, está separado por tabuladores y contiene el nombre de la muestra en la primera columna. Las columnas siguientes pueden contener rutas a datos de lectura ensamblados o sin procesar (el tipo se deducirá automáticamente comprobando la presencia de puntuaciones de calidad)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Agrupar los genomas\n",
    "\n",
    "El comando para agrupar sus genomas es el siguiente:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!poppunk_assign --db GPS_v9 --query poppunk_input.tsv --output poppunk_clusters --external-clustering GPS_v9_external_clusters.csv "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La explicación de este comando es la siguiente:\n",
    "\n",
    "`poppunk_assign`: es la herramienta/programa/script\n",
    "\n",
    "`--db GPS_v9`: especifica la base de datos\n",
    "\n",
    "`--query /data/poppunk_input.tsv`: recibe la entrada del list.txt\n",
    "\n",
    "`--output poppunk_clusters`: especifica el archivo de salida\n",
    "\n",
    "`--external-clustering GPS_v9_external_clusters.csv`: directorio que contiene los clusters GPS para las referencias\n",
    "\n",
    "\n",
    "Al finalizar, se generará una nueva carpeta \"poppunk_clusters\". Navega hasta esta carpeta y explora su contenido.\n",
    "\n",
    "Los archivos de salida:\n",
    "\n",
    "**poppunk_clusters_clusters.csv**: clusters popPUNK con nomenclatura específica del conjunto de datos\n",
    "\n",
    "**poppunk_clusters_external_clusters.csv**: Designaciones de esquemas GPSC v9\n",
    "\n",
    "Podemos explorar \"poppunk_clusters_clusters.csv\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">Nota: Si a una cepa ya se le ha asignado un cluster, por favor cambie el nombre para ejecutar popPUNK (esto es para evitar aplastar la herramienta) A los nuevos clusters se les asigna NA en el archivo _external_clusters.csv ya que no han sido definidos en el conjunto de datos v6 utilizado para designar los GPSCs. Por favor, envíe un correo electrónico a: globalpneumoseq@gmail.com para que se añadan nuevos clusters a la base de datos y se les asigne un nombre de cluster GPSC después de haber comprobado que no hay contaminación de bajo nivel que pueda contribuir a distancias accesorias sesgadas."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Este módulo popPUNK ha sido desarrollado a partir de las siguientes páginas web:\n",
    "\n",
    "https://poppunk.readthedocs.io/en/latest/query_assignment.html\n",
    "\n",
    "https://www.pneumogen.net/gps/training_command_line.html\n"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}