En esta actividad vamos a ensamblar un genoma bacteriano contra un genoma de referencia.
Vamos a descargar las secuencias fastq producidas con un equipo Illumina desde el repositorio SRA del NCBI. Para hacer esto tenemos que instalar primero el paquete sratoolkits que produce el NCBI.
Desde la línea de comando ir a HOME y descargar los archivos para instalar sratoolkits
cd ~/
Ubuntu 64-bit (no hay versión para 32 bits)
wget --output-document sratoolkit.tar.gz http://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/current/sratoolkit.current-ubuntu64.tar.gz
Para OS X
wget --output sratoolkit.tar.gz http://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/current/sratoolkit.current-mac64.tar.gz
Descomprimir el archivo descargado:
tar -vxzf sratoolkit.tar.gz
Cambiar el nombre del directorio recién creado a sratoolkit
mv sratoolkit.2.6.3-ubuntu64 sratoolkit
Ejecutar el siguiente comando y además copiarlo al final del archivo .bashrc
export PATH=$PATH:$HOME/sratoolkit/bin
Para verificar que está todo en orden correr este comando
fastq-dump
Existen diferentes programas para hacer ensamblados contra una referencia. Bowtie es un programa que existe desde hace tiempo y es bastante popular. Para genomas grandes hoy existen otras alternativas, pero necesitan máquinas más potentes. Para el ejercicio que vamos hacer nosotros, que es ensamblar un genoma bacterianos, va a funcionar bien.
En Ubuntu
sudo apt-get install bowtie2
Para OS X:
Intentar primer con Homebrew. Si Bowtie2 no está disponible, ir a la carpeta de descarga de Bowtie2 (https://sourceforge.net/projects/bowtie-bio/files/bowtie2/) y descargar el archivo binario correspondiente a la última versión
Después de instalarlo, agregar la dirección de la carpeta con los archivos ejecutables (bowtie2, bowtie2-build y bowtie2-inspect) a $PATH, y luego crear la variable de entorno BT2_HOME que apunta al directorio que contiene los mismos archivos anteriores.
En Ubuntu
sudo apt-get install samtools
Para OS X:
Intentar primer con Homebrew. Si no funciona, descargar desde http://www.htslib.org/ e instalar.
sudo apt-get install igv
Creamos un directorio llamado taller_4 y hacerlo nuestro directorio de trabajo
mkdir ~/taller_4
cd ~/taller_4
A continuación descargarmos desde el repositorio SRA del NCBI las secuencias fastq obtenidas a partir de la secuenciación de un aislamiento de Mycobacterium bovis con tecnología Illumina y estrategia paired-end (http://www.ncbi.nlm.nih.gov/sra/?term=SRR3135071).
fastq-dump --split-files SRR3135071
Mientras esto se descarga, vamos a ir con el navegador al link anterior para ver que información proporciona el NCBI para estos registros y entender como se almacena la información de este tipo de proyectos.
Con fastq-dump descargamos dos archivos fastq, cada uno con la secuencias de uno de los extremos de los fragmentos paired end.
El próximo paso es hacer un alineamiento utilizando un genoma de referencia con Bowtie2. Para eso necesitamos la secuencia fasta del genoma de referencia y debido a que no podemos acceder directamente por ftp a los repositorios del NCBI, vamos a ir a la página web de la secuencia nucleotídica de Mycobacterium bovis AF2122/97:
https://www.ncbi.nlm.nih.gov/nuccore/NC_002945.3
Y desde ahí vamos a descargar el archivo de la secuencia en formato fasta. La secuencia que descargamos va a tener como nombre nombre sequence.fasta, lo cambiamos a NC_002945.fasta.
Dentro de taller_3 creamos un directorio ref y nos movemos ahi
mkdir ref
cd ref
Movemos a esta carpeta el archivo NC_002945.fasta, y corremos el comando que sigue para construir el índice del genoma de referencia:
bowtie2-build NC_002945.fasta M_bovis_ref
Volvemos al directorio superior
cd ..
Y ya podemos construir el ensamblado contra la referencia:
bowtie2 -p 2 -x ref/M_bovis_ref -1 SRR3135071_1.fastq -2 SRR3135071_2.fastq -S ensamble.sam
El archivo ensamble.sam es un archivo de texto con toda la información sobre el mapeo de nuestros reads paired-end sobre el genoma de referencia. Hacemos un backup de este archivo:
cp ensamble.sam ensamble.bak
Los archivos SAM suelen ser muy grandes. Por ejemplo, el que construimos tiene 473 Mb, es común transformarlos a su versión binaria, que ocupan menos espacio. Tienen una extensión .bam, y se crean usando una de las herramientas de samtools:
samtools view -bST ref/NC_002945.fasta ensamble.sam > ensamble.bam
Por ejemplo, para ver el encabezado del alineamiento:
samtools view -H ensamble.bam
Las operaciones que se pueden hacer con un archivo SAM, también se pueden hacer con uno BAM. Pero hay algunas operaciones que sólo se pueden hacer con archivos BAM ordenados.
samtools sort ensamble.bam ensamble_sorted
samtools view -H ensamble_sorted.bam
También podríamos eliminar todos los reads no mapeados del ensamble ordenado:
samtools view -h -F 4 -b ensamble_sorted.bam > ensamble_map.bam
Para visualizar archivos BAM, por ejemplo con igv, normalmente hay que construir primero un índice:
samtools index ensamble_map.bam ensamble_map.bai
Podemos generar una secuencia consenso y descargarla como archivo fasta. Para esto primero hay que crear una variable de entorno y aagregarla a .bashrc:
export PATH=$PATH:/usr/share/samtools/
Luego corremos:
samtools mpileup -uf ref/NC_002945.fasta ensamble_map.bam | bcftools view -cg - | vcfutils.pl vcf2fq > consenso.fastq
Para convertir este archivo a formato fasta, instalamos seqtk
sudo apt-get install seqtk
seqtk seq -A consenso.fastq > consenso.fasta
Esta secuencia consenso fasta se puede anotar con herramientas on-line como RAST (http://rast.nmpdr.org/) o la herramienta de anotación del NCBI (http://www.ncbi.nlm.nih.gov/genome/annotation_prok/) o con programas que se pueden instalar localmente como Prokka (http://www.vicbioinformatics.com/software.prokka.shtml).
Especificación de archivos SAM/BAM (https://samtools.github.io/hts-specs/SAMv1.pdf)
Una explicación breve sobre le formato SAM (http://genome.sph.umich.edu/wiki/SAM)
Operaciones con samtools (http://davetang.org/wiki/tiki-index.php?page=SAMTools)
Los archivos SAM tienen un encabezado que tiene varias líneas que comienzan con @. Cada una de estas líneas contiene información específica del mapeo que realizamos, y en conjunto se llaman el encabezado del alineamiento:
encabezado @HD: primera línea del alineamiento
encabezado @SQ: contiene información sobre las secuencias de referencia. La clave SN contiene la identificación de la referencia, y la clave LN el largo de esa secuencia de referencia.
encabezado @PG: datos sobre el programa usado para hacer el ensamblado
encabezado @RG: información sobre grupos de secuencias
encabezado @CO: comentarios
A continuación está la sección con información de los alineamientos:
samtools view ensamble_sorted.bam | tr ' \t' '\n' | head -n 11
SRR3135071.155741
163
gi|31791177|ref|NC_002945.3|
1
23
2M13I236M
=
189
426
GGGAGATACGTCGT ....
11A>AAFBFAFFGGG....
1. QNAME: El nombre de la consulta, esto es el identificador del read.
2. FLAG: Un byte representado como un numero decimal con información del alineamiento
3. RNAME: nombre de la referencia
4. POS: posición del read sobre la referencia
5. MAPQ: calidad del mapeo
6. CIGAR: el string CIGAR describe características del alineamiento (inserciones, deleciones. coincidencias, etc.)
8. RNEXT y PNEXT: nombre y posición del read compañero del par paired-end.
9. TLEN: largo de la cadena molde para las lecturas paired-end.
10. SEQ: secuencia del read
11. QUAL: calidad de la secuencia del read.
El byte FLAG
En el listado anterior vimos que el punto 2, almacena información sobre el alineamiento. Cada número representa una o más de las condiciones que se listan en la tabla de más abajo.
Para el read de ejemplo que analizamos más arriba vemos que el FLAG es 163. Para entender como funciona esto tenemos que convertir este número a binario (se puede hacer fácil en Google entrando ""163 to binary", esto nos da:
0b10100011.
Esta secuencia la tenemos que leer de atrás para adelante y vemos que los bits "encendidos" son el 1,2 , 6 y 8. Es decir, este es un read que forma un par, que está mapeado en el orden correcto, que es el read de la cadena reversa y que es el segundo read del par.
Una herramienta aún más sencilla para entender esto: http://broadinstitute.github.io/picard/explain-flags.html
Es raro que necesitemos ver los bytes FLAG uno por uno. Una situación más usual es calcular un resumen de la información contenida en todos los reads:
samtools flagstat ensamble.bam
samtools flagstat ensamble_map.bam
El string CIGAR
El punto 6, CIGAR, describe características más detalladas del alineamiento del read con respecto a la referencia. Los códigos para describir coincidencias, inserciones, deleciones, cortes, etc se muestran a continuación:
El CIGAR de nuestro ejemplo es 2M13I236M, esto quiere decir que hay dos nucleótidos alineados (pueden matchear los dos o no), una inserción de 13 bases con respecto a la referencia y 236 bases alineadas.