11.Modulo_11_mapeo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Módulo 11: Mapeo\n",
    "\n",
    "## Descripción general\n",
    "\n",
    "El mapeo de lecturas cortas con respecto a un genoma de referencia suele ser el primer paso para analizar los datos de secuenciación de nueva generación, y debe ser lo más preciso posible. Debido al elevado número de lecturas que hay que manejar, se han desarrollado numerosos algoritmos sofisticados para abordar este problema y en la actualidad existen muchas herramientas de mapeo.\n",
    "\n",
    "### Selección de una secuencia de referencia para el mapeo\n",
    "\n",
    "Se designa un único genoma de referencia que represente a una especie para el análisis comparativo. Un genoma de referencia completo debe tener una anotación de alta calidad y reunir el máximo nivel de apoyo experimental para la anotación estructural y funcional.\n",
    "\n",
    ">**Nota**: Es importante tener en cuenta que, aunque un único genoma de referencia arbitrario es un enfoque frecuentemente utilizado en genómica microbiana, la elección de una referencia puede representar una fuente de errores que pueden afectar a los análisis posteriores, como la detección de polimorfismos de nucleótido único (SNP) y la inferencia filogenética.\n",
    "\n",
    "*Lecturas adicionales*: \n",
    "\n",
    "https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3375638/\n",
    "\n",
    "https://academic.oup.com/nar/article/42/D1/D553/1066302 \n",
    "\n",
    "https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1008678 "
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "### Instalar programas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instalar bwa\n",
    "!conda install -c bioconda bwa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Install samtools\n",
    "!conda config --add channels bioconda\n",
    "!conda config --add channels conda-forge\n",
    "!conda install -c bioconda samtools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Descargar datos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://zenodo.org/records/13750987/files/Module_11.tar.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extraer el archivo .tar.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!tar xvf Module_11.tar.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parte 1: Evaluar archivos\n",
    "\n",
    "Navegue hasta esta carpeta del módulo y explore su contenido. Debería ver los archivos de lectura de extremos emparejados FASTQ comprimidos:\n",
    "\n",
    "- ERR2667737_1.fastq.gz\n",
    "\n",
    "- ERR2667737_2.fastq.gz\n",
    "\n",
    "En este módulo analizaremos las secuencias que  corresponden al run accession [ERR2667737](https://www.ebi.ac.uk/ena/browser/view/ERR2667737) del proyecto [PRJEB3084](https://www.ebi.ac.uk/ena/browser/view/PRJEB3084).\t\n",
    "\n",
    "Algunos datos importantes de la muestra:\n",
    "\n",
    "- País de origen: Argentina\n",
    "- Organismo: *Streptococcus pneumoniae*\n",
    "- Instrumento: ILLUMINA\n",
    "- Modelo del instrumento: Illumina HiSeq X Ten  \n",
    "- Conteo de reads: 2148156\n",
    "- Conteo de bases: 324371556 \n",
    "- Nombre del centro: Wellcome Sanger Institute;SC \n",
    "- Diseño de librerías:  PAIRED\n",
    "- Estrategia de librerías: WGS\n",
    "\n",
    "Y también dos archivos FASTA de secuencias de referencia:\n",
    "\n",
    "- Reference_sequence_GPSC1.fa\n",
    "\n",
    "- Reference_sequence_GPSC33.fa\n",
    "\n",
    "Alinearemos las lecturas de extremos emparejados con las dos secuencias de referencia. Las secuencias de referencia representan los linajes GPSC1 y GPSC33 de *Streptococcus pneumoniae*, y utilizaremos los resultados de la alineación para determinar a qué linaje pertenecen las lecturas.\n",
    "\n",
    "Más información: https://pubmed.ncbi.nlm.nih.gov/31003929/\n",
    "\n",
    "Como comprobación rápida, cuente el número de líneas en cada uno de los archivos de lectura y compruebe que tienen el mismo número. Como se trata de lecturas de extremos pareados, debería haber una lectura de cada par de lecturas en cada uno de los archivos y, por lo tanto, el mismo número de líneas (y, por lo tanto, de lecturas) en cada archivo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd Module_11\n",
    "!ls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!gzip -cdf ERR2667737_1.fastq.gz | wc -l"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!gzip -cdf ERR2667737_2.fastq.gz | wc -l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A continuación se ofrece una explicación de estos comandos:\n",
    "\n",
    "`gzip`: Comando para comprimir o descomprimir archivos\n",
    "\n",
    "`cdf`: Opciones donde c significa escribir el archivo en la salida estándar, manteniendo los archivos originales sin cambios; d significa descomprimir (a .fastq); y f significa forzar la sobreescritura de los archivos de salida y el uso de enlaces simbólicos. (En este caso, los archivos fastq.gz son enlaces simbólicos)\n",
    "\n",
    "`|`: Canaliza el archivo .fastq de salida al siguiente comando\n",
    "\n",
    "`wc`: comando para un programa llamado word count, proporcionando el flag `-l` le decimos a word count que cuente el número de líneas en un archivo.\n",
    "\n",
    "Los archivos de lecturas pareadas siempre deben tener el mismo número de líneas/lecturas (el orden de las lecturas en cada archivo también es crítico), por lo que si sus dos archivos emparejados tienen un número diferente de líneas, algo ha ido mal (por ejemplo, el filtrado/recorte ha ido mal y ha corrompido la salida, o quizás se están utilizando archivos de muestras diferentes)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parte 2: Creación de un alineamiento\n",
    "\n",
    "### Paso 1: Indexar el genoma de referencia\n",
    "\n",
    "Primero crearemos un índice del archivo de secuencia de referencia \"Reference_sequence_GPSC1.fa\". Para referencias grandes, esto puede llevar un tiempo, pero una vez creado el índice puede reutilizarlo para múltiples muestras. Por ejemplo, si tiene 100 muestras de *Streptococcus pneumoniae*, sólo tendrá que crear el índice una vez y podrá utilizarlo para las 100 muestras.\n",
    "\n",
    "Ejecute el comando en el terminal para ejecutar bwa:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Indexar el genoma de referencia\n",
    "!bwa index Reference_sequence_GPSC1.fa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La explicación de este comando es la siguiente:\n",
    "\n",
    "`bwa`: es la herramienta\n",
    "\n",
    "`index`: indica a bwa que indexe el archivo de entrada\n",
    "\n",
    "`Reference_sequence_GPSC1.fa`: es el archivo de entrada\n",
    "\n",
    "Si lista `ls` el contenido del directorio, ahora debería ver los archivos de índice bwa, tendrán el prefijo \"Reference_sequence_GPSC1\" y tendrán extensiones como .amb, .ann, .bwt, .pac y .sa."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Listar los archivos en el directorio\n",
    "!ls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Paso 2: Alineamiento\n",
    "\n",
    "Para alinear las lecturas con la secuencia de referencia, escriba este comando:\n",
    "\n",
    ">Nota: Este proceso puede tardar varios minutos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Alinear las secuencias contra el genoma de referencia\n",
    "!bwa mem Reference_sequence_GPSC1.fa ERR2667737_1.fastq.gz ERR2667737_2.fastq.gz > GPSC1bwa.sam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La explicación de este comando es la siguiente:\n",
    "\n",
    "`bwa`: es la herramienta\n",
    "\n",
    "`mem`: indica a bwa que utilice el algoritmo mem para alinear los archivos de lectura con la referencia\n",
    "\n",
    "`ERR2667737_1.fastq.gz`: archivo de entrada de lecturas forward\n",
    "\n",
    "`ERR2667737_1.fastq.gz`: archivo de entrada de lecturas inversas\n",
    "\n",
    "`Reference_sequence_GPSC1.fa`: es el archivo de entrada\n",
    "\n",
    "`> GPSC1bwa.sam`: dirige `>` los resultados de la alineación al fichero GPSC1bwa.sam\n",
    "\n",
    "Cuando bwa se ejecute, imprimirá mensajes en la pantalla de su terminal:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cuando bwa haya terminado, comprueba que el fichero SAM ha sido creado usando `ls`\n",
    "\n",
    "Ahora debería haber un archivo \"GPSC1bwa.sam\" en el directorio. Normalmente, un archivo SAM contiene una sola línea para cada lectura del conjunto de datos, y esta línea almacena los resultados de la alineación de cada lectura (nombre de referencia, ubicación de la alineación, cadena CIGAR, la propia secuencia de lectura, calidad, etc.).\n",
    "\n",
    "Los archivos SAM están en formato de texto. Para visualizarlos, escriba:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Listar los archivos en el directorio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ver las primeras líneas del archivo SAM\n",
    "!head GPSC1bwa.sam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Obtendrá el siguiente resultado:\n",
    "\n",
    "![bwa](images/bwa.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Paso 3: Convertir el archivo SAM a BAM\n",
    "\n",
    "Es una buena práctica convertir los archivos SAM en archivos BAM (Mapa de Alineamiento Binario), que son versiones binarias comprimidas de los mismos datos, y pueden ser ordenados e indexados fácilmente para hacer las búsquedas más rápidas. Utilizaremos samtools para convertir nuestro SAM a BAM, y ordenar e indexar el archivo BAM.\n",
    "\n",
    "Para convertir un archivo sam en un archivo bam, escriba este comando en el terminal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convertir el archivo SAM a BAM\n",
    "!samtools sort GPSC1bwa.sam -o GPSC1bwa.bam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La explicación de este comando es la siguiente:\n",
    "\n",
    "`samtools`: es la herramienta\n",
    "\n",
    "`sort`: indica a samtools que ordene el archivo SAM (GPSC1bwa.sam)\n",
    "\n",
    "`GPSC1bwa.sam`: archivo de entrada\n",
    "\n",
    "`-o GPSC1bwa.bam`: flag para el archivo de salida llamado GPSC1bwa.bam que son los datos ordenados en formato BAM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Paso 4: Indexar el archivo BAM\n",
    "\n",
    "la indexación, que se basa en datos ordenados, permite realizar búsquedas más rápidas. Escriba este comando en el terminal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Indexar el archivo BAM\n",
    "!samtools index GPSC1bwa.bam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La explicación de este comando es la siguiente:\n",
    "\n",
    "`samtools`: es la herramienta\n",
    "\n",
    "`index`: indica a samtools que indexe el archivo de entrada (GPSC1bwa.bam)\n",
    "\n",
    "`GPSC1bwa.bam`: archivo de entrada\n",
    "\n",
    "Ahora debería haber dos nuevos archivos llamados GPSC1bwa.bam y GPSC1bwa.bam.bai (el archivo índice BAM) en el directorio. Ahora vamos a listar `ls -alh` el contenido del directorio para comprobar que tenemos nuestros nuevos ficheros, y también comprobar sus tamaños.\n",
    "\n",
    ">Nota: Si su archivo SAM es 0B (es decir, 0 bytes) entonces algo salió mal con el paso de alineación bwa, así que reinicie desde allí. Si su archivo SAM está bien (es decir, >0), pero su archivo BAM es 0B (es decir, vacío), entonces algo salió mal con su conversión de SAM a BAM, así que vuelva a hacer esa sección.\n",
    "\n",
    "Ya no necesitamos nuestro archivo SAM original (ya que ahora tenemos el archivo BAM). Así que eliminamos `rm` el archivo SAM (GPSC1bwa.sam)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Paso 5: Evaluación de la alineación\n",
    "\n",
    "Samtools es uno de los programas clave en los análisis de datos de secuenciación de alto rendimiento (HTS). Tiene una amplia gama de funciones y se combina fácilmente con herramientas relacionadas como vcftools (para la llamada de variantes). Una cosa que se suele comprobar es cuántas lecturas se han alineado con la referencia (mapeadas) y cuántas no (no mapeadas). Samtools puede informarnos de esto fácilmente, utilizando las banderas SAM del alineador que aprendió en la sección 2 (formatos de archivos NGS).\n",
    "\n",
    "La segunda columna del archivo SAM contiene la bandera para la alineación de la lectura. Si el flag incluye el número 4 en su composición entonces la lectura no está mapeada, si no incluye el número 4 entonces está mapeada.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Número de alineaciones de lectura mapeadas**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ver las primeras líneas del archivo BAM\n",
    "!samtools view -c -F4 GPSC1bwa.bam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La explicación de este comando es la siguiente\n",
    "\n",
    "`samtools view`: para ver el archivo bwa.bam\n",
    "\n",
    "`-c`: contar la alineación de lectura\n",
    "\n",
    "`-F4`: omitir la alineación de lecturas que tienen el Flag 4 unmapped"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Número de lecturas no mapeadas**\n",
    "\n",
    "Esta vez utilizamos `-f4`, que sólo incluye la alineación de lectura que tiene la bandera no mapeada 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ver las primeras líneas del archivo BAM\n",
    "!samtools view -c -f4 GPSC1bwa.bam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Técnicamente, el comando anterior da el número de alineaciones de lectura mapeadas, no de lecturas. Una lectura puede estar mapeada igualmente bien en múltiples posiciones (una se llamará alineamiento primario, y otras alineamiento secundario [sam flag 256], o una lectura puede estar dividida en dos partes (por ejemplo, empalme) siendo una parte el alineamiento primario y las otras suplementarias [sam flag 2048].\n",
    "\n",
    "Por lo tanto, para obtener el número real de lecturas mapeadas es necesario contar los alineamientos que no tienen las banderas 4 (no mapeado), 256 (no primario) y 2048 (suplementario) = 4+256+2048 = 2308"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Número de lecturas mapeadas**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!samtools view -c -F2308 GPSC1bwa.bam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!samtools view -c -F4 -F256 -F2048 GPSC1bwa.bam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Estadísticas**\n",
    "\n",
    "Puede generar las estadísticas de mapeo utilizando el comando:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!samtools stats GPSC1bwa.bam > GPSC1bwa_bamstats.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Abra el archivo \"GPSC1bwa_bamstats.txt\" para ver las estadísticas de todas las lecturas."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parte 3: Visualizacion de mapeo en IGV\n",
    "\n",
    "Descargue [IGV](https://igv.org/doc/desktop/#DownloadPage/) \n",
    "\n",
    "Necesitará los siguientes archivos para ver las lecturas mapeadas en IGV:\n",
    "\n",
    "- Secuencia de referencia que utilizó para mapear sus lecturas\n",
    "\n",
    "- Archivo BAM de lecturas mapeadas ordenadas e indexadas\n",
    "\n",
    "Inicie IGV\n",
    "\n",
    "1. Cargue la secuencia de referencia: En la barra de herramientas, haga clic en Genoma > Cargar genoma desde archivo > Buscar y seleccione Reference_sequence_GPSC1.fa.\n",
    "\n",
    "2. Cargue el archivo BAM: Vaya a Archivo > Cargar desde archivo > Seleccione GPSC1bwa.bam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ahora es libre de investigar diferentes áreas y las alineaciones en el genoma."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parte 4: Ahora realize un nuevo alienamiento seguiendo los pasos anteriores utilizando como genoma de referencia `Reference_sequence_GPSC33`\n",
    "\n",
    "> Nota: Esencialmente, usted tendrá que indexar el nuevo genoma de referencia, cambiar los nombres por el nuevo genoma de referencia en el comando bwa y todos los nombres de archivo SAM/BAM en los comandos bwa y samtools de GPSC1 a GPSC33. Ejemplo:\n",
    "\n",
    "`bwa mem Reference_sequence_GPSC33.fa ERR2667737_1.fastq.gz ERR2667737_2.fastq.gz > GPSC33bwa.sam`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Adaptado de:*\n",
    "\n",
    "- Advanced Bioinformatics Course developed for the GPS and JUNO projects - Wellcome Sanger Insitute\n",
    "    \n",
    "*Modificado por Luisa Sacristán (Universidad de los Andes-CABANA)*"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}