En esta actividad veremos como ensamblar el genoma del virus causante del COVID19, el SARS-CoV-2; este virus tiene un tamaño de 29,903 pb. La secuenciación de este virus es a partir de cDNA el cual se amplifica con múltiples juegos de primers en dos pooles de PCR, los amplicones obtenidos por estos dos pooles se unen y se secuencian ya sea en plataforma Illumina o Nanopore; aquí lo haremos con secuencias Illumina obtenidas en nuestro laboratorio.

Este pipeline se basa en https://github.com/kkosaki/covid19_genome_analysis

Software

Para este ejercicio necesitaremos los siguientes programas:

• fastp A tool designed to provide fast all-in-one preprocessing for FastQ files.

• bwa BWA is a software package for mapping DNA sequences against a large reference genome.

• samtools Tools for manipulating next-generation sequencing data.

• bedtools A powerful toolset for genome arithmetic.

• ivar A computational package that contains functions broadly useful for viral amplicon-based sequencing.

Estos programas ya se encuentran instalados en la imágen virtual AnvioVM, todos los set de datos se pueden bajar directamente desde la terminal en la imágen virtual: $ download_data descargará los datos a la carpeta Documents/SARSCoV2/

Archivos necesarios

Si no se bajaron los datos como se mencionó anteriormente, entonces bajar el set de datos de Figshare a su archivo a su directorio de trabajo y descomprimirlo

$ tar xzf SARSCoV2.tar.gz

• En este set de datos tenemos tres set de datos de variantes de SARS-CoV-2:

  • Alpha

  • Delta

  • Omicron)

• Además necesitaremos un genoma de referencia (cepa Wuhan, NC_045512.fasta) y en formato gff (NC_045512.gff) y secuencias de los primers usados (Artic.v3.fasta); estos archivos ya están en la imágen virtual en /opt/COVID

• Para este ejercicio usaremos las secuencias pair-end Illumina comprimidas (.gz) de la muestra clasificada como variante Alpha con una cobertura de 362X que se encuentran en el subdirectorio Alpha:

  • Alpha_S001_L001_R1_001.fastq.gz

  • Alpha_S001_L001_R2_001.fastq.gz

Análisis bioinformático

Si estamos usando la imágen virtual AnvioVM, los programas necesarios están instalados en CONDA, por lo que habrá que activarla:

$ conda activate anvio6

Limpieza de secuencias

Lo primero que tenemos que hacer es limpiar las secuencias de nuestra muestra (aquí denominada Alpha), esto lo podemos hacer con fastp y usaremos ambas secuencias pair-end de Illumina:

$ fastp -i Alpha_S001_L001_R1_001.fastq.gz -o alpha.R1.fq -I Alpha_S001_L001_R2_001.fastq.gz -O alpha.R2.fq -q 30 -l 20 --low_complexity_filter --adapter_fasta /opt/COVID/Artic.v3.fasta --thread 2 -h alpha-fastp.html -R "alpha fastp trimming report"

Aquí estamos usando los dos archivos fastq de entrada (con -i y -I) pair-end y como salida (-o y -O) archivos tipo fastq pero con terminación .fq sólo para distinguirlos de los de entrada.

Limpiaremos secuencias con calidad menor a Q30 (-q 30), eliminaremos secuencias cortas (-l 20) y filtrar para baja complejidad (--low_complexity_filter).

Eliminaremos los primers (--adapter_fasta /opt/COVID/Artic.v3.fasta) ya que éstos pueden tener bases mal llamadas que no son reales.

Usaremos dos núcleos (--thread 2) y crearemos un archivo de reporte tipo html (-h alpha-fastp.html) y ponerle un título a ese reporte (-R "Alpha fastp trimming report").

Mapeo de secuencias a una referencia

Ahora debemos mapear esas secuencias limpias a un genoma de SARS-CoV-2 de referencia (NC_045512.fasta); claro que tenemos que tener la secuencia de la referencia en el subdirectorio en el que estemos, para ello primero la copiaremos del directorio al subdirectorio Alpha:

$ cp ../NC_045512.fasta .

$ bwa index NC_045512.fasta

$ bwa mem -t 4 NC_045512.fasta alpha.R1.fq alpha.R2.fq | samtools sort -O BAM - > assembly.sorted.bam

$ samtools index assembly.sorted.bam

$ samtools faidx NC_045512.fasta

Secuencia consenso

Ya te sabemos en donde van cada una de las secuencias usadas, podemos crear ahora una secuencia consenso (ensamble) de nuestra muestra (alpha) que en este ejemplo, la secuencia consenso en formato fasta se llamará alpha-consensus.fa

$ samtools mpileup -d 50000 --reference NC_045512.fasta -a -Q 30 assembly.sorted.bam | ivar consensus -t 0 -m 2 -n N -p alpha-consensus.fa

Como vemos en el resultado de Nextclade (Fig arriba), la secuencia consenso fue clasificada como Alpha, con 36 mutaciones, 136 nucleótidos no resueltos (Ns), 9 gaps y un codón Stop prematuro (SC).

Cobertura

Podemos también obtener información para saber en que lugares del genoma quedó una baja cobertura o profundidad con bedtools:

$ bedtools genomecov -bga -ibam assembly.sorted.bam | grep -w '0$\|1$\|2$\|3$' > alpha-lowcoverage.bed

Ahora vamos a obtener la cobertura o profundidad de las secuencias a lo largo del genoma de referencia.

$ bedtools genomecov -d -split -ibam assembly.sorted.bam | cut -f2,3 | sed 's/\t/,/' | sed '1 i\position,coverage' > alpha-coverage.csv

Este archivo delimitado por comas tiene el valor de profundidad para cada nucleótido y puede ser graficado fácilmente con R (ver abajo).

Variantes

Por último podemos obtener las variaciones que tiene nuestra muestra con respecto al genoma de referencia (cepa original de Wuhan)

$ samtools mpileup -A -d 0 --reference NC_045512.fasta -B -Q 0 assembly.sorted.bam | ivar variants -p variants -q 20 -t 0.25 -r NC_045512.fasta -g NC_045512.gff

Esto nos crea un archivo (variants.tsv) con la información de cambios que encontró en nuestra muestra con respecto a la referencia.

Por último podemos limpiar archivos intermedios que ya no nos son útiles:

$ rm *.fasta.* fastp.json assembly.sorted.bam* *.txt

Gráficas

Para graficar la profundidad con R, podemos en RStudio

importar el archivo en RStudio, seleccionar del menú:

Import dataset > From Text (base) ... y seleccionar alpha-coverage.csv

cargar ggplot2:

> library(ggplot2)

Crear una variable:

> df <- alpha-coverage

Construir la gráfica con R:

> ggplot(data = df, aes(position, coverage)) +

ggtitle("alpha") +

scale_y_continuous(trans = 'log10', name = "Profundidad (X)") + scale_x_continuous(name = "Genoma SARS-CoV-2") +

geom_line(colour="grey50") + geom_area(fill="grey", alpha=0.6) + geom_hline(yintercept = 10, colour = "grey60", size = 0.5)

Se puede tardar unos segundos en generar la gráfica