Metagenoma Nanopore: Guía de comandos para análisis funcional

Tienes secuencias provenientes de una secuenciación masiva Nanopore, esta guía te servirá para realizar un análisis de funcionalidad. Esta guía también te sirve para secuencias llumina pues realizo una comparación de ambas técnicas..

No es una guía de conceptos y procedimientos
Sólo son los comandos con algunas explicaciones útiles para quienes saben del tema o tiene un jomaldon para guiarlos.

Glosario de extensiones de archivos

Archivo .gz = Cualquier archivo comprimido con el programa gzip
Archivo .tar =Paquete de archivos (Un archivo que en su interior tiene varios archivos y carpetas) pero no está comprimido.
Archivo .tar.gz o .tgz = Es un archivo .tar comprimido con gzip.
Archivos .fastq o .fq = Archivos de secuencias en un formato que almacena secuencia y calidad.
Archivos .fasta o .fa = Archivos de secuencia en un formato que no almacena calidad.
Archivos .qza: Archivos comprimidos (.zip) que contienen toda la info que sale desde un determinado proceso de Qiime2. Se puede descomprimir (winzip, winrar).
Archivos .qzv: Archivo comprimido (.zip) que contiene información para "visualizar" los resultados de un archivo .qza. Se puede abrir arrastrándolo a la interfáz web view.qiime2.org o, se puede descomprimir (winzip, winrar) y revisar su contenido. De hecho, dentro de la carpeta "data" encontrarás un archivo index.html el cual puedes abrir y se desplegará en tu navegador web toda la info y gráficos correspondientes al tipo de .qza de origen.
Archivos .tsv: Archivos de texto (.txt) donde las columnas está separadas por tabuladores. Se puede abrir en Excel como un archivo .csv.

1. Archivos de ejemplo

1.1 Archivos de secuencias en formato fasta o fasta.gz

Vamos a utilizar como ejemplo un set de secuencias de metagenomas obtenidos desde superficie marina del Océano Atlántico disponibles en NCBI las cuales fueron secuenciadas mediante Nanopore y mediante Illumina con el fin de identificar virus marinos. Se trata de tres muestras que representan el día 7, 9 y 10 de una expedición de la cual no tenemos mayores datos en NCBI.

Este es el link al bioproject para obtener mas información: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA1130046

Para facilitar el tutorial y usar menos recursos de cómputo sólo usaremos la muestra del día 7 pero todos los comandos se pueden replicar con las muestras restantes.

El set de datos es el siguiente con sus respectivos códigos en la base de datos SRA de NCBI:

Atlantic Ocean water, Día 7, Nanopore metagenome
SRR29656297
Atlantic Ocean water, Día 7, Illumina metagenome
SRR29656296

Lo primero que haremos será descargar estas secuencias en nuestro disco duro utilizando un ambiente conda.

conda create -n nanopore1 -c conda-forge -c bioconda -c cyclus parallel-fastq-dump porechop fastp kraken2 seqtk flye unicycler canu NanoPlot megahit NanoFilt java-jdk fastqc multiqc prokka pigz

# Activamos el ambiente conda
conda activate nanopore1

# Archivos para descarga "abreviada" (sólo 100 mil secuencias)
wget http://genius.bio.puc.cl/workshop/SRR29656295.100k.fastq.gz; mv SRR29656295.100k.fastq.gz SRR29656295.fastq.gz
wget http://genius.bio.puc.cl/workshop/SRR29656297.100k.fastq.gz; mv SRR29656297.100k.fastq.gz SRR29656297.fastq.gz
wget http://genius.bio.puc.cl/workshop/SRR29656299.100k.fastq.gz; mv SRR29656299.100k.fastq.gz SRR29656299.fastq.gz

# Archivos para descarga "abreviada" (sólo 1 millón de secuencias)
wget http://genius.bio.puc.cl/workshop/SRR29656295.1M.fastq.gz; mv SRR29656295.1M.fastq.gz SRR29656295.fastq.gz
wget http://genius.bio.puc.cl/workshop/SRR29656297.1M.fastq.gz; mv SRR29656297.1M.fastq.gz SRR29656297.fastq.gz
wget http://genius.bio.puc.cl/workshop/SRR29656299.1M.fastq.gz; mv SRR29656299.1M.fastq.gz SRR29656299.fastq.gz

# Si Ud. descargó los archivos de descarga "abreviada" siga al paso 1.2

# Descargamos y descomprimimos la herramienta de NCBI que permite manejar archivos de secuencias
wget https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/current/sratoolkit.current-ubuntu64.tar.gz
tar xfz sratoolkit.current-ubuntu64.tar.gz

# Ejecutamos el comando para descargar los archivos
# El parámetro -pe indica el número de hilos de procesamiento a utilizar. Deben cambiar ese parámtro según corresponda.

sratoolkit.3.1.1-ubuntu64/bin/prefetch --max-size 30GB SRR29656297 && sratoolkit.3.1.1-ubuntu64/bin/fasterq-dump -pe 4 ./SRR29656297
sratoolkit.3.1.1-ubuntu64/bin/prefetch --max-size 30GB SRR29656296 && sratoolkit.3.1.1-ubuntu64/bin/fasterq-dump -pe 4 ./SRR29656296

La secuencia Illumina (SRR29656296) descargará dos archivos fastq, SRR29656296_1.fastq que corresponde a la secuencias forward y SRR29656296_2.fastq que corresponde a la secuencia reverse.

# comprimimos todos los archivos para usar menos recursos de disco
pigz *.fastq

1.2 Filtrado/limpieza de secuencias

Ahora debemos remover las secuencias adaptadoras y filtrar calidad de las secuencias lo cual realizaremos en el ambiente nanopore1.

# Activamos el ambiente conda
conda activate nanopore1

## Nanopore
# Utilizamos el programa porechop para remover adaptadores de las secuencias Nanopore
# El parámetro -t indica el número de hilos de procesamiento a utilizar. Deben cambiar ese parámtro según corresponda.

porechop -i SRR29656297.fastq.gz --check_reads 1000 -t 4 -v 1 | gzip > SRR29656297.wo_adapt.fastq.gz

gzip SRR29656297.wo_adapt.fastq.gz | NanoFilt -q 10 -l 500 --logfile SRR29656297.nanofilt.log | gzip > SRR29656297.wo_adapt.q10.l500.fastq.gz

## Illumina
# Utilizamos el programa Trimmomatic para remover adaptadores de las secuencias Illumina
# http://www.usadellab.org/cms/index.php?page=trimmomatic
# El parámetro -threads indica el número de hilos de procesamiento a utilizar. Deben cambiar ese parámtro según corresponda.

# Descargar y descomprimir el programa Trimmomatic
wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.39.zip
unzip Trimmomatic-0.39.zip

# Copiar el archivo TruSeq3-PE.fa en el directorio de trabajo
cp Trimmomatic-0.39/adapters/TruSeq3-PE.fa .

# Ejecutar trimmomatic
java -jar Trimmomatic-0.39/trimmomatic-0.39.jar PE -threads 4 SRR29656296_1.fastq SRR29656296_2.fastq -baseout SRR29656296_trimmed.fastq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:20 TRAILING:20 SLIDINGWINDOW:10:30 MINLEN:50 AVGQUAL:25;

1.4 Revisión de calidad de secuencias

# El parámetro -t indica el número de hilos de procesamiento a utilizar. Deben cambiar ese parámtro según corresponda.

## Nanopore
# Usaremos el programa NanoPlot para obtener gráficos que representan la calidad de la secuencia.
NanoPlot --fastq SRR29656297.fastq.gz -t 4 --plots kde --maxlength 10000 -o Plot10000.raw.SRR29656297
NanoPlot --fastq SRR29656297.wo_adapt.q10.l500.fastq.gz -t 4 --plots kde --maxlength 10000 -o Plot1000.trim.SRR29656297

## Illumina
# Usaremos el programa fastqc
fastqc -t 4 SRR29656296_1.fastq.gz SRR29656296_2.fastq.gz -o fastqc_results
fastqc -t 4 SRR29656296_trimmed_1P.fastq.gz SRR29656296_trimmed_2P.fastq.gz -o fastqc_results

# Juntamos los resultados fastqc en un único informe usando el comando multiqc
multiqc fastqc_results/*.zip -o fastqc_results

# El reporte conjunto estará en el earchivo fastqc_results/multiqc_report.html

2. Análisis de secuencias

Ahora que tenemos secuencias "limpias" podemos proceder de dos formas:

Mapear los reads contra una base de datos para obtener información de las funciones que están presentes en el DNA
Ensamblar los reads para obtener secuencias de genomas (contigs y scaffolds) y luego hacer predicción génica. Con esta predicción génica podemos realizar la anotación funcional para conocer las funciones moleculares presentes en el DNA.

2.1 Método 1: Mapeo de reads contra referencia.

Si queremos obtener el máximo de información desde nuestras secuencias entonces lo recomendable sería usar toda la base de datos no redundante de NCBI (NR) pero esto es muy costoso computacionalmente y requiere mucho tiempo. Por lo anterior, vamos a elegir un subconjunto de la base de datos NR que sólo contenga las secuencias de virus.

2.1.1 Descarga de base de datos de virus ya preparada

He preparado una base de datos centrada en secuencias de virus, llegar y usar en el paso 2.1.2. Si deseas saber cómo se prepara la base de datos puedes saltar al punto 2.1.3

# Base de datos de Virus (1,9 Gb comprimido, 5 Gb descomprimido)
wget http://genius.bio.puc.cl/workshop/nr.virus.dmnd.gz 
pigz -d nr.virus.dmnd.gz

# Base de datos de Virus, Bacteria, Arquea y Hongos (80 Gb comprimido, 150 Gb descomprimido)
# wget http://genius.bio.puc.cl/workshop/nr.virus-bacteria-archaea-fungi.dmnd.gz

2.1.2 Mapeo de reads contra la base de datos.

Ahora vamos a buscar cada secuencia de nuestros archivos de metagenoma dentro de la base de datos de virus.
A este proceso lo llamamos "mapeo", el cual básicamente es una búsqueda por homología ultra rápida.

# Descarga e instalación del programa diamond
wget https://github.com/bbuchfink/diamond/releases/download/v2.1.9/diamond-linux64.tar.gz
tar xfz diamond-linux64.tar.gz

# Mapeo de secuencias

## Nanopore
./diamond blastx -d nr.virus.dmnd -q SRR29656297.wo_adapt.q10.l500.fastq.gz -o SRR29656297.wo_adapt.q10.l500.daa -f 100 -p 4 -b 12 -c 1 --parallel-tmpdir /projects4/tmp/ -t /projects4/tmp/ --long-reads

## Illumina
./diamond blastx -d nr.virus.dmnd -q SRR29656296_trimmed_1P.fastq.gz -o SRR29656296_trimmed_1P.daa -f 100 -p 4 -b 12 -c 1 --parallel-tmpdir /projects4/tmp/ -t /projects4/tmp/

./diamond blastx -d nr.virus.dmnd -q SRR29656296_trimmed_2P.fastq.gz -o SRR29656296_trimmed_2P.daa -f 100 -p 4 -b 12 -c 1 --parallel-tmpdir /projects4/tmp/ -t /projects4/tmp/

# Descargamos el programa MEGAN para linux y lo instalamos en la carpeta "megan"
wget https://software-ab.cs.uni-tuebingen.de/download/megan6/MEGAN_Community_unix_6_25_10.sh
bash MEGAN_Community_unix_6_25_10.sh

# Descargamos la base de datos de MEGAN (1.8 Gb que se descomprimen en 8.7 Gb)
wget https://software-ab.cs.uni-tuebingen.de/download/megan6/megan-map-Feb2022.db.zip
unzip megan-map-Feb2022.db.zip

# Transformamos el archivo de mapeo en un formato que podremos abrir con el programa Megan
~/megan/tools/daa-meganizer -i *.daa -mdb megan-map-Feb2022.db -t 4 -tsm

# Ahora será posible abrir y analizar los archivos en el software MEGAN (que tiene versión para Windows).

2.1.3 Preparación de base de datos de virus (NO HACER)

# Preparación de base de datos de virus

# Instalamos un ambiente conda que nos ayudará en el proceso
conda create -n NRtoTax -c conda-forge -c bioconda taxonkit csvtk seqkit aria2c
conda activate NRtoTax

# Vamos a obtener el ID de todos los organismos catalogados como "virus" en la base de datos taxonomy de NCBI
# El ID de virus es 10239. Si queremos obtener el ID de bacterias debemos reemplazar el ID 10239 por el ID 2.
taxonkit list -i 10239 -I "" -o virus.taxid.txt

# El archivo virus.taxid.txt contiene todos los IDs de virus que están contenidos en NCBI.

# Descargamos la base de datos NCBI (mas de 134 Gb!)
aria2c -s10 -x10 ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz
wget ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz.md5

# Usando los IDs de virus vamos a buscar todos los IDs de secuencias que pertenecen a esas taxas
zcat prot.accession2taxid.gz | 
                      csvtk -t grep -f taxid -P virus.taxid.txt | 
                      csvtk -t cut -f accession.version > virus.taxid.acc.txt

# El archivo virus.taxid.acc.txt tendrá todos los códigos de genes de NCBI que pertenecen a virus.

# Ahora vamos a rescatar la secuencia de esos genes desde el archivo de secuencias de la base de datos nr.
gzip -c -d nr.gz | seqkit grep -f virus.taxid.acc.txt | sed -e 's/^>gb|/>/' -e  's/|//' > nr.virus.fa

# El archivo nr.virus.fa contiene todas las secuencias de genes de NCBI que pertenecen a virus.

# Finalmente vamos a crear la base de datos usando el programa diamond (puede demorar un par de horas).
wget https://github.com/bbuchfink/diamond/releases/download/v2.1.9/diamond-linux64.tar.gz
tar xfz diamond-linux64.tar.gz
./diamond makedb --threads 4 --in nr.virus.fa --db nr.virus --taxonnodes /home/jomaldon/.taxonkit/nodes.dmp

# El archivo nr.virus.dmnd es la base de datos indexada de los genes de virus obtenidos desde de NCBI

2.2 Método 2: Ensamblar genomas.

Existen distintas estrategias para realizar esta tarea. En estas líneas usaré el programa Flye.

## Nanopore
flye --nano-raw SRR29656297.wo_adapt.q10.l500.fastq.gz --out-dir SRR29656297.flye --meta --threads 4 1>SRR29656297.flye.log1 2>SRR29656297.flye.log2

# Vamos a renombrar el archivo de ensamble para tener claro su origen
cd SRR29656297.flye
mv contigs.fa SRR29656297.flyeMeta.fasta

## Illumina
megahit -1 SSRR29656296_trimmed_1P.fastq.gz -2 SRR29656296_trimmed_2P.fastq.gz -o SRR29656297.megahit -t 4 1>SRR29656297.megahit.log1 2>SRR29656297.megahit.log2

cd SRR29656297.megahit
mv final.contigs.fa SSRR29656296.MegaHit.fasta

2.3 Anotación estructural (predicción de genes) y anotación funcional.

Vamos a usar el programa Prokka para "predecir" genes desde nuestros metagenomas y para anotarlos funcionalmente.

## Nanopore
cd SRR29656297.flye
prokka SRR29656297.flyeMeta.fasta --outdir PROKKA_SRR29656297 --prefix SRR29656297 --cpus 4--metagenome --norrna --notrna

## Illumina
cd SRR29656297.megahit
prokka SSRR29656296.MegaHit.fasta --outdir PROKKA_SSRR29656296 --prefix SSRR29656296 --cpus 4--metagenome --norrna --notrna

El resultado de la anotación estructural y funcional son archivos con diferentes formatos lo cual es se pueden revisar por separado o en programas como IGV o UGENE.