Predicción de perfiles: Tax4Fun
Los comandos usados y la ciencia detrás de ellos se pueden encontrar en el artículo original.
Preparación de los datos
Para poder predecir la función de la comunidad se requieren las funciones runRefBlast y makeFunctionalPrediction del paquete Tax4Fun2. Para usar estas funciones se necesitan dos tablas que se pueden obtener de los resultados del pipeline de DADA:
- tabla con las secuencias de los ASV encontrados en las muestras (nombres de las columnas del
otu_tabledel objeto phyloseq generado) - tabla con las abundancias de dichas secuencias (
otu_tabledel objeto phyloseq generado)
Lo primero que se tiene que hacer es obtener los datos que fueron guardados como un archivo RDS (la ubicación del archivo depende de donde se haya guardado previamente):
library(Tax4Fun2)
library(phyloseq)
set.seed(1)
ps <- readRDS("~/Documents/Curso_R/data/Ejercicio_DADA")
ps## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 112 taxa and 12 samples ]
## sample_data() Sample Data: [ 12 samples by 2 sample variables ]
## tax_table() Taxonomy Table: [ 112 taxa by 7 taxonomic ranks ]En caso de que sea necesario eliminar una muestra en particular (un control, por ejemplo) es mejor hacerlo desde el principio a través del comando prune_samples y posteriormente se eliminan los ASV que solo estaban presentes en las muestras eliminadas usando el comando prune_taxa:
En este ejercicio no es necesario, por lo que se extraerá la información de abundancias de todas las muestras y se convertirá en una tabla que llamaremos foo
como se mencionó antes, en el nombre de las columnas de la otu_table se encuentran las secuencias de los ASV por lo que se extraen directamene con el comando colnames, posteriormente reemplezaremos los nombres de las columnas de la tabla de abundancias (foo) con el objetivo de que sean más fáciles de manejar:
Las tablas de secuencias y abundancias son leidas por Tax4Fun a través de un archivo externo, por lo que es necesario guardarlas en archivos independientes. Para guardar las secuencias en formato fasta empleamos la función write.fasta del paquete seqinr
seqinr::write.fasta(sequences = as.list(seqs),
names = nams,
open = 'w',
file.out = '~/Documents/Curso_R/data/Ejercicio_seqs.fasta')La tabla de abundancias es guardada con la función write.table usando el argumento sep = '\t' para indicar que esta delimitado por tabulaciones
Una vez que se tienen los archivos necesarios, ya se pueden usar los comandos para la predicción del perfil funcional, por último guardaremos en la memoria la dirección en la que se encuentra la base de datos de referencia (que se creó al instalar el paquete Tax4Fun -Tax4Fun2_ReferenceData_v2-) para evitar escribirla constantemente.
PATH = "~/Documents/Tax4Fun2_ReferenceData_v2/"Creación de la tabla de referencia
El primer paso para predecir el perfil funcional de una comunidad es alinear las secuencias de nuestras muestras con la base de datos de referencia del paquete a través de BLAST. Para ello usamos la función runRefBlast, en este caso usaremos opción database_mode = "Ref100NR" para una mayor precisión en el alineamiento, es posible usar Ref99NR como base de datos altenartiva para hacer un análisis más laxo de los datos:
runRefBlast(path_to_otus = "~/Documents/Curso_R/data/Ejercicio_seqs.fasta",
path_to_reference_data = PATH,
path_to_temp_folder = "~/Documents/Curso_R/data/Ejercicio_Ref100NR",
database_mode = "Ref100NR",
use_force = T, num_threads = 4)La función además de crear la base de datos de referencia indica cuántas de los ASV y secuencias fueron usados para crear dicha base de referencia. Esta información se puede encontrar en el archivo llamado logfileN.txt, el cual estará ubicado en la dirección que asignamos como path_to_temp_folder, en este ejemplo la ubicación del archivo es “~/Documents/Curso_R/data/Ejercicio_Ref100NR/logfile2.txt”.
Predicción del perfil funcional
Una vez que se tiene la base de datos de referencia se emplea la función makeFunctionalPrediction para obtener la predicción de la abundancia de las funciones en la comunidad
Para que la función funcione es necesario abrir el archivo de las abundancias de otus previamente guardado (Ejercicio_otu_table.txt), agregar una tabulación al principio de la primera fila y guardar el archivo
makeFunctionalPrediction(path_to_otu_table = "~/Documents/Curso_R/data/Ejercicio_otu_table.txt",
path_to_reference_data = PATH,
path_to_temp_folder = "~/Documents/Curso_R/data/Ejercicio_Ref100NR",
database_mode = "Ref100NR",
normalize_by_copy_number = T,
min_identity_to_reference = 0.97,
normalize_pathways = F)Como resultado de este comando, se generan dos archivos, uno en el cual se indican las predicciónes de las abundancias relativas de cada enzima para cada muestra de acuerdo a su composición (functional_prediction.txt) y otro que indica las abundancias de las vias metabólicas de cada muestra acorde a su composición (pathway_prediction.txt).
Ambos archivos están en formato de tabla que es facilmente leído y manejado por R para su posterior análisis. En este ejercicio analizaremos los resultados de la predicción de las vías metabólicas (pathway_prediction.txt).
El archivo pathway_prediction.txt está estructurado de la siguiente manera:
| pathway | X5622_P10_R1_MS515F | X5622_P10_R2_MS515F | level1 | level2 | level3 |
|---|---|---|---|---|---|
| ko00010 | 0.012261197 | 0.011472961 | Glycolysis/Gluconeogenesis | metabolism | Metabolism |
Como se puede observar en la tabla, las vías metabólicas están organizadas en tres niveles de jerarquías cada uno con diferentes categorías (level1-3); para revisar las categorías de los niveles de organización de las rutas metabólicas, se puede visitar la página de KEGG correspondiente.
Algunas de las categorías del primernivel son:
## [1] "level1"## [1] "Mineral absorption"
## [2] "HIF-1 signaling pathway"
## [3] "Tropane, piperidine and pyridine alkaloid biosynthesis"
## [4] "p53 signaling pathway"
## [5] "D-Glutamine and D-glutamate metabolism"Del segundo nivel son:
## [1] "level2"## [1] "Folding, sorting and degradation" "Translation"
## [3] "Cellular community - prokaryotes" "Cardiovascular diseases"
## [5] "Cell motility"En cuanto a las categorías del tercer nivel, son únicamente las siguientes:
## [1] "Cellular Processes"
## [2] "Environmental Information Processing"
## [3] "Genetic Information Processing"
## [4] "Human Diseases"
## [5] "Metabolism"
## [6] "Organismal Systems"Se puede observar que algunas de las categorías están directamente ligadas a los humanos, como es de esperarse por la naturaleza de las bases de datos empleadas (NCBI, KEGG), pero es importante mencionar que no es una buena idea eliminar de inmediato todas las rutas metabólicas que están relacionadas con enfermedades humanas, ya que en ocasiones existen rutas de señalización inter especies, de quimiotaxis y producción de metabolitos secundarios que están agrupados bajo la categoría ‘Human diseases’ que son de relevancia ecológica en el contexto de ingenieria ambiental.
foo <- read.delim(
"~/Documents/Curso_R/data/Ejercicio_Ref100NR/pathway_prediction.txt"
)
sample(levels(foo[foo$level3 == 'Human Diseases', 'level1']),25)## [1] "Colorectal cancer"
## [2] "PPAR signaling pathway"
## [3] "beta-Alanine metabolism"
## [4] "Fluorobenzoate degradation"
## [5] "Prodigiosin biosynthesis"
## [6] "Microbial metabolism in diverse environments"
## [7] "Dopaminergic synapse"
## [8] "Biofilm formation - Pseudomonas aeruginosa"
## [9] "Galactose metabolism"
## [10] "Antifolate resistance"
## [11] "Biosynthesis of amino acids"
## [12] "cGMP-PKG signaling pathway"
## [13] "Glycolysis / Gluconeogenesis"
## [14] "Bile secretion"
## [15] "Necroptosis"
## [16] "Arachidonic acid metabolism"
## [17] "Small cell lung cancer"
## [18] "Phenazine biosynthesis"
## [19] "Betalain biosynthesis"
## [20] "MAPK signaling pathway - yeast"
## [21] "Aminoacyl-tRNA biosynthesis"
## [22] "Histidine metabolism"
## [23] "Biosynthesis of ansamycins"
## [24] "Apoptosis"
## [25] "Photosynthesis"Actualmente existen esfuerzos que intentan formar bases de datos especializadas en en tratamiento de agua, pero estos esfuerzos se encuentran lejos de sustituir las bases de datos actuales
Predicción de la redundancia funcional
También es posible calcular la redundancia de cada función para cada muestra a través del comando calculateFunctionalRedundancy
Este comando puede ser extremadamente demandante de recursos, por lo que es mejor tomar las precauciones necesarias correspondientes antes de ejecutarlo, como cerrar todos los demás programas y guardar cualquier cambio
calculateFunctionalRedundancy(path_to_otu_table = "~/Documents/Curso_R/data/Ejercicio_otu_table.txt",
path_to_reference_data = PATH,
path_to_temp_folder = "~/Documents/Curso_R/data/Ejercicio_Ref100NR",
database_mode = "Ref100NR",
min_identity_to_reference = 0.97)Resultados
Con la información que se ha generado, es posible hacer diferentes comparación acerca de la diversidad metabólica de nuestras muestras, por ejemplo podemos analizar si existen diferencias en las abundancia de las vías metabólica entre las muestras correspondientes a dos reactores:
O las diferencias entre puntos subsecuentes en el tiempo de cada reactor (circulos R1, diamantes R2):
También es posible enfocarse en vías metabólicas en particular, por ejemplo reacción para la formación de hidrógeno molecular tiene el número de reacción R00019 en KEGG, en la página se puede ver que la enzima responsable tiene el ID 1.12.7.2 correspondiente a los ortólogos K00532, K00533, K00534, K06441, K18016, K18017, K18023, podemos usar esta información para determinar la redundancia funcional de esta vía en las muestras de los reactores 1 y 2.
Al igual que los resultados de la secuenciación masiva del 16S, los resultados de la predicción de perfiles funcionales son extensos por los que se requiere de métodos estadísticos adecuados para obtener la información más relevante.