Bioconductor

Bioconductor es un proyecto de código abierto con una colección de paquetes de R para la recolección y el análisis de datos genómicos de alto rendimiento. Tiene herramientas y paquetes para microarrays, preprocesamiento, diseño experimental, enfoques integradores y reproducibles de la bioinformática tareas, etc. El sitio web de Bioconductor (http://www.bioconductor.org) proporciona detalles sobre instalación, repositorio de paquetes, ayuda y otra documentación.

Este tutorial pretende introducir a los estudiantes con los datos mas comunes en bioinformática con los paquetes disponibles en Bioconductor.

Instalación de paquetes de Bioconductor

Bioconductor, como se mencionó anteriormente, tiene una colección de paquetes y software que sirven diferentes propósitos en el análisis de datos biológicos u otros objetivos en bioinformática. En orden para poder usarlos, debe instalarlos y cargarlos en su entorno R o sesión R.

Para instalar paquetes desde Bioconductor se deben seguir los siguientes pasos:

  1. Utilice el script biocLite.R para instalar paquetes de bioconductor:
# instalar instalador
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")

# instalar bioconductor
BiocManager::install(version = "3.11")
  1. Instalar paquetes específicos, por ejemplo, GenomicFeatures y Biostrings:
BiocManager::install(c("GenomicFeatures", "Biostrings"))
BiocManager::available(pattern = "genet")

Biostrings

Empezaremos aprendiendo a manipular cadenas de caracteres (character strings) usando el paquete Biostrings. Empecemos por instalar y cargar el paquete (si aun no lo han hecho):

BiocManager::install("Biostrings")

library(Biostrings)

R proporciona funciones que generan letras mayúsculas y minúsculas.

letters
##  [1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l" "m" "n" "o" "p" "q" "r" "s"
## [20] "t" "u" "v" "w" "x" "y" "z"
LETTERS
##  [1] "A" "B" "C" "D" "E" "F" "G" "H" "I" "J" "K" "L" "M" "N" "O" "P" "Q" "R" "S"
## [20] "T" "U" "V" "W" "X" "Y" "Z"

Una forma de generar las letras del código genético “ACGT” es hacer:

sample(LETTERS[c(1, 3, 7, 20)], size = 10, replace = TRUE)
##  [1] "T" "T" "A" "G" "C" "C" "A" "T" "G" "C"

Biostrings facilita esto un poco con la variable DNA_ALPHABET:

DNA_ALPHABET
##  [1] "A" "C" "G" "T" "M" "R" "W" "S" "Y" "K" "V" "H" "D" "B" "N" "-" "+" "."
seq <- sample(DNA_ALPHABET[1:4], size = 10, replace = TRUE)
seq
##  [1] "A" "G" "T" "G" "A" "C" "T" "T" "G" "A"

Esto lo podemos convertir en una cadena:

seq <- paste(seq, collapse = "")

seq
## [1] "AGTGACTTGA"

Las clases (X)String

La clase XString nos permite crear, almacenar y trabajar con diferentes tipos de cadenas. Solo podemos trabajar con subclases de XString: BString, DNAString, RNAString y AAString.

Primero creemos un objeto BString:

bstring <- BString("Soy un objeto BString")

bstring
## 21-letter BString object
## seq: Soy un objeto BString

Evaluemos el largo de la cadena con length():

length(bstring)
## [1] 21

Ahora creemos un objeto DNAString:

dnastring <- DNAString("TTGAAA-CTC-N")

dnastring
## 12-letter DNAString object
## seq: TTGAAA-CTC-N

Tambien podemos medir el largo con length():

length(dnastring)
## [1] 12

Las diferencias con la clase BString son:

  • Solo admite letras del alfabeto extendido genetico IUPAC y el caracter de huecos (-)
  • Cada letra es codificada antes de ser guardada

Podemos ver al estructura del objeto XString:

str(dnastring)
## Formal class 'DNAString' [package "Biostrings"] with 5 slots
##   ..@ shared         :Formal class 'SharedRaw' [package "XVector"] with 2 slots
##   .. .. ..@ xp                    :<externalptr> 
##   .. .. ..@ .link_to_cached_object:<environment: 0x55f4cd5cefa0> 
##   ..@ offset         : int 0
##   ..@ length         : int 12
##   ..@ elementMetadata: NULL
##   ..@ metadata       : list()
slotNames(dnastring)
## [1] "shared"          "offset"          "length"          "elementMetadata"
## [5] "metadata"

Hay algunas otras funciones útiles para manipular estos objetos:

alphabetFrequency(dnastring)
## A C G T M R W S Y K V H D B N - + . 
## 3 2 1 3 0 0 0 0 0 0 0 0 0 0 1 2 0 0
alphabetFrequency(dnastring, baseOnly = TRUE, as.prob = TRUE)
##        A        C        G        T    other 
## 0.250000 0.166667 0.083333 0.250000 0.250000
letterFrequency(dnastring, "A", as.prob = TRUE)
##    A 
## 0.25
reverseComplement(dnastring)
## 12-letter DNAString object
## seq: N-GAG-TTTCAA

Podemos hacer subconjuntos usando indexación:

bstring[3]
## 1-letter BString object
## seq: y
dnastring[2]
## 1-letter DNAString object
## seq: T
dnastring[7:12]
## 6-letter DNAString object
## seq: -CTC-N

Si bien este método para obtener subcadenas es fácil e intuitivo, puede ser lento con estos objetos, especialmente para cadenas grandes. Biostrings proporciona una función subseq() que aprovecha los componentes internos de XString.

subseq(bstring, start = 3, end = 3)
## 1-letter BString object
## seq: y
subseq(bstring, 3, 3)
## 1-letter BString object
## seq: y
subseq(dnastring, 7, 12)
## 6-letter DNAString object
## seq: -CTC-N

Podemos convertir cualquier objeto XString en un vector de caracteres con toString()

toString(bstring)
## [1] "Soy un objeto BString"
toString(dnastring)
## [1] "TTGAAA-CTC-N"

Comparing XStrings

Podemos comparar objetos XStrings con cadenas de caracteres y con objetos XStrings:

bb <- subseq(bstring, 3, 6)

bb == "am a"
## [1] FALSE
bb == BString("am a")
## [1] FALSE

Traten de no comparar diferentes clases de XString entre sí.

bstring == dnastring
## Error in bstring == dnastring: comparison between a "BString" instance and a "DNAString" instance is not supported
bstring != dnastring
## Error in e1 == e2: comparison between a "BString" instance and a "DNAString" instance is not supported

Las clases (X)StringViews

Esta es una forma útil de ver múltiples subsecuencias de un XString.

view <- Views(dnastring, start = 3:0, end = 5:8)

view
## Views on a 12-letter DNAString subject
## subject: TTGAAA-CTC-N
## views:
##       start end width
##   [1]     3   5     3 [GAA]
##   [2]     2   6     5 [TGAAA]
##   [3]     1   7     7 [TTGAAA-]
##   [4]     0   8     9 [ TTGAAA-C]

Para revisar el largo:

length(view)
## [1] 4

Y nuevamente podemos ver subconjuntos usando índices:

view[4:2]
## Views on a 12-letter DNAString subject
## subject: TTGAAA-CTC-N
## views:
##       start end width
##   [1]     0   8     9 [ TTGAAA-C]
##   [2]     1   7     7 [TTGAAA-]
##   [3]     2   6     5 [TGAAA]

El objeto devuelto sigue siendo un objeto XStringViews, incluso si seleccionamos solo un elemento. Para extraer cualquier vista dada como un objeto XString, use el operador [[]].

view[[2]]
## 5-letter DNAString object
## seq: TGAAA

Para ver todas las vistas menos la tercera, use:

view[-3]
## Views on a 12-letter DNAString subject
## subject: TTGAAA-CTC-N
## views:
##       start end width
##   [1]     3   5     3 [GAA]
##   [2]     2   6     5 [TGAAA]
##   [3]     0   8     9 [ TTGAAA-C]

XStringSet

Los objetos XStringSet proporcionan una forma conveniente de almacenar varios objetos XString juntos:

set <- NULL

for (i in 1:4) set <- c(set, paste(sample(DNA_ALPHABET[1:4], 10, replace = T), collapse = ""))

set
## [1] "TCGCCCTACT" "GTAGGCTCTC" "CGTGGTAGTA" "CATTAACGCG"
length(set)
## [1] 4

Ahora convirtámoslo en un conjunto de DNAStrings:

set <- DNAStringSet(set)

set
## DNAStringSet object of length 4:
##     width seq
## [1]    10 TCGCCCTACT
## [2]    10 GTAGGCTCTC
## [3]    10 CGTGGTAGTA
## [4]    10 CATTAACGCG
length(set)
## [1] 4

Las funciones reverseComplement(), alphabetFrequency(), width() y subseq() funcionan en cada una de las XStrings en un XStringSet:

reverseComplement(set)
## DNAStringSet object of length 4:
##     width seq
## [1]    10 AGTAGGGCGA
## [2]    10 GAGAGCCTAC
## [3]    10 TACTACCACG
## [4]    10 CGCGTTAATG
alphabetFrequency(set, baseOnly = T, as.prob = T)
##        A   C   G   T other
## [1,] 0.1 0.5 0.1 0.3     0
## [2,] 0.1 0.3 0.3 0.3     0
## [3,] 0.2 0.1 0.4 0.3     0
## [4,] 0.3 0.3 0.2 0.2     0
letterFrequency(set, "A", as.prob = T)
##        A
## [1,] 0.1
## [2,] 0.1
## [3,] 0.2
## [4,] 0.3
width(set)
## [1] 10 10 10 10
subseq(set, 1, 3)
## DNAStringSet object of length 4:
##     width seq
## [1]     3 TCG
## [2]     3 GTA
## [3]     3 CGT
## [4]     3 CAT

Se puede nombrar cada XString:

names(set) = paste("name", 1:4, sep = "")

Ejercicio 1


Examinaremos el genoma de Escherichia coli APEC O1 (NC 008563). Esto se puede encontrar en el paquete BSgenome (junto con genomas de otros organismos modelo). Descargue e instale los datos del paquete:

BiocManager::install("BSgenome.Ecoli.NCBI.20080805")

Ahora podemos cargar el genoma con:

library(BSgenome.Ecoli.NCBI.20080805)

eco <- Ecoli$NC_008563


1.1 Muestree el genoma generando 1000 vistas (Views) aleatorias con largos aleatorios entre 50 y 100 bases.


1.2 Utilice las funciones alphabetFrequency() y reverseComplement() en estas vistas


1.3 Repita 1.1 y 1.2 pero esta vez con solo 100 muestras


1.4 Compare los resultados de alphabetFrequency() ¿Cómo se comparan ambos con el genoma completo?

 

IRanges

Al igual que las vistas (Views), la función IRanges() se puede utilizar para ver sub-secuencias. Por ejemplo, una tarea común es especificar una lista de posiciones iniciales en un cromosoma y luego mirar sub-secuencias de una longitud determinada con ese punto de inicio:

ir1 <- IRanges(start = 1:10, width = 10:1)

ir2 <- IRanges(start = 1:10, end = 11)

ir3 <- IRanges(end = 11, width = 10:1)

Exploremos la estructura del resultado de la función:

str(ir1)
## Formal class 'IRanges' [package "IRanges"] with 6 slots
##   ..@ start          : int [1:10] 1 2 3 4 5 6 7 8 9 10
##   ..@ width          : int [1:10] 10 9 8 7 6 5 4 3 2 1
##   ..@ NAMES          : NULL
##   ..@ elementType    : chr "ANY"
##   ..@ elementMetadata: NULL
##   ..@ metadata       : list()

Podemos extraer valores del inicio, fin y ancho con estas funciones:

start(ir1)
##  [1]  1  2  3  4  5  6  7  8  9 10
end(ir1)
##  [1] 10 10 10 10 10 10 10 10 10 10
width(ir1)
##  [1] 10  9  8  7  6  5  4  3  2  1

Asi como con las vistas (Views), podemos hacer subconjuntos de los objetos IRanges para obtener otro objecto IRanges:

ir1[1:4]
## IRanges object with 4 ranges and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]         1        10        10
##   [2]         2        10         9
##   [3]         3        10         8
##   [4]         4        10         7

Esto tambien se puede hacer utilizando expresiones lógicas:

ir1[start(ir1) <= 6]
## IRanges object with 6 ranges and 0 metadata columns:
##           start       end     width
##       <integer> <integer> <integer>
##   [1]         1        10        10
##   [2]         2        10         9
##   [3]         3        10         8
##   [4]         4        10         7
##   [5]         5        10         6
##   [6]         6        10         5
sep <- 0.5
height <- 1
xlim <- c(min(start(ir1)) - sep, max(end(ir1)) + sep)
bins <- disjointBins(IRanges(start(x), end(x) + 1))

plot.new()
plot.window(xlim, c(0, max(bins) * (height + sep)))
ybottom = bins * (sep + height) - height
rect(start(x) - 0.5, ybottom, end(x) + 0.5, ybottom + height, col = "#35B779FF")
axis(1)

Tambien podemos hacerlo usando el paquete para graficación ggplot2:

sep <- 0.5
height <- 1
if (is(xlim, "Ranges")) xlim = c(min(start(xlim)) - sep, max(end(xlim)) + sep)
bins = disjointBins(IRanges(start(x), end(x) + 1))
ybottom = bins * (sep + height) - height
df = data.frame(ybottom = ybottom, xleft = start(x) - 0.5, xright = end(x) + 0.5, 
    ytop = ybottom + height)
ggplot(df) + geom_rect(aes(xmax = xright, xmin = xleft, ymax = ytop, ymin = ybottom))

Ejercicio 2


2.1 Cree una función que genere el gráfico de objetos IRanges utilizando R básico como en el código tras-anterior


2.2 Haga que la función pueda crear el gráfico con R básico o con ggplot2. Para esto agregue un argumento que permita al usuario decidir que opción usar


 

Islas CpG

El contenido de GC es el porcentaje del genoma (o fragmento de ADN) que es “G” o “C”. Para calcular el contenido de GC, contamos las apariciones de los alfabetos “G” y “C” y lo dividimos por la longitud de la cadena en cuestión. Las islas CpG regiones donde existe una gran concentración de pares de citosina y guanina enlazados por fosfatos. La “p” en CpG representa que están enlazados por un fosfato. La definición formal de una isla CpG es una región de tamaño igual o superior a 500 pb y con un porcentaje de GC mayor al 50% y con un razón de CpG observado a esperado mayor a 0.6. Gracias a estas regiones podemos identificar dos tipos de genes: genes constitutivos (70% de los genes del genoma), cuyo promotor tiene un promedio de CpG en torno al %61 y genes que se expresan en tejidos específicos (30% de los genes) que tienen un promedio de CpG en torno al 23%. Las islas CpG marcan regiones importantes en el genoma porque más del 65% de las regiones promotoras de genes se pueden encontrar dentro de estas “islas”.

Para explorar esta caraterística de las secuencias genéticas usaremos datos del cromosoma 8 de la versión 19 del genoma humano. Podemos obtener los datos leyéndolos directamente de mediación virtual:

seqChr8 <- readDNAStringSet("seqChr8.fasta")[[1]]

A continuación, vamos a descargar las ubicaciones de las islas CpG en esta versión del genoma de esta dirección web: http://rafalab.jhsph.edu/CGI/model-based-cpg-islands-hg19.txt. El archivo ‘model-based-cpg-islands-hg19.txt’ también se encuentra en mediación virtual. Guárdelo en su directorio de trabajo y léalo con los siguientes comandos:

cpglocs <- read.table("model-based-cpg-islands-hg19.txt", header = TRUE)

Podemos explorar el archivo así:

head(cpglocs)
chr start end length CpGcount GCcontent pctGC obsExp
chr10 93098 93818 721 32 403 0.559 0.572
chr10 94002 94165 164 12 97 0.591 0.841
chr10 94527 95302 776 65 538 0.693 0.702
chr10 119652 120193 542 53 369 0.681 0.866
chr10 122133 122621 489 51 339 0.693 0.880
chr10 180265 180720 456 32 256 0.561 0.893

Ahora debemos extraer el inicio y fin de las islas CpG para el cromosoma 8:

cpglocs8 <- cpglocs[which(cpglocs[, 1] == "chr8"), ]

# revisar dimiensiones
dim(cpglocs8)
## [1] 2855    8
# ver las primeras filas
head(cpglocs8)
chr start end length CpGcount GCcontent pctGC obsExp
56961 chr8 32437 33708 1272 116 795 0.625 0.934
56962 chr8 40354 40656 303 29 180 0.594 1.093
56963 chr8 44536 46203 1668 96 1045 0.626 0.618
56964 chr8 99457 100721 1265 108 958 0.757 0.595
56965 chr8 155954 156761 808 85 549 0.679 0.912
56966 chr8 179033 179756 724 50 424 0.586 0.824

 

La columna ‘length’ indica el largo de la secuencia, la columna ‘GCcontent’ el número de pares de bases GC y la columna ‘obsExp’ .

Ejercicio 3


3.1 Calcule el porcentaje de pares de bases GC para la secuencia completa del chromosoma 8


3.2 Extraiga las secuencias de islas CpG en un objeto “DNAStringSet”


3.3 Extraiga las secuencias que no corresponden a islas CpG en otro objeto “DNAStringSet”


3.4 Calcule el porcentaje de bases CG para ambos subconjuntos de datos

 

Información de la sesión

## R version 4.0.1 (2020-06-06)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] BSgenome.Ecoli.NCBI.20080805_1.3.1000 BSgenome_1.58.0                      
##  [3] rtracklayer_1.50.0                    GenomicFeatures_1.42.0               
##  [5] AnnotationDbi_1.52.0                  Biobase_2.50.0                       
##  [7] GenomicRanges_1.42.0                  GenomeInfoDb_1.26.0                  
##  [9] Biostrings_2.58.0                     XVector_0.30.0                       
## [11] IRanges_2.24.0                        S4Vectors_0.28.0                     
## [13] BiocGenerics_0.36.0                   BiocManager_1.30.10                  
## [15] kableExtra_1.1.0                      knitr_1.29                           
## [17] ggplot2_3.3.2                         RColorBrewer_1.1-2                   
## 
## loaded via a namespace (and not attached):
##  [1] MatrixGenerics_1.2.0        httr_1.4.1                 
##  [3] bit64_4.0.5                 viridisLite_0.3.0          
##  [5] assertthat_0.2.1            askpass_1.1                
##  [7] highr_0.8                   BiocFileCache_1.14.0       
##  [9] blob_1.2.1                  GenomeInfoDbData_1.2.4     
## [11] Rsamtools_2.6.0             yaml_2.2.1                 
## [13] progress_1.2.2              lattice_0.20-41            
## [15] pillar_1.4.4                RSQLite_2.2.1              
## [17] glue_1.4.1                  digest_0.6.25              
## [19] rvest_0.3.5                 colorspace_1.4-1           
## [21] Matrix_1.2-18               htmltools_0.5.0            
## [23] XML_3.99-0.5                pkgconfig_2.0.3            
## [25] biomaRt_2.46.0              zlibbioc_1.36.0            
## [27] purrr_0.3.4                 scales_1.1.1               
## [29] webshot_0.5.2               BiocParallel_1.24.0        
## [31] tibble_3.0.1                openssl_1.4.1              
## [33] farver_2.0.3                generics_0.0.2             
## [35] ellipsis_0.3.1              SummarizedExperiment_1.20.0
## [37] withr_2.2.0                 magrittr_1.5               
## [39] crayon_1.3.4                memoise_1.1.0              
## [41] evaluate_0.14               xml2_1.3.2                 
## [43] tools_4.0.1                 prettyunits_1.1.1          
## [45] hms_0.5.3                   formatR_1.7                
## [47] matrixStats_0.57.0          lifecycle_0.2.0            
## [49] stringr_1.4.0               munsell_0.5.0              
## [51] DelayedArray_0.16.0         compiler_4.0.1             
## [53] rlang_0.4.8                 grid_4.0.1                 
## [55] RCurl_1.98-1.2              rstudioapi_0.11            
## [57] rappdirs_0.3.1              labeling_0.3               
## [59] bitops_1.0-6                rmarkdown_2.3              
## [61] gtable_0.3.0                DBI_1.1.0                  
## [63] curl_4.3                    R6_2.4.1                   
## [65] GenomicAlignments_1.26.0    dplyr_1.0.2                
## [67] bit_4.0.4                   readr_1.3.1                
## [69] stringi_1.4.6               Rcpp_1.0.4.6               
## [71] vctrs_0.3.4                 dbplyr_1.4.4               
## [73] tidyselect_1.1.0            xfun_0.15