-
Notifications
You must be signed in to change notification settings - Fork 0
/
Bioinformatica_SyT8087_MUCsS.Rmd
371 lines (265 loc) · 18.7 KB
/
Bioinformatica_SyT8087_MUCsS.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
---
title: |
| BioInformática (Salud y Tecnología)
| Máster Universitario en Ciencias de la Salud
author: |
| **Antonio Canepa, Ph.D.**
| *[email](mailto:[email protected])* // [twitter](https://twitter.com/OMDataScience)
date: "`r format(Sys.time(), '%B, %Y')`"
output:
html_document:
df_print: paged
toc: yes
toc_float: yes
always_allow_html: true
---
```{r, echo=FALSE, out.width="20%", fig.align = "right", fig.cap = '', warning=FALSE, message=FALSE}
knitr::include_graphics("Images/ubu.jpg")
```
# 1.- Introducción
Esta es una actividad complementaria, pero altamente recomendable a seguir ya que nos introduce en el mundo de la bioinformática con códigos muy sencillos y usando herramientas disponibles para el análisis de secuencias que nos facilitarán mucho la vida.
Comenzaremos con el trabajo sobre una secuencia "ficticia" de ADN, recordando de dónde provienen tenemos el siguiente esquema:
```{r, echo=FALSE, out.width="90%", fig.align = "center", fig.cap = 'Esquema transcripción', warning=FALSE, message=FALSE}
knitr::include_graphics("Images/Bioinformatica_Schema.png")
```
Por lo que lo primero que veremos cuando tengamos una lectura de ADN, será una secuencia nucleotídica del tipo:
```{r}
seq_dna <- c("ATGTCACCACAAACAGAGACT")
seq_dna
```
Vemos que nuestro ejemplo de cadena `seq_dna` no es más que un conjunto de letras. Dentro de R (y otros lenguajes de programación) a esto se le conoce como **character string** o **cadena de caracteres**. Para entender cómo lo está leyendo R, le pediremos que nos muestre su *estructura* usando la función `str()`.
```{r}
str(seq_dna)
```
Las tres letras señaladas **chr** hacen referencia a que es un *character*, es decir una cadena de caracteres que contienen los nucleótidos de nuestra secuencia.
Como primer ejercicio de programación, vamos a crear una función muy sencilla llamada `contar_nucleotidos()`, que nos permitirá, como dice su nombre, contar el número de cada uno de los nucleótidos que tenemos en nuestra secuencia:
```{r}
contar_nucleotidos <- function(secuencia_dna = seq_dna){
conteo_nucleotidos = unlist(strsplit(x = secuencia_dna, split = "", fixed = TRUE))
print(table(conteo_nucleotidos))
}
```
Cuando ejecutamos la caja de código superior vemos que nada sucede y esto es correcto, no vemos que suceda nada porque cuando creamos una función no estamos solicitando un resultado, solamente estamos creando una **receta** para trabajar con nuestros datos.
Tampoco importa si no entendemos al 100&% la función que estamos creando. Comentaros que el paso fundamental de la función es el uso de `unlist()` y de `strsplit()`.
Ahora sí, con nuestra función recién creada vamos a "*contar*" cuántos nucleótidos de cada tipo hay:
```{r}
contar_nucleotidos(secuencia_dna = seq_dna)
```
A continuación, revisaremos algunos de los pasos más fundamentales en *BioInformática* y para ello nos ayudaremos de funciones que ya existen (¡no tendremos que programarlas nosotros! ;P ).
# 2.- Cadena complementaria
Para poder generar una cadena complementaria de ADN, deberemos seguir las reglas de apareamiento de bases nitrógenedas, donde las correspondencias son *Adenina-Timina* (**`A-T`**) y *Citosina-Guanina* (**`C-G`**). Por lo tanto, cada vez que un nucleótido de nuestra cadena sea `A`, en la cadena complementaria correspondiente, el nucleótido que le tocaría sería el `T`.
Podemos crear una función que realice esta acción a la que llamaremos `Complementaria()`. No es problema si no entiendes lo que sucede dentro de la función, pero es recomendable que intentes *mirar las tripas* de la función:
```{r}
Complementaria <- function(secuencia_dna = seq_dna){
complementaria_adn = chartr("ATGC","TACG", secuencia_dna)
complementaria_adn = paste(complementaria_adn, collapse = "")
return(complementaria_adn)
}
```
Para ver si la función `Complementaria()` está haciendo lo que toca, imprimeremos ambas secuencias:
```{r}
seq_dna
Complementaria(seq_dna)
```
Ahora bien, para facilitarnos la vida, usaremos el paquete de R llamado `bioseq`, que como su nombre nos indica es *una caja de herramientas para manipular secuencias biológicas*. Para ver el paquete clicar [aquí](https://cran.r-project.org/web/packages/bioseq/index.html) y para la publicación del mismo, clicar [aquí](https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.13490).
Para ello necesitaremos instalarlo (descargarlo desde la web) usando la función `install.packages()`, y una vez descargado, para poder usarlo debemos *activarlo*/*encenderlo*, con la función `library()`.
```{r warning=FALSE, message=FALSE}
install.packages("bioseq")
library("bioseq")
```
Sin entrar en conceptos específicos de programación, lo primero que deberemos hacer será crear un objeto llamado `x_dna` que usando la función `dna()` y nuestra secuencia de nucleótidos (`seq_dna`), creará un objeto de clase `"bioseq_dna"`.
```{r}
# library(bioseq)
x_dna <- dna(seq_dna)
x_dna
```
Al imprimir el objeto `xdna` vemos que *luce diferente* y esta es una de las particularidades y ayudas que representa trabajar con paquetes específicos para trabajos bioinformáticos, en este caso.
Lo importante es que tenemos una sola secuencia y los nucleótidos se pintan de un color específico cada uno.
Si te interesa mirar y revisar en detalle la ayuda de este paquete la información la puedes encontrar [aquí](https://cran.r-project.org/web/packages/bioseq/vignettes/intro-bioseq.html).
Entonces, para poder crear la cadena complementaria vamos a usar la función `seq_complement()` a la que le pasamos nuestro objeto `x_dna`. Para una mejor comparación, vamos a imprimir en pantalla ambas secuencias:
```{r}
x_dna
seq_complement(x_dna)
```
Como podemos ver, hemos sido capaces de recrear la cadena complementaria y gracias al uso de funciones *ad-hoc* para este tipo de análisis el resultado luce bastante mejor, aunque en escencia es lo mismo que hemos creado a mano más arriba.
# 3.- Transcripción
Cuando vamos de camino a generar una proteína, necesitamos sacar la información desde dentro del núcleo al exterior y para ello usamos el ARN mensajero `ARNm`. Es decir, necesitamos **transcribir** la información desde la cadena de ADN a una cadena de ARN.
```{r, echo=FALSE, out.width="90%", fig.align = "center", fig.cap = 'Detalle del proceso de la transcripción', warning=FALSE, message=FALSE}
knitr::include_graphics("Images/Transcripcion_es_0.jpg")
```
Como ya sabréis las correspondencias en este caso cambian y tenemos que la *Timina* es reemplazada por el *Uracilo*, por lo que las parejas ahora son **Adenina-Uracilo* (**`A-U`**) y *Citosina-Guanina* (**`C-G`**).
En términos de programación, lo único que tenemos que hacer es repetir el paso de cambiar las letras, en este caso solo intercambiando la T por la U:
```{r}
Transcripcion <- function(secuencia_dna = seq_dna){
seq_arn = chartr("T","U", secuencia_dna)
seq_arn = paste(seq_arn, collapse = "")
return(seq_arn)
}
```
Para comprobar que nuestro código está funcionando, volveremos a imprimir la secuencia de ADN original (objeto `seq_dna`) y el resultante de aplicar nuestra función `Transcripcion()` sobre este ADN original:
```{r}
print(seq_dna)
Transcripcion(seq_dna)
```
Otra manera, más sencilla si se quiere, es el uso de la función `seq_transcribe()` del paquete `bioseq`. Fundamentalmente hace lo mismo que nuestra función, pero resaltando los resultados de una manera mucho más visual. En este caso el objeto que deberemos pasar es el objeto `x_dna` que posee la clase `"bioseq_dna"`. Tal que:
```{r}
x_rna <- seq_transcribe(x_dna)
x_rna
```
Para evaluar los cambios imprimiremos en pantalla tanto el ADN como el ARN.
```{r}
x_dna
x_rna
```
# 4.- Traducir secuencias en aminoácidos
Estando ahora fuera del núcleo, la información transportada como `ARNm` es procesada por el ribosoma y el ARN d etransfeencia `ARNt` para generar una cadena de polipéptidos, que dará a lugar a un proteína que será funcional. Este poceso denominado *traducción* puede verse en el siguiente esquema:
```{r, echo=FALSE, out.width="90%", fig.align = "center", fig.cap = 'Esquema de la traducción', warning=FALSE, message=FALSE}
knitr::include_graphics("Images/Traduccion_es_0.jpg")
```
Los ribosomas y el `ARNt` solo pueden leer trozos o **codones** de `ARNm` que consisten en 3 aminoácidos. Estos tres aminoácidos son lueg convertidos a polipéptidos según la siguiente tabla de códigos genéticos
```{r, echo=FALSE, out.width="90%", fig.align = "center", fig.cap = 'Tabla de Aminoácidos según ARN', warning=FALSE, message=FALSE}
knitr::include_graphics("Images/Codigo-genetico_es_0.jpg")
```
Como ejemplo, usaremos la función `seq_translate()` del paquete `bioseq` para traducir una secuencia de ARNm en Proteínas.
```{r}
seq_translate(x_rna)
```
Veremos con solo un ejemplo más, que las funciones de este paquete `bioseq` nos permiten traducir a proteína pasándole directamente unas cadenas de AND, realizando internamente y sin que veamos el esultado, la transformación a ARNm y su traducción, tal que:
```{r}
x <- dna(c("ATGCAGA", "GGR","TTGCCTAGKTGAACC", "AGGNGC", "NNN"))
seq_translate(x)
```
### 4.1 Visualización de proteínas
Por último, revisaremos algunas heramientas para visualizar proteínas, ya que gracias a los nuevos algoritmos, procesadores y lenguajes de programación esta tarea resulta ahora mucho más sencilla.
La herramienta que revisaremos es un un **paquete** de R llamado "*NGLVieweR*" y su repositorio en github puede visitarse clicando [aquí](https://github.com/nvelden/NGLVieweR/).
```{r message=FALSE, warning=FALSE}
install.packages("NGLVieweR")
library("NGLVieweR")
```
Para visualizar las proteínas necesitamos su código del "**Protein Data Bank**" o "**Banco de datos de proteínas**", disponible en [https://www.rcsb.org/](https://www.rcsb.org/).
Buscaremos a modo de ejemplo la proteína **7CID**, catalogada como una estructura cristalina de *Pseudomonas aeruginosa* "**LpxC**" en complejo con inhibidor PAO1. Puede visitarse clicando [aquí](https://www.rcsb.org/structure/7CID).
```{r}
#Load protein by PDB code
NGLVieweR("7CID") %>%
addRepresentation("cartoon")
```
Este paquete también posee su propia aplicación **shiny** ("*shinyNGLVieweR*"), que se programa en su totalidad usando el lenguaje de programación R. El repositorio en GitHub de esta *APP* puede visitarse clicando [aquí](https://github.com/nvelden/shinyNGLVieweR); o bien para acceder e interactuar con la **APP Shiny**, deberemos acceder al siguiente [enlace](https://niels-van-der-velden.shinyapps.io/shinyNGLVieweR/).
Como tarea deberéis interactuar con la aplicación Shiny ("**shinyNGLVieweR**") y probar a cargar la estructura de la glicoproteína espiga (*spike*) del SARS-CoV-2 , usando el código del PBD: **6VXX**. Puedes ver esta protéina directament en la página del [Protein Data Bank ](https://www.rcsb.org/structure/6vxx).
Acá os dejo una animación de esta glicoproteína.
```{r message=TRUE, warning=FALSE, echo=FALSE}
NGLVieweR("6VXX") %>%
stageParameters(backgroundColor = "white", zoomSpeed = 1) %>%
addRepresentation("cartoon",
param = list(name = "cartoon", colorScheme = "residueindex")
) %>%
setSpin()
```
# 5.- Caso Práctico
Este caso práctico está basado en una de las ayudas del paquete `bioseq` que se puede visitar [aquí](https://cran.r-project.org/web/packages/bioseq/vignettes/ref_database.html).
```{r warning=FALSE, message=FALSE}
library(bioseq)
library(tidyverse)
```
Los datos genéticos suelen entregarse en formato FASTA. En este ejemplo utilizaremos secuencias `rbcL` de especies de *Fragilaria* recuperadas del NCBI y disponibles directamente como FASTA sin analizar en el paquete. Podemos utilizar la función `read_fasta()` para leer directamente un archivo FASTA en R. Aquí leemos un vector de caracteres ya disponible en R, pero la mayoría de las veces el argumento file se utilizará para pasar una ruta a un archivo disponible localmente o a través de la red.
```{r}
data(fragilaria, package = "bioseq")
fra_data <- read_fasta(fragilaria)
fra_data
```
Según el resultado, tenemos un vector de ADN con 41 secuencias. Estas secuencias se obtuvieron directamente del NCBI. Por lo tanto, no están alineadas y tienen diferentes longitudes. Es fácil obtener el rango de longitudes con el siguiente comando:
```{r}
seq_nchar(fra_data) %>% range()
```
A continuación, utilizaremos la función as_tibble para convertir el vector ADN en un tibble con dos columnas (una para los nombres y otra para las secuencias). También limpiaremos los nombres y los dividiremos en dos columnas. Finalmente creamos una columna extra con la longitud (número de caracteres) de cada secuencia
```{r}
fra_data <- tibble(label = names(fra_data), sequence = fra_data)
fra_data <- fra_data %>%
mutate(genbank_id = str_extract(label, "([^\\s]+)"),
taxa = str_extract(label, "(?<= ).*")) %>%
select(genbank_id, taxa, sequence)
fra_data <- fra_data %>%
mutate(n_base = seq_nchar(sequence))
fra_data
```
### 5.1 Recorte de secuencias
En este ejemplo sólo estamos interesados en una región concreta de las secuencias, un fragmento corto utilizado en muchos estudios de metabarcodificación de diatomeas. Para extraer la región de interés podemos proporcionar la posición de dos nucleótidos, pero esta solución requiere que las secuencias estén alineadas. La solución alternativa consiste en extraer utilizando patrones de nucleótidos. En este caso, utilizaremos los cebadores directo e inverso que delimitan la región del código de barras para recortar las secuencias. En primer lugar, creamos dos vectores de ADN, uno para los cebadores directo e inverso:
```{r}
FWD <- dna("AGGTGAAGTAAAAGGTTCWTACTTAAA",
"AGGTGAAGTTAAAGGTTCWTAYTTAAA",
"AGGTGAAACTAAAGGTTCWTACTTAAA")
REV <- dna("CAGTWGTWGGTAAATTAGAAGG",
"CTGTTGTWGGTAAATTAGAAGG")
```
Se habrá dado cuenta de que algunos caracteres de los cebadores son ambiguos, lo que significa que pueden representar varios nucleótidos. Para obtener todas las combinaciones representadas por esas secuencias ambiguas de cebadores podemos utilizar la función seq_disambiguate_IUPAC(). Por ejemplo, para los cebadores delanteros:
```{r}
seq_disambiguate_IUPAC(FWD)
```
Podemos ver que las tres secuencias ambiguas utilizadas como cebadores delanteros representan en realidad ocho secuencias diferentes. Afortunadamente, las funciones de bioseq pueden reconocer automáticamente los nucleótidos ambiguos e interpretarlos como tales. Por ello, siempre que sea posible, se recomienda construir vectores de ADN/ARN/AA en lugar de vectores de caracteres simples.
Ahora recortamos las secuencias utilizando los cebadores directo e inverso como delimitadores. Podemos esperar que las secuencias devueltas se alineen porque corresponden a una región muy específica con una longitud fija.
```{r}
fra_data <- fra_data %>%
mutate(barcode = seq_crop_pattern(sequence,
pattern_in = list(FWD),
pattern_out = list(REV)))
fra_data
```
Como era de esperar, las secuencias del código de barras parecen estar alineadas. También podemos ver que esta operación ha producido varios NA que corresponden a secuencias en las que no se detectaron los cebadores, probablemente porque la región del código de barras no estaba incluida en la secuencia. Podríamos simplemente excluir los valores NA, pero es más seguro filtrar sólo las secuencias cuya longitud coincide exactamente con la longitud esperada del código de barras (es decir, 312 pb).
```{r}
fra_data <- fra_data %>%
filter(seq_nchar(barcode) == 312)
fra_data
```
### 5.2 Secuencias de consenso y filogenia
Ahora queremos reconstruir un árbol filogenético para las secuencias. Podemos incluir todas las secuencias en el árbol, pero para mayor claridad queremos mantener sólo una secuencia por taxón. En lugar de seleccionar una secuencia al azar, calculamos una secuencia de consenso para cada taxón utilizando la función seq_consensus(). Hay diferentes métodos para calcular una secuencia de consenso implementados en la función, puede consultar el manual para más detalles. El método por defecto utilizado a continuación utiliza la regla de la mayoría.
```{r}
fra_consensus <- fra_data %>%
group_by(taxa) %>%
summarise(consensus_barcode = seq_consensus(barcode))
fra_consensus
```
A continuación, utilizaremos ape para hacer una rápida reconstrucción filogenética de estas secuencias. Podemos utilizar la función as_DNAbin para convertir un tibble en un objeto DNAbin y luego utilizar las funciones de ape para reconstruir y trazar un árbol filogenético con el algoritmo BIONJ basado en distancias.
```{r message=FALSE, warning=FALSE}
fra_consensus %>%
as_DNAbin(consensus_barcode, taxa) %>%
ape::dist.dna() %>%
ape::bionj() %>%
plot()
```
### 5.3 Agrupación de secuencias
El árbol anterior sugiere que varias secuencias son exactamente iguales. Esto puede evaluarse rápidamente:
```{r warning=FALSE, message=FALSE}
duplicated(fra_consensus$consensus_barcode)
```
De hecho, varias secuencias están duplicadas y es posible que queramos conservar sólo las secuencias únicas. Podríamos seleccionar sólo las secuencias únicas, pero en este caso utilizaremos la función seq_cluster para agrupar las secuencias en función de su similitud.
```{r}
fra_consensus <-
fra_consensus %>%
mutate(cluster = seq_cluster(consensus_barcode,
threshold = 0.001))
fra_consensus
```
Por último, volvemos a calcular una secuencia consenso para cada grupo y reconstruimos el árbol filogenético.
```{r warning=FALSE, message=FALSE}
fra_consensus <-
fra_consensus %>%
group_by(cluster) %>%
summarise(taxa_group = paste(taxa, collapse = "/"),
consensus_barcode = seq_consensus(consensus_barcode))
fra_consensus
```
```{r message=FALSE, warning=FALSE}
fra_consensus %>%
as_DNAbin(consensus_barcode, taxa_group) %>%
ape::dist.dna() %>%
ape::bionj() %>%
plot()
```
Es importante señalar que este árbol no es una reconstrucción filogenética exacta, pero puede ayudar a visualizar las relaciones jerárquicas entre las secuencias. Esto es especialmente útil con conjuntos de datos más grandes para detectar clados polifiléticos no resueltos por la región del código de barras.
### 5.4 Exportar datos
Cuando estamos contentos con los resultados, normalmente queremos exportar los datos. En este caso, podríamos exportar directamente el marco de datos utilizando write_csv(). Sin embargo, podría ser más conveniente exportar secuencias en FASTA para abrirlas en otro software más tarde. Esto se puede conseguir utilizando la función write_fasta().
```{r warning=FALSE, message=FALSE}
fra_consensus %>%
select(taxa_group, consensus_barcode) %>%
deframe() %>%
write_fasta("my_sequences.fasta")
```