---
title: 'Anotacion Funcional'
disqus: Miguel Tripp
---
Modulo 6: Anotación funcional: homología de genes, ontología, enriquecimiento y rutas metabolicas
===
## Tabla de contenido
[TOC]
Dentro de nuestro análisis transcriptomico, el siguiente paso es asignar la identidad de cada transcrito, asi como la información biológica del mismo.
La anotación funcional de un gen puede realizarse por las siguientes vias:
+ busqueda de similitud
+ identificación de dominios transmembranales en la secuencia de proteínas
+ ontología génica
+ información de las rutas metabolicas
Para esto, se realiza una busqueda de similitudes entre nuetras secuencias y las secuencias existentes en las bases de datos.
# 6.1 BLAST
BLAST (_Basic Local Alignment Search Tool_) es una herramienta que realiza la busqueda de similitudes entre secuencias biológicas. El programa compara secuencias de nucleótidos o proteínas con las secuencias presentes en las bases de datos, asignandoles una significancia estadística <sup>[[1]](https://blast.ncbi.nlm.nih.gov/Blast.cgi) .
EL tipo de BLAST que se ejecute dependerá del tipo de datos que se van a utilizar.
| Tipo de Blast | Secuencia desconocida (query)| Base de datos |
| -------- | -------- | -------- |
| Blastn | Nucleótidos | Nucleótidos |
|Blastp | Proteínas (AA) | Proteínas|
|Blastx | Nuleótido traducido (AA) | Proteínas |
|tBlastn | Proteínas | Nucleótido traducido |
tBlastx | Nucleótido traducidos | Nucleótidos traducidos |
## 6.1.1 Blastx:
Este proceso es uno de los que ocupa mas tiempo, dependiendo de la base de datos con la que se realiza la anotación, siendo la base `nr` la que ocupa varios días.
### 6.1.1.1 Prueba de Blast en web y Blast en linea de comando
Dentro de nuestra carpeta principal, generamos una nueva carpeta llamada `ANOTACION`. Una vez dentro, creamos un link simbolico de nuestro ensamble.
::: success
Para mantener consistencia con los análisis, vamos a utilizar el mismo ensamble que se utilizó para la expresión diferencial: _good.Trinity.fasta_
:::
```
mkdir ANOTACION
cd ANOTACION
ln -s /LUSTRE/bioinformatica_data/genomica_funcional/Tripp/Curso_2022/ENSAMBLE/good.Trinity.fasta .
```
Ahora, vamos a extraer las primeras 10 secuencias de nuestro ensamble y las importamos a nuetra computadora local
```
awk "/^>/ {n++} n>10 {exit} {print}" good.Trinity.fasta > Trinity.fasta.N10.fasta
```
Ahora entramos a la pagina de [Blast](https://blast.ncbi.nlm.nih.gov/Blast.cgi)
Realiza un Blastx de las secuencias utilizando:
+ non-redundant protein sequences (nr)
+ Swiss-prot
Discute los resultados.
A continuación haremos un Blastx utilizando la linea de comando de estas diez secuencias. Para esto, utilizaremos el siguiente script al que llamaremos `prueba_blast.slrm`:
```
#!/bin/bash
### Directivas
#SBATCH -p cicese
#SBATCH --job-name=Blastx
#SBATCH --ntasks-per-node=24
#SBATCH -N 1
#SBATCH --mem=20GB
#SBATCH --output=blastx-%j.log
#SBATCH --error=blastx-%j.err
#SBATCH -t 06-00:00:00
BLAST=/LUSTRE/bioinformatica_data/genomica_funcional/bin/ncbi-blast-2.4.0+/bin/
DBS=/LUSTRE/bioinformatica_data/genomica_funcional/Tripp/DATABASES/TRINOTATE_DB
cd ${SLURM_SUBMIT_DIR}
$BLAST/blastx -query Trinity.fasta.N10.fasta -db $DBS/uniprot_sprot.pep -num_threads 8 -max_target_seqs 1 -evalue 1E-05 -outfmt 6 > Prueba_Blastx.outfmt6
exit 0
```
La aplicación de Blast en la linea de comando puede generar los resultados en diferentes formatos dependiendo de las aplicaciones que se usen en analisis posteriores (por ejemplo _Blast2GO_ solo acepta resultados en formato fmt5/xml).
El formato fmt consta de seis columnas:
| Encabezado | Contenido |
| -------- | -------- |
| 1. qseqid | query |
| 2. sseqid | referencia |
| 3. pident | Porcentaje de identidad |
| 4. lenght | longitud del alineamiento |
| 5. mismatch | número de mismatches |
| 6. gapopen | número de gaps |
|7. qstart | inicio del alineamiento en query |
| 8. qend | fin del alineamiento en query |
| 9. sstart | inicio del alineamiento en referencía |
| 10. send | fin del alineamiento en referencia |
|11. evalue | valor E |
|12. bitscore | bit score|
---
# 6.2: Anotación con Trinotate
_Trinotate_ es un paquete de anotación diseñado para la anotación funcional automática de un transcriptoma, particularmente de ensambles _de novo_ de organismos modelo y no modelo.
Trinotate utiliza una serie de métodos bien referenciados para la anotación, incluyendo:
* Busqueda de homología en bases de secuencias conocidas (BLAST+/SwissProt).
* Identificación de dominios proteicos (HMMER/PFAM).
* Predicción de peptidos señal y dominios transmembranales (signalP/tmHMM).
* Uso de distintas bases de datos de anotaciones (eggNOG/GO/KEGG)
Trinotate puede realizar la anotación de manera automatica, como se describe [aqui](https://github.com/Trinotate/Trinotate.github.io/wiki/Automated-Execution-of-Trinotate:-Running-computes-and-loading-results).
Para fines de este curso, se realizara la anotación semi-automatica.
El script realiza los sigueintes pasos:
1. **Blastx** usando el ensamble `good.Trinity.fasta`contra la base de datos de **SwissProt**
2. **Transdecorder** para identificar regiones codificantes dentro del ensamble. El proceso consiste en dos pasos:
* Identificación de los ORFs largos con _TransDecoder.LongOrfs_
* Predicción de los ORF que pueden ser codificantes con _TransDecoder.Predict_.
Los detalles del funcionamiento de Transdecoder se pueden encontrar [aqui](https://github.com/TransDecoder/TransDecoder/wiki).
3. __Blastp__ para la busqueda de homología utilizando las secuencias de proteinas identificadas con TransDecoder.
4. __HMMER__ para la dentificación de dominios conservados con la base de datos _Pfam_. Los dominios conservados pueden ser un indicativo de la función de la proteína.
5. __signalP__ y __tmhmm__ para la identificación de peptido señal y dominios transmembranales, respectivamente.
::: warning
Si bien es posible hacer Blast con otras bases de datos (Nr), Trinotate depende fuertemente de las bases de **SwissProt** y **Pfam**
:::
## Paso 1: Homologia y busqueda en bases de datos.
Debido a que la busqueda de homología es un proceso tardado (approx 5 hr para estas secuencias), se usara el siguiente script para realizar la anotación:
```
#!/bin/bash
### Directivas
#SBATCH -p cicese
#SBATCH --job-name=Anotacion
#SBATCH --ntasks-per-node=8
#SBATCH -N 1
#SBATCH --mem=20GB
#SBATCH --output=anotacion-%j.log
#SBATCH --error=anotacion-%j.err
#SBATCH -t 06-00:00:00
BLAST=/LUSTRE/bioinformatica_data/genomica_funcional/bin/ncbi-blast-2.4.0+/bin/
DBS=/LUSTRE/bioinformatica_data/genomica_funcional/Tripp/DATABASES/TRINOTATE_DB
#paqueterias
PATH=$PATH:/LUSTRE/apps/bioinformatica/Trinotate/PerlLib
export PATH=$PATH:/LUSTRE/apps/bioinformatica/Trinotate/auto
export TRINOTATE=/LUSTRE/apps/bioinformatica/Trinotate/
#Trinotate_admin contiene el script para generar el boiler template
export PATH=$PATH:/LUSTRE/apps/bioinformatca/Trinotate/admin
UTILS=/LUSTRE/apps/bioinformatica/trinityrnaseq/util/support_scripts
#programas
export HMMSEARCH_BINARY=$BIN/RNAMMER2_TRINOTATE/hmmer-2.3.bin.intel-linux/binaries/
export PATH=/LUSTRE/apps/bioinformatica/transdecoder/:$PATH
export PATH=$PATH:/LUSTRE/bioinformatica_data/genomica_funcional/bin/signalp-4.1
export PATH=$PATH:/LUSTRE/bioinformatica_data/genomica_funcional/bin/tmhmm-2.0c/bin
export PATH=$PATH:/LUSTRE/bioinformatica_data/genomica_funcional/bin/Trinotate/util/rnammer_support
export PATH=$PATH:/LUSTRE/bioinformatica_data/genomica_funcional/bin/rnammer
export PATH=$PATH:/LUSTRE/bioinformatica_data/genomica_funcional/bin/hmmer-3.1b2-linux-intel-x86_64/binaries
echo "Ejecuntandose la anotacion funcional"
echo "Fecha de inicio: `date`"
v1=`date`
cd ${SLURM_SUBMIT_DIR}
echo "Blastx"
echo "inicio BlastX: `date`"
#Cambiar de nombre el ensamble para evitar nombres de archivos muy largos
mv good.Trinity.fasta Trinity.fasta
$BLAST/blastx -query Trinity.fasta -db $DBS/uniprot_sprot.pep -num_threads 8 -max_target_seqs 1 -evalue 1E-05 -outfmt 6 > Blastx_swissprot.outfmt6
echo "Resultados de Blastx se guardaron como BlastX_swissprot.outfmt6"
#relacion de isoformas y genes
echo "generacion de gene_trans_map"
$UTILS/get_Trinity_gene_to_trans_map.pl Trinity.fasta > Trinity.gene_trans_map
#identificacion de regiones codificantes
echo "Identificacion de regiones codificantes con transdecoder"
echo "inicio Transdecoder: `date`"
TransDecoder.LongOrfs -t Trinity.fasta --gene_trans_map Trinity.gene_trans_map > transdecoder.log
#Predecir probables ORF codificantes
echo "Predecir probables ORF codificantes: `date`"
TransDecoder.Predict -t Trinity.fasta --cpu 8 >> transdecoder.log
#BlastP
echo "Inicio de BlastP: `date`"
$BLAST/blastp -query Trinity.fasta.transdecoder.pep -db $DBS/uniprot_sprot.pep -num_threads 8 -max_target_seqs 1 -outfmt 6 -evalue 1e-5 > Blastp_swissprot.outfmt6
echo "se genero el archivo BlastP_swissprot.outfmt6 `date`"
#Identificacion de dominios proteicos con HMMSCAN
hmmscan --cpu 8 --domtblout TrinotatePFAM.out $DBS/Pfam-A.hmm Trinity.fasta.transdecoder.pep
echo "hmmscan `date`"
#Prediccion de senhales peptidicas
signalp -f short -n signalp.out Trinity.fasta.transdecoder.pep
#Prediccion de dominios transmembranales
tmhmm --short < Trinity.fasta.transdecoder.pep > tmhmm.out
echo "tmhmm.out: `date`"
###final e imprimir fecha###
echo "Tiempo de inicio: $v1 "
echo "Tiempo final: `date`"
exit 0
```
## Paso 2: Preparar y generar el reporte de la anotación
En principio, Trinotate hace uso de una base de datos SQLite las cual permite una busqueda rapida y eficiente de un térmico especifico.
Para generar el reporte de la anotación con Trinotate, primero es necesario cargar todos los resultados bioinformaticos dentro de esta base SQLite. Trinotate proporciona un molde precargado con información relacionada con __SwissProt__ y __Pfam__.
::: info
Para fines de este curso, ya se tiene preparado un _boilerplate_ el cual será necesario copiar a nuestro directorio. Sin embargo, este se puede prepara facilmente siguiente las instrucciones descritas [aqui](https://github.com/RJEGR/Transcriptomics/blob/master/markdown/trinotate.md)
:::
Vamos a realizar una copia del _boilerplate_ con el siguiente ruta:
```
cp /LUSTRE/bioinformatica_data/genomica_funcional/Tripp/DATABASES/TRINOTATE_DB/Trinotate.sqlite ./
```
Ahora vamos a generar un archivo llamado `exports_Trinotate`el cual incluye algunas rutas de los programas y utilizades que se van a requerir. Estos Ya deberian estar contenidos en sus respectivos _bash_profile_, pero se usaran en caso de que alguno tenga errores.
```
export PATH=$PATH:/LUSTRE/apps/bioinformatica/Trinotate/PerlLib
export PATH=$PATH:/LUSTRE/apps/bioinformatica/Trinotate/auto
export PATH=$PATH:/LUSTRE/apps/bioinformatica/Trinotate/
#Trinotate_admin contiene el script para generar el boiler template
export PATH=$PATH:/LUSTRE/apps/bioinformatica/Trinotate/admin
export PATH=$PATH:/LUSTRE/apps/bioinformatica/Trinotate/util
```
Y ahora cargamos estas rutas con `source exports_Trinotate`
Iniciamo el _boilerplate_
```
Trinotate Trinotate.sqlite init --gene_trans_map Trinity.gene_trans_map --transcript_fasta Trinity.fasta --transdecoder_pep Trinity.fasta.transdecoder.pep
```
Ahora cargamos las bases de datos de __BLastx__ y __Blastp__:
```
Trinotate Trinotate.sqlite LOAD_swissprot_blastx Blastx_swissprot.outfmt6
Trinotate Trinotate.sqlite LOAD_swissprot_blastp Blastp_swissprot.outfmt6
```
Ahora cargamos los resultados de __Pfam__, __signalP__ y __tmhmm__
```
Trinotate Trinotate.sqlite LOAD_pfam TrinotatePFAM.out
Trinotate Trinotate.sqlite LOAD_signalp signalp.out
Trinotate Trinotate.sqlite LOAD_tmhmm tmhmm.out
```
Y finalmente generamos el reporte, al que llamaremos `Trinotate_report.xls`:
```
Trinotate Trinotate.sqlite report > Trinotate_report.xls
```
Finalmente, vamos a extraer la anotación de Gene Ontology de cada transcrito ejecutando la siguiente linea
```
extract_GO_assignments_from_Trinotate_xls.pl --Trinotate_xls Trinotate_report.xls -G --include_ancestral_terms > Trinotate_GO_annotations.txt
```
Mueve reporte y las anotaciones GO a tu computadora personal para analizarlo utilizando `TrinotateR`en R.
# 6.3 Enriquecimiento funcional con DAVID
El enriquecimiento funcional consiste en aplicar una prueba estadística para verificar si los genes de interes estan mas asociados a una cierta función biológica que lo que se esperaria en un set de genes aleatorios. Para ello se recurre a diferentes bases de datos como GO o KEGG, para obtener la información ontológica-funcional relacionada a dichos genes candidatos.
De manera general, para realizar un análisis de enriquecimiento funcional es necesario contar con:
* Un set de genes de interes (e.g. expresados diferencialmente): **_study set_**
* Un set con todos los genes anotados son considerados en el análisis y que deben contener el set de estudio: **_background_**
## 6.3.1: Extracción de los sets de genes
Dentro de la carpeta `ANOTACION` crea una nueva carpeta llamada `DAVID_SUBSET`. Dentro de esta carpeta, copia los transcritos que fueron expresados diferencialmente en cada una de las comparaciones (Ctrl-T1; Ctrl-T2; Ctrl-T3).
```
cp ../EXP_DIF/diffExpr/*Ctrl_vs*.DE.subset ./
```
Posteriormente, vamos a extraer la primera columna (ID del transcrito)
```
cut RSEM.isoform.counts.matrix.Ctrl_vs_T1.edgeR.DE_results.P0.01_C1.DE.subset -f 1 > Ctrl_vs_T1.txt
cut RSEM.isoform.counts.matrix.Ctrl_vs_T2.edgeR.DE_results.P0.01_C1.DE.subset -f 1 > Ctrl_vs_T2.txt
cut RSEM.isoform.counts.matrix.Ctrl_vs_T3.edgeR.DE_results.P0.01_C1.DE.subset -f 1 > Ctrl_vs_T3.txt
```
Ya que tenemos la lista de todos los transcritos sobre expresados, vamos a filtrar los resultados de Blastx de la anotación. Para esto, vamos a copiar nuestro archivo `Blastx_swissprot.outfmt6` y ejecutamos el siguiente comando para filtrar los transcritos DE anotados y quedarnos con la columna con el `uniprot_ID`:
```
grep -Fwf Ctrl_vs_T1.txt Blastx_swissprot.outfmt6 | cut -f 2 > Ctrl_vs_T1_Blasted.txt
grep -Fwf Ctrl_vs_T2.txt Blastx_swissprot.outfmt6 | cut -f 2 > Ctrl_vs_T2_Blasted.txt
grep -Fwf Ctrl_vs_T3.txt Blastx_swissprot.outfmt6 | cut -f 2 > Ctrl_vs_T3_Blasted.txt
```
Este archivo `Blasted_DE.txt`constituye nuestro set de genes de interese **_study set_**. Ahora necesitamos la lisa de genes **_background_**.
```
cut -f 2 Blastx_swissprot.outfmt6 > Blast_background.txt
```
Descargamos ambos archivos a nuestra computadora local.
:::info
Puedes utilizar el comando anterior para filtrar el reporte de Trinotate_report.xls con los transcritos que fueron expresados diferencialmente. Por ejemplo:
```
grep -Fwf Ctrl_vs_T1.txt ../Trinotate_report.xls > Ctrl_vs_T1_DE_Anotados.txt
```
:::
***
## 6.3.2 Análisis en DAVID
DAVID (Base de datos para Anotación, Visualización y Descubrimiento Integrado[Database for Annotation, Visualization and Integrated Discover])
1. Para analizar nuestros datos, entramos a la pagina de [DAVID](https://david.ncifcrf.gov/)
2. Posteriormente, ingresamos a la opcion `start analysis (1)`
3. Ingresamos nuestro archivo con nuestras secuencias anotadas que fueron sobre-expresadas `(2)`
4. Seleccionamos el tipo de identificador, que en este caso es `uniprot id (3)`
5. Seleccionamos que esta lista corresponde a un `gene list (4)` y presionamos `submit list`
6. Posteriormente, repetimos los pasos del 2 al 4 pero seleccionado la lista de secunecias anotadas del background.
**Anotación funcional (_Functional Annotation Chart_)**
En el panel de resultados, selecciona la opcion `Functional Annotation Chart`. Esto abrira una nueva ventana llamada **Annotation Summary Results**. La tabla muestra una lista de los terminos que fueron enriquecidos estadisticamente basados en el set de genes background. Esta tabla puede descargarse en un formato compatible con Excel.
Las columnas _Count_ y _%_ indican el número de genes en la lista anotados con ese término y su porcentage.
**Functional Annotation Clustering**
Dentro de esta opción, en lugar de ver la lista de terminos individuales, podras observar grupos (_clusters_) de terminos que estan relacionados entre ellos.
Para ver un tutorial mas detallado, visita [aqui](https://biochem.slu.edu/bchm628/handouts/2015/DavidTutorial_2015.pdf)
:::info
**Filtros en la tabla de anotación:**
* **Count threshold:** Número minimo de cuenta de genes correspondientes a un determinado termino de anotación (_annotation term_).
* **EASE Score threshold:** Valor P de la prueba exacta de Fisher modificado para enriquecimiento. Va desde 0 a 1, donde un valor P = 0 representa un enriquecimiento perfecto.
* **Fold enrichment:** Definido como el ratio entre dos proporciones. Por ejemplo, si 40/400 (10%) de tus genes estan involucrados en un determinado procesos (eg. "actividad quinasas") and la información del _background_ es 300/30000 (1%), entonces tendras un 10%/1% = 10 fold enrichment.
* **Benjamin-Hochberg, Bonferroni, FDR:** Corrección estadistica para comparaciones multiples. Corrigen el valor P para hacerlo mas conservador y reducir el número de falsos positivos.
* **LT (_list total_):** Número de genes en tu lista que fueron anotados con cualquier termino
* **PH (_population hits_):** Número de genes con este termino GO en el universo _background_
* **PT (_population total_):** Número de genes totales en el universo (_background_)
:::
Modulo 7: Análisis de redes de asociación de proteías
===
[String](https://string-db.org/)
Las proteínas y sus interacciones funcionales constituyen la medula de la maquinaria celular. Sus redes de conectividad requieren ser analizadas para lograr un entendimiento mas completo de un fenómeno biológico [[1]](https://academic.oup.com/nar/article/47/D1/D607/5198476.)
STRING es una base de datos interacciones proteina-protein conocidas o predichas. La interacción puede ser mediante asociación directa (física) o indirecta (funcional).
## 7.1: ¿Para que usar String?
* _Predicción de la función de un gen_: Examinar genes (proteínas) en el contexto de una red muestra conecciones en un conjunto de genes/proteínas involucradas en el mismo proceso biológico.
* _Identificación de enfermedades/indicadores:_ Mediante el análisis de genes/proteínas dentro de una red de asociación, es posible identificar una red o sub-red que pueda resultar como indicador de la respuesta ante una enfermedad o condición célular/ambiental.
## 7.2: Evidencias utilizadas en String
Las interacciones entre proteínas mostradas en String provienen de cinco fuentes principales:
La base de datos de string cubre 24,584,628 proteínas de 5,090 organismos
## 7.3: Análisis con String
Para ver un tutorial completo, visita [aqui](https://string-db.org/cgi/help.pl?sessionId=LvoIktpzlX2Q)
### 7.3.1: Asociación usando una lista de genes de ejemplo
La siguiente lista corresponde a genes que se sobre expresaron en un molluco bentónico tras una exposición crónica a temperaturas elevadas.
1. Abre los [datos](https://drive.google.com/open?id=1Y2XjKeZKkk3lLoMPWwJykDiz-jZLFaKw&authuser=1) y copia la lista.
2. Abre [string](https://string-db.org/) dentro de la pestaña de _Multiple proteins_ y pega la lista de genes
3. Selecciona _Homo sapies_
4. En el siguiente paso, no hagas cambios en la lista de genes y solo presiona _Continue_.
5. Dentro de la pestaña de _Settings_ selecciona un nivel de confianza alta (0.70) y selecciona la opción de ocultar nodos desconectados.
6. Dentro de la pestaña _Analysis_ identifica que procesos biológicos fueron enriquecidos en esta red.
### 7.3.2: Asociación usando la lista de genes expresados diferencialmente con **Ctrl-T1; Ctrl-T2; Ctrl-T3**.
1 . Para poder utilizar nuestras listas de genes expresados diferencialmente en String, debemos remover el identificador del taxa para cada anotación, esto se logra ejecutando el siguiente script:
```
awk -F"_" '$1=$1' OFS="\t" Ctrl_vs_T1_Blasted.txt | cut -f 1 > Ctrl_vs_T1_Blasted_String.txt
awk -F"_" '$1=$1' OFS="\t" Ctrl_vs_T2_Blasted.txt | cut -f 1 > Ctrl_vs_T2_Blasted_String.txt
awk -F"_" '$1=$1' OFS="\t" Ctrl_vs_T3_Blasted.txt | cut -f 1 > Ctrl_vs_T3_Blasted_String.txt
```
Una vez que tengamos nuestras listas, podemos importarlas a nuestra computadora local o copiar la lista, para posteriormente pegarla en String.
2. Repite los pasos de 1 a 6 de la sección 7.3.1.
###### tags: `Clase NGS`
Sign in with Wallet
Connect another wallet