PICRUSt

PICRUSt2 (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) nos permite hacer una predicción funcional a partir de solo de datos de 16S. Ojo, no es equivalente a una análisis de metagenómica funcional pero nos puede dar una buena idea de que funciones están presentes en nuestras muestras.

Procedimiento

Es necesario crear dos archivos con un formato específico para ser usado por PICRUSt2:

1. Secuencias centroides en formato fasta de nuestras muestras.

2. Numero de secuencias homólogas a estas secuencias centroides en cada una de las muestras a analizar.

Limpieza de secuencias

Si no tenemos ya las secuencias limpias y en formato fasta, debemos primero hacerlo, ver el script pair-end_cleaner. Este script nos creará una nueva carpeta con las secuencias en formato fasta llamada assembled... (añade la fecha y hora al nombre de la carpeta).

Entremos en la carpeta
$ cd assembled...

Creación de archivos

Teniendo en una carpeta nuestros archivos fasta, podemos ejecutar el script centroids_assignments entrando en la carpeta:

$ centroids_assignments *.fasta

éste nos generará dos archivos: centroid_sequences.fna y num_sequences.tsv.

PICRUSt2

Teniendo nuestros dos archivos de entrada, podemos clasificar ya con PICRUSt2 activando antes el ambiente conda picrust2:

$ conda activate picrust2

$ picrust2_pipeline.py -s centroid_sequences.fna -i num_sequences.tsv -o picrust2 -p 8 --verbose

Este análisis dura bastante minutos por lo que es recomendable correrlo con screen.

Una vez terminado tendremos una nueva carpeta llamada picrust2 con los resultados, podemos éstos mejorarlos añadiéndole descripciones a las funciones encontradas ejecutando los siguientes comandos desde dentro de esa carpeta:

$ cd picrust2

$ add_descriptions.py -i EC_metagenome_out/pred_metagenome_unstrat.tsv.gz -m EC -o EC_metagenome_out/pred_metagenome_unstrat_descrip.tsv.gz

$ add_descriptions.py -i KO_metagenome_out/pred_metagenome_unstrat.tsv.gz -m KO -o KO_metagenome_out/pred_metagenome_unstrat_descrip.tsv.gz

$ add_descriptions.py -i pathways_out/path_abun_unstrat.tsv.gz -m METACYC -o pathways_out/path_abun_unstrat_descrip.tsv.gz

Algunos de los comandos que realiza el script centroids_assignments son los siguientes:

Secuencias centroides

Son secuencias únicas presentes en todas nuestras muestras, aquellas con una similitud igual o superior a 97%. Para obtenerlas podemos usar vsearch, primero concatenemos todos los archivos fasta de nuestras muestras:

$ cat *.fasta > all_sequences.fa

ahora obtengamos las secuencias centroides con una similitud del 97% (0.97):

$ vsearch --cluster_size all_sequences.fa --id 0.97 --strand plus --centroids centroides.temp --uc clusters.uc --threads 4

renombremos las secuencias centroides:

$ awk '/^>/ {print ">centroid_" ++i; next} {print}' centroides.temp > centroides.fa

Asignación de secuencias

Teniendo ya un archivo que clasifica todas las secuencias de todas las muestras a cada secuencia centroide (clusters.uc) generado en el paso de vsearch (opción --uc clusters.uc), podemos extraer de éste la información para calcular cuantas secuencias de cada muestra son asignadas a cada secuencia centroide. 

en construcción