Introducción a Bioconductor - Biostrings

Bioconductor es simplemente una serie de paquetes para hacer varios tipos de análisis de relevancia biológica dentro de R. Pueden leer más sobre este proyecto en: www.bioconductor.org

En esta sesión veremos un paquete de Biocondudctor dedicado a la manipulación de secuencias biológicas. Biostrings está diseñado para poder trabajar con cadenas de caracteres muy largos (secuencias de DNA o RNA) de una manera sencilla y rápida.

Instalando Bioconductor y Biostrings

Es posible que en sus computadoras aún no esté instalado Bioconductor. Para verificarlo pueden hacer lo siguiente:

library(Biobase)

Si les marca un error, deben terminar el resto de esta sección para instalarlo.

Para instalar Bioconductor, deben hacerlo dentro de R. Pueden consultar la descripción completa del proceso de instalación en la página de Bioconductor: install

A partir de la versión 3.5.0 de R, la instalación de Bioconductor cambió. Si tienen una versión anterior, lo siguiente debe bastar:

source("http://bioconductor.org/biocLite.R")
biocLite()

biocLite("Biostrings")  # Para instalar un paquete en particular

Para las versiones más nuevas de R (de la 3.5.0 en adelante), pueden instalar Bioconductor de la siguiente manera:

if (!requireNamespace("BiocManager"))
    install.packages("BiocManager")
BiocManager::install()

BiocManager::install("Biostrings")  # Para instalar un paquete en particular

Esto puede tomar algo de tiempo, dependiendo de la velocidad de conexión al internet.

La lista completa de paquetes que pueden instalar de la versión más reciente de Bioconductor está disponible en: packages

Cabe mencionar que las versiones de Bioconductor están forzosamente ligadas a ciertas versiones de R. Por lo tanto, para tener la versión más nueva de un paquete de Bioconductor hay que primero instalar la versión más nueva de R.

Manipulación de secuencias con Biostrings

Empecemos a ver que nos ofrece Biostrings. Carguemos el paquete, incluyendo el paquete base de Bioconductor:

library(Biobase)
library(Biostrings)

Biostrings define una nueva clase de objeto XString, el cual puede ser de tipo BString (cadenas genéricas), DNAString, RNAString o AAString. Veamos como trabajar con un DNAString.

dnaSeq = DNAString("TTCAGATCTAGTTCGTGTGTGACTGATGATCTGTCACACGTTTTTCTGATCTTCTGACTAGTCGAT")
dnaSeq

##   66-letter "DNAString" instance
## seq: TTCAGATCTAGTTCGTGTGTGACTGATGATCTGTCACACGTTTTTCTGATCTTCTGACTAGTCGAT

class(dnaSeq)

## [1] "DNAString"
## attr(,"package")
## [1] "Biostrings"

Podemos accesar a parte de la secuencia utilizando [ ]:

dnaSeq[4]

##   1-letter "DNAString" instance
## seq: A

dnaSeq[4:10]

##   7-letter "DNAString" instance
## seq: AGATCTA

Existen muchas funciones definidas (o re-definidas) para trabajar con objetos tipo XString:

length(dnaSeq)

## [1] 66

reverseComplement(dnaSeq)

##   66-letter "DNAString" instance
## seq: ATCGACTAGTCAGAAGATCAGAAAAACGTGTGACAGATCATCAGTCACACACGAACTAGATCTGAA

alphabetFrequency(dnaSeq)

##  A  C  G  T  M  R  W  S  Y  K  V  H  D  B  N  -  +  . 
## 12 13 14 27  0  0  0  0  0  0  0  0  0  0  0  0  0  0

También podemos leer archivos FASTA directamente. Para este ejercicio pueden descargar un archivo que contiene las secuencias de algunos genes de la bacteria E. coli: ecoliORFs.fa

Procuren guardar el archivo en la carpeta donde están trabajando con R.

dnaSeqs = readDNAStringSet("ecoliORFs.fa")
dnaSeqs

##   A DNAStringSet instance of length 11
##      width seq                                         names               
##  [1]    66 ATGAAACGCATTAGCACCAC...CAGGTAACGGTGCGGGCTGA thrL
##  [2]  2463 ATGCGAGTGTTGAAGTTCGG...CATGGAAGTTAGGAGTCTGA thrA
##  [3]   933 ATGGTTAAAGTTTATGCCCC...CACGAGTACTGGAAAACTAA thrB
##  [4]  1287 ATGAAACTCTACAATCTGAA...TGATGATGAATCATCAGTAA thrC
##  [5]   504 ATGCCGGGCAACAGCCCGCA...ACATAAAACACTATCAATAA insB
##  ...   ... ...
##  [7]   606 ATGGCAGAGAAATTTATCAA...AACCTGCGTTTATGAATTAA leuD
##  [8]  1401 ATGGCTAAGACGTTATACGA...ACATTCGCAACATTAAATAA leuC
##  [9]  1092 ATGTCGAAGAATTACCATAT...ATGTAGCAGAAGGGGTGTAA leuB
## [10]  1572 ATGAGCCAGCAAGTCATTAT...ACAACAAGGAAACCGTGTGA leuA
## [11]    87 ATGACTCACATCGTTCGCTT...TGAGCGGCATCCAGCATTAA leuL

Un DNAStringSet es simplement una lista de objetos DNAString. Algunas funciones útiles para manipularlas son las siguientes (consulten la ayuda si quieren saber más):

length(dnaSeqs)

## [1] 11

names(dnaSeqs)

##  [1] "thrL" "thrA" "thrB" "thrC" "insB" "insA" "leuD" "leuC" "leuB" "leuA"
## [11] "leuL"

width(dnaSeqs)

##  [1]   66 2463  933 1287  504  276  606 1401 1092 1572   87

subseq(dnaSeqs,11,30)

##   A DNAStringSet instance of length 11
##      width seq                                         names               
##  [1]    20 TTAGCACCACCATTACCACC                        thrL
##  [2]    20 TGAAGTTCGGCGGTACATCA                        thrA
##  [3]    20 TTTATGCCCCGGCTTCCAGT                        thrB
##  [4]    20 ACAATCTGAAAGATCACAAC                        thrC
##  [5]    20 ACAGCCCGCATTATGGGCGT                        insB
##  ...   ... ...
##  [7]    20 AATTTATCAAACACACAGGC                        leuD
##  [8]    20 CGTTATACGAAAAATTGTTC                        leuC
##  [9]    20 ATTACCATATTGCCGTATTG                        leuB
## [10]    20 AAGTCATTATTTTCGATACC                        leuA
## [11]    20 TCGTTCGCTTTATCGGTCTA                        leuL

dinucleotideFrequency(dnaSeqs)

##        AA  AC  AG  AT  CA  CC  CG  CT  GA  GC  GG  GT TA  TC  TG  TT
##  [1,]   3  10   2   5  11   7   3   1   2   4   4   2  4   1   3   3
##  [2,] 176 113 102 161 130 143 199 143 160 220 170 142 86 139 221 157
##  [3,]  52  40  52  50  52  50  80  49  59  94  85  56 31  47  77  58
##  [4,] 103  66  65  70  74  64 105  73  93 119  79  73 34  67 115  86
##  [5,]  36  31  21  33  34  18  45  29  29  50  43  27 22  27  40  18
##  [6,]  21  20   9  14  25  16  18  21   7  25  19  15 12  19  19  15
##  [7,]  50  37  29  33  36  34  45  35  44  55  37  26 19  24  51  50
##  [8,] 131  86  54  74  90 101 133  62  80 133 110  72 44  66  98  66
##  [9,]  69  60  50  73  75  73 100  50  70 119  67  52 38  46  91  58
## [10,] 138  91  78 103 111  97 138  65 100 136  97  89 61  87 109  71
## [11,]   3   7   4   6   6   1   8   8   4   7   3   4  7   8   3   7

reverseComplement(dnaSeqs)

##   A DNAStringSet instance of length 11
##      width seq                                         names               
##  [1]    66 TCAGCCCGCACCGTTACCTG...GTGGTGCTAATGCGTTTCAT thrL
##  [2]  2463 TCAGACTCCTAACTTCCATG...CCGAACTTCAACACTCGCAT thrA
##  [3]   933 TTAGTTTTCCAGTACTCGTG...GGGGCATAAACTTTAACCAT thrB
##  [4]  1287 TTACTGATGATTCATCATCA...TTCAGATTGTAGAGTTTCAT thrC
##  [5]   504 TTATTGATAGTGTTTTATGT...TGCGGGCTGTTGCCCGGCAT insB
##  ...   ... ...
##  [7]   606 TTAATTCATAAACGCAGGTT...TTGATAAATTTCTCTGCCAT leuD
##  [8]  1401 TTATTTAATGTTGCGAATGT...TCGTATAACGTCTTAGCCAT leuC
##  [9]  1092 TTACACCCCTTCTGCTACAT...ATATGGTAATTCTTCGACAT leuB
## [10]  1572 TCACACGGTTTCCTTGTTGT...ATAATGACTTGCTGGCTCAT leuA
## [11]    87 TTAATGCTGGATGCCGCTCA...AAGCGAACGATGTGAGTCAT leuL

translate(dnaSeqs)

##   A AAStringSet instance of length 11
##      width seq                                         names               
##  [1]    22 MKRISTTITTTITITTGNGAG*                      thrL
##  [2]   821 MRVLKFGGTSVANAERFLRV...TAAGVFADLLRTLSWKLGV* thrA
##  [3]   311 MVKVYAPASSANMSVGFDVL...EGFVHICRLDTAGARVLEN* thrB
##  [4]   429 MKLYNLKDHNEQVSFAQAVT...SHNLPADFAALRKLMMNHQ* thrC
##  [5]   168 MPGNSPHYGRWPQHDFTSLK...SVELHDKVIGHYLNIKHYQ* insB
##  ...   ... ...
##  [7]   202 MAEKFIKHTGLVVPLDAANV...LQHDDAIAAYEAKQPAFMN* leuD
##  [8]   467 MAKTLYEKLFDAHVVYEAEN...AMAAAAAVTGHFADIRNIK* leuC
##  [9]   364 MSKNYHIAVLPGDGIGPEVM...AVSTDEMGDIIARYVAEGV* leuB
## [10]   524 MSQQVIIFDTTLRDGEQALQ...VEKELQRKAQHNENNKETV* leuA
## [11]    29 MTHIVRFIGLLLLNASSLRGRRVSGIQH*               leuL

Y al final siempre podemos accesar a la secuencia como una cadena de texto normal:

as.character(dnaSeqs[[1]][1:20])

## [1] "ATGAAACGCATTAGCACCAC"

Por último, vale la pena saber que casi todos los paquetes de Bioconductor tienen manuales o vignettes que los acompañan. Para accesar la lista de manuales disponibles, usen la función:

openVignette()

Eligiendo el número del manual que quieren ver, en principio lo debe abrir como PDF. Por ejemplo, pueden buscar el que se llama “Biostrings - Biostrings Quick Overview” para ver una descripción de la mayoría de las funciones nuevas que están incluidas en Biostrings.

Ejercicios:

Nota: todos estos ejercicios se pueden solucionar exclusivamente con las funciones vistas dentro de las prácticas.

Cuántas secuencias miden más de 1000?
Cuál es el %GC de la secuencia más corta? (tip: ?alphabetFrequency)
Obtén el primer codón (los primeros tres nucleótidos) de cada secuencia. Para qué aminoácido codifican? (tip: prueba la función table para agrupar o tabular tus resultados)
Ahora obtén el último codón (los últimos tres nucleótidos) de cada secuencia. Intenta hacerlo de más de una forma. Cómo puedes demostrar que tu resultado está bien?
Cuál es la secuencia que tiene más dinucleótidos GC? y si tomas en cuenta la longitud (esto es, por frecuencia en lugar de por número total)?
Muestra la frecuencia de nucleótidos de todas las secuencias gráficamente (barplot)