Ejercicio 1

a) ¿Cuáles son las principales diferencias entre el aprendizaje supervisado y no supervisado?

El aprendizaje supervisado está destinado a la realización de predicciones, mientras que el aprendizaje supervisado busca descrubrir patrones en un conjunto de datos.

El aprendizaje supervisado utiliza conjuntos de datos etiquetados, lo que significa que cada entrada posee una etiqueta o resultado asociado. El aprendizaje no supervisado, sin embargo, se aplica a conjuntos de datos no etiquetados o no procesados.

El aprendizaje supervisado, a su vez, utiliza un conjunto de datos de entrenamiento y otro de prueba. Se realiza un entrenamiento de algoritmos para predecir los resultados en el conjunto de prueba y se compara con las etiquetas reales de dicho conjunto.

b) Cita un ejemplo de caso real del ámbito biosanitario, en el que sea necesario el uso de conceptos del aprendizaje supervisado.

El aprendizaje supervisado se podría aplicar en el ámbito sanitario a conjuntos de datos con biomarcadores para la inferencia del riesgo de infección urinaria. Estos datos estarían etiquetados según el paciente desarrollara una infección posteriormente o no. Esto contribuiría a evitar diagnosis y administraciones de antibióticos prematuras. Esta aplicación ya ha sido implementada por Hernandez et al., 2017.

c) Cita un ejemplo de caso real del ámbito biosanitario, en el que sea necesario el uso de conceptos del aprendizaje no supervisado.

Se podría aplicar el aprendizaje no supervisado a conjuntos de datos de expresión génica de pacientes con cáncer. Usando técnicas de agrupamiento, se podría distinguir qué subtipo de cáncer presenta un paciente, para realizar una estimación del tiempo de vida. Esta aplicación ya ha sido implementada por Von Heydebreck et al., 2001.

Ejercicio 2

a) Realizar un agrupamiento jerárquico aglomerativo, especificando los valores de distancia y dibujando el dendograma asociado. Comentad los resultados.

# Estandarizamos los datos
swiss <- scale(swiss)

# Calculamos la matriz de distancias
dist_matrix <- dist(swiss, method = "euclidean")

# Calculamos número recomendado de clústeres (como orientación)
# Creamos un vector vacío para incorporar las sumas de cuadrados
set.seed(999)
sc_cluster <- vector()

# Generamos suma de las varianzas para cada número de clústeres y representamos (método del codo)
for (i in 1:10) sc_cluster[i] <- sum(kmeans(swiss, i)$withinss)
plot(1:10,
     sc_cluster,
     type = 'b', 
     main = paste('Método del codo'),
     xlab = 'Número de clústeres',
     ylab = 'Suma de cuadrados')

Representando la suma de las varianzas en el clúster a medida que aumentamos el número de clústeres, podemos observar que a partir de 3 clústeres la reducción de la varianza disminuye significativamente. Se obtienen los mismos resultados con el método de la silueta (silhouette).

# Calculamos grupos con hclust utilizando distintos métodos
hcw <- as.dendrogram(hclust(dist_matrix, method = "ward.D2"))
hcc <- as.dendrogram(hclust(dist_matrix, method = "complete"))
hca <- as.dendrogram(hclust(dist_matrix, method = "average"))

# Representamos los dendrogramas
par(mfrow=c(3,1), mar=c(2, 3, 2, 1))

plot(hcw, main = "Dendrograma de clúster Ward.D2", leaflab = "none")
title(xlab = "Región", line = 0, cex.lab = 1.2)
title(ylab = "Altura", line = 2, cex.lab = 1.2)
abline(h = 8, col="red")
abline(h = 12, col="blue")
legend(x = "topright", legend = c("k=2", "k=3"), col = c("blue", "red"), lwd = 2) 

plot(hcc, main = "Dendrograma de clúster Complete", leaflab = "none")
title(xlab = "Región", line = 0, cex.lab = 1.2)
title(ylab = "Altura", line = 2, cex.lab = 1.2)
abline(h = 5.3, col="red")
abline(h = 7, col="blue")
legend(x = "topright", legend = c("k=2", "k=3"), col = c("blue", "red"), lwd = 2) 

plot(hca, main = "Dendrograma de clúster Average", leaflab = "none")
title(xlab = "Región", line = 0, cex.lab = 1.2)
title(ylab = "Altura", line = 2, cex.lab = 1.2)
abline(h = 4, col="red")
abline(h = 5, col="blue")
legend(x = "topleft", legend = c("k=2", "k=3"), col = c("blue", "red"), lwd = 2)

Podemos observar que los métodos Ward D2 y Complete producen clústeres más equilibrados, que el método Average. En dicho método, si escogiéramos un punto de corte para obtener 3 clústeres, 2 de ellos contendrían una región y el resto de regiones se agruparía en un único clúster. Por tanto, el método Average no es conveniente para realizar la agrupación.

En caso de considerar 3 clústeres, el método Ward D2 sería el más conveniente, ya que el método Complete produciría un clúster con un único elemento.

# Calculamos el dendrograma con método Ward D2
hcw <- hclust(dist_matrix, method = "ward.D2")
hcw_den <- as.dendrogram(hcw)

# Representamos el dendrograma
fviz_dend(hcw_den, k = 3, cex = 0.55, k_colors = c("#00AFBB", "#E7B800", "#FC4E07"), 
color_labels_by_k = TRUE, rect = TRUE, xlab = "Región", y = "Altura", 
main = "Dendrograma de clúster", horiz = TRUE)

Si observamos qué regiones se agrupan en los clústeres, encontramos un clúster pequeño (clúster amarillo) que engloba las regiones más importantes, como Lausanne, Neuchatel, Geneve y alrededores. Este clúster se caracteriza por tener una menor fertilidad y ocupación agrícola, y un nivel educativo y una evaluación en el servicio militar mayor que los otros clústeres. El resto de regiones, a su vez, puede dividirse en el clúster 1 y clúster 2, en los que se observa un alto y un bajo porcentaje de católicos, respectivamente.

# Elegimos el nivel de corte del dendrograma
nivel_corte <- 3

cluster_asignado <- cutree(hcw, k = nivel_corte)

# Guardamos resultados de agrupaciones en un dataframe  
(head (resultado <- data.frame(Región = 1:nrow(swiss), Clúster = cluster_asignado)))    
##              Región Clúster
## Courtelary        1       1
## Delemont          2       2
## Franches-Mnt      3       2
## Moutier           4       1
## Neuveville        5       1
## Porrentruy        6       2

b) Realizar un agrupamiento jerárquico aglomerativo, utilizando la función agnes(), dibujar el dendograma y calcular el coeficiente de aglomeración. Comentad los resultados.

# Calculamos coeficiente de aglomeración para distintos métodos de agnes
vector <- c("ward", "complete", "average")

for (i in vector){
  print(i)
  print(agnes(swiss, method = i)$ac)
}
## [1] "ward"
## [1] 0.9148672
## [1] "complete"
## [1] 0.864543
## [1] "average"
## [1] 0.7958782

El coeficiente de aglomeración nos indica la representación de las distancias iniciales, calculadas por la matriz de distancias, en el cálculo de las distancias finales, representadas en el dendrograma. El método con el coeficiente de aglomeración más alto es el método Ward. El coeficiente es 0.91, muy cercano a 1, por lo que la agrupación es fuerte. Por tanto, continuaremos con dicho método para calcular el dendrograma.

# Realizamos agrupamiento con agnes con método Ward
agw_den <- as.dendrogram(agnes(swiss, method="ward")) 

# Representamos el dendrograma
fviz_dend(agw_den, k = 3, cex = 0.55, k_colors = c("#FC4E07", "#E7B800", "#00AFBB"), 
color_labels_by_k = TRUE, rect = TRUE, xlab = "Región", y = "Altura", 
main = "Dendrograma de agnes", horiz = TRUE)

# Calculamos correlación entre dendrograma de clúster y agnes
lista_dend <- dendlist(hcw_den, agw_den)

cor.dendlist(lista_dend)
##      [,1] [,2]
## [1,]    1    1
## [2,]    1    1

Observamos que las regiones agrupadas con la función de clúster y agnes son las mismas.

c) Realizar un agrupamiento divisional jerárquico, utilizar la función diana(), dibujar el dendongrama y calcular el coeficiente de división. Comentad los resultados.

# Realizamos agrupamiento divisional jerárquico
hc_div <- diana(swiss)

# Representamos el dendrograma
fviz_dend(hc_div, cex = 0.6, xlab = "Región", k = 3, k_colors = c("#FC4E07", "#E7B800", "#00AFBB", "#2E4FDF"), 
color_labels_by_k = TRUE, rect = TRUE, y = "Altura", 
main = "Dendrograma de Diana", horiz = TRUE)

Podemos observar que el agrupamiento divisional jerárquico produce resultados distintos a los obtenidos con agnes y hclust, en el caso de establecer 3 clústeres. Este agrupamiento establece un clúster con una sola región, V. De Geneve, otro clúster que incluye otras regiones importantes (Neuchatel, Laussane, alrededores de Geneve), y un clúster con el resto de regiones.

Sería conveniente establecer 4 clústeres en este caso, para tener en cuenta las diferencias entre las regiones con menor importancia socioeconómica, separadas según el porcentaje de católicos. Podemos observar los 4 clústeres en la siguiente representación.

# Asociamos los agrupamientos a los datos y representamos
grupos <- cutree(hc_div, k = 4)
fviz_cluster(list(data = swiss, cluster = grupos), main = "Representación de clústeres")

# Calculamos coeficiente de división
hc_div$dc
## [1] 0.862368

El coeficiente de división es relativamente próximo a 1, por lo que la agrupación de los datos también adecuada, aunque peor que la realizada con el agrupamiento aglomerativo.

Ejercicio 3

Responde a la cuestión y proporciona un ejemplo.

a) ¿Para qué sirve ANOVA y qué información nos proporciona?

El ANOVA es una técnica de análisis de varianza que nos indica si existen diferencias significativas (no debidas al azar) entre las medias de dos o más grupos. La técnica consiste en la comparación de la varianza entre los grupos con la varianza dentro de los grupos. Comparando dichos valores, podemos determinar si las diferencias entre las medias son estadísticamente significativas y, por tanto, determinar si el factor o factores que estamos estudiando tienen un efecto en nuestra población de estudio.

b) A nivel de ejemplo: elige un conjunto de datos del paquete MASS y especifica qué variables considerarás para realizar un ANOVA. Justifica tu elección y haz el ANOVA.

Para ejemplificar cómo realizaríamos un análisis de ANOVA, podemos utilizar el dataset anorexia. Podemos analizar si los distintos tratamientos (variable categórica) tienen un efecto significativo en la diferencia de peso (variable numérica), con respecto al control.

# Convertimos a factor y calculamos diferencia de pesos
tratamiento <- factor(anorexia$Treat)
diff_peso <- anorexia$Postwt - anorexia$Prewt

# Creamos dataframe con diferencia de pesos
datos_anorexia <- data.frame(Tratamiento = tratamiento, Diferencia = diff_peso)

# Comprobamos normalidad con test de Kolmogorov-Smirnov, ya que n > 50
by(data = datos_anorexia, INDICES = tratamiento, FUN = function(x){lillie.test(diff_peso)})
## tratamiento: CBT
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  diff_peso
## D = 0.096197, p-value = 0.0959
## 
## ------------------------------------------------------------ 
## tratamiento: Cont
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  diff_peso
## D = 0.096197, p-value = 0.0959
## 
## ------------------------------------------------------------ 
## tratamiento: FT
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  diff_peso
## D = 0.096197, p-value = 0.0959
# Comprobamos homocedasticidad
fligner.test(diff_peso ~ tratamiento)   
## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  diff_peso by tratamiento
## Fligner-Killeen:med chi-squared = 0.59304, df = 2, p-value = 0.7434

Aceptamos la hipótesis de normalidad, ya que el valor p es mayor al 0.05. Como los valores están muy cercanos al nivel de significación, realizamos la comprobación de la homocedasticidad mediante el test de Fligner-Killeen. El resultado indica que no se detecta falta de homocedasticidad en los datos. Por tanto, se podría proceder con el ANOVA.

# Media de diferencias de peso según tratamiento
aggregate(diff_peso ~ tratamiento, data = datos_anorexia, FUN = mean)   
##   tratamiento diff_peso
## 1         CBT  3.006897
## 2        Cont -0.450000
## 3          FT  7.264706
# Realizamos el ANOVA
resultado_anova <- aov(diff_peso ~ tratamiento, data = datos_anorexia)

# Resumen del ANOVA
summary(resultado_anova)
##             Df Sum Sq Mean Sq F value Pr(>F)   
## tratamiento  2    615  307.32   5.422 0.0065 **
## Residuals   69   3911   56.68                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El valor p es menor que un nivel de significación de 0.05, por lo que podemos concluir que existe una diferencia significativa en la diferencia de peso entre al menos dos de los tratamientos. De acuerdo al test de Tukey, únicamente observamos diferencias significativas entre el tratamiento Fam Tr y el control.

# Identificamos qué tratamientos difieren
TukeyHSD(resultado_anova)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = diff_peso ~ tratamiento, data = datos_anorexia)
## 
## $tratamiento
##               diff       lwr       upr     p adj
## Cont-CBT -3.456897 -8.327276  1.413483 0.2124428
## FT-CBT    4.257809 -1.250554  9.766173 0.1607461
## FT-Cont   7.714706  2.090124 13.339288 0.0045127

Ejercicio 4

Haced una aplicación con Shiny que permita introducir una cadena ADN y muestre el tamaño de esta cadena. Podeis seleccionar vosotros mismos la cadena ADN a introducir, por ejemplo:

cadADN1: ATCGATCGATCGATCGATCGATCGATCGATCG

A partir de la cadena definida y de la medida, comentad brevemente el resultado.

# Definimos la interfaz de usuario con fluidPage()
ui <- fluidPage(
  
  # Título de la aplicación
  titlePanel("Contador de nucleótidos"),
  
  # Entrada de texto con textInput()
  sidebarLayout(
    sidebarPanel(
      textInput("secuencia", "Cadena de ADN"),
      actionButton("contador", "Contar número")
    ),
  
  # Resultados panel principal
  mainPanel(
    verbatimTextOutput("resultado")
  )
 )
)

# Definimos los parámetros del servidor 
server <- function(input, output) { 
  
  # Definimos el botón   
  observeEvent(input$contador, {  
    
    # Eliminamos espacios de la string y contamos número de carácteres    
    ADN_sec <- gsub(" ", "", input$secuencia) 
    num_char <- nchar(ADN_sec)

    # Mostramos el resultado final     
    output$resultado <- renderText({       
      paste("El número de nucleótidos es:", num_char)     
      
    })   
  }) 
} 
  
# Ejecutamos la aplicación Shiny 
shinyApp(ui, server) 
## 
## Listening on http://127.0.0.1:8219

Cuando introducimos nuestra secuencia, cadADN1, el resultado es 32 nucleótidos.

Ejercicio 5

En el LAB6: Paquetes de R para la bioinformática, se ha trabajado con Bioconductor (https://www.bioconductor.org/) y varios paquetes (BiocManager, …), en este ejercicio se pide especificar y explicar el código en R de una aplicación práctica que utilice uno de estos paquetes, así como justificar la utilidad y objetivo de este código. Si el código mostrado no es de realización propia, es imprescindible citar la fuente.

En este ejercicio, comentamos una propuesta de uso de Bioconductor para analizar una muestra de Chip-seq, con la que se busca mapear los sitios de unión de factores de transcripción en un genoma de referencia. Realizamos el pre-procesamiento de las lecturas y el alineamiento con el genoma de ratón (mm10). El código ha sido realizado por la Universidad de Rockefeller (referencia disponible en la bibliografía).

file_dir <- "C:/Users/adrii/Downloads/PEC3/"

# Tomamos una muestra aleatoria para analizar la calidad de los datos
fqSample <- FastqSampler("chipseqexample.fastq", n = 10^6)

fastq <- yield(fqSample)

# Calculamos para cada lectura la suma de las probabilidades de las bases
readQuality <- quality(fastq)
readQualities <- alphabetScore(readQuality)
readQualities[1:10]
##  [1] 1103 1014   72  387  239   72  762   72 1081 1155
# Distribución de los scores de calidad 
toPlot <- data.frame(ReadQ = readQualities)
ggplot(toPlot, aes(x = ReadQ)) + geom_histogram() + theme_minimal()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

En el histograma, podemos observar que una parte de las lecturas presenta una calidad baja. Por tanto, es conveniente filtrar dichas lecturas.

# Calculamos frecuencia de cada nucleótido para cada lectura
readSequences <- sread(fastq)
readSequences_AlpFreq <- alphabetFrequency(readSequences)
readSequences_AlpFreq[1:3, ]
##       A  C  G  T M R W S Y K V H D B N - + .
## [1,]  8 12  8  8 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 12  2 10 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [3,]  8  7 18  3 0 0 0 0 0 0 0 0 0 0 0 0 0 0
# Sumamos frecuencias de cada nucleótido para todas las lecturas
summed__AlpFreq <- colSums(readSequences_AlpFreq)
summed__AlpFreq[c("A", "C", "G", "T", "N")]
##        A        C        G        T        N 
## 10259243  7714271  7439041 10442637   144808
# Comparamos frecuencia de bases para cada ciclo
readSequences_AlpbyCycle <- alphabetByCycle(readSequences)

AFreq <- readSequences_AlpbyCycle["A", ]
CFreq <- readSequences_AlpbyCycle["C", ]
GFreq <- readSequences_AlpbyCycle["G", ]
TFreq <- readSequences_AlpbyCycle["T", ]
toPlot <- data.frame(Count = c(AFreq, CFreq, GFreq, TFreq), Cycle = rep(1:36, 4),
    Base = rep(c("A", "C", "G", "T"), each = 36))

ggplot(toPlot, aes(y = Count, x = Cycle, colour = Base)) + geom_line() + theme_bw()

Podemos observar que las frecuencias de cada base se mantienen estables en los distintos ciclos.

# Creamos un archivo FASTQ, filtramos lecturas con calidad < 300 y contenido N > 10
fqStreamer <- FastqStreamer("chipseqexample.fastq", n = 1e+05)

TotalReads <- 0
TotalReadsFilt <- 0
while (length(fq <- yield(fqStreamer)) > 0) {
    TotalReads <- TotalReads + length(fq)
    filt1 <- fq[alphabetScore(fq) > 300]
    filt2 <- filt1[alphabetFrequency(sread(filt1))[, "N"] < 10]
    TotalReadsFilt <- TotalReadsFilt + length(filt2)
    writeFastq(filt2, paste0(file_dir, "filtered_chipseq.fastq.gz"), mode = "a")
}

# Descargamos manualmente el archivo con el genoma mm10 y lo instalamos
install.packages("BSgenome.Mmusculus.UCSC.mm10_1.4.3.tar.gz", repos = NULL, type = "source")
mm10_mchrs <- "BSgenome.Mmusculus.UCSC.mm10.mainChrs"
mm10_mchrs_fa <- "BSgenome.Mmusculus.UCSC.mm10.mainChrs.fa"

# Usamos los cromosomas principales del genoma y excluimos contigs que no estén mapeadas
mainChromosomes <- paste0("chr", c(1:19, "X", "Y", "M"))

mainChrSeq <- lapply(mainChromosomes, function(x) BSgenome.Mmusculus.UCSC.mm10[[x]])

names(mainChrSeq) <- mainChromosomes

# Creamos un objeto DNAStringSet con las secuencias recuperadas
mainChrSeqSet <- DNAStringSet(mainChrSeq)
mainChrSeqSet
## DNAStringSet object of length 22:
##          width seq                                          names               
##  [1] 195471971 NNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNN chr1
##  [2] 182113224 NNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNN chr2
##  [3] 160039680 NNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNN chr3
##  [4] 156508116 NNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNN chr4
##  [5] 151834684 NNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNN chr5
##  ...       ... ...
## [18]  90702639 NNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNN chr18
## [19]  61431566 NNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNN chr19
## [20] 171031299 NNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNN chrX
## [21]  91744698 NNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNN chrY
## [22]     16299 GTTAATGTAGCTTAATAACAA...TACGCAATAAACATTAACAA chrM
# Creamos un archivo fasta con las secuencias con las que vamos a alinear nuestras lecturas
writeXStringSet(mainChrSeqSet, mm10_mchrs_fa)

# Construimos un índice para nuestro genoma de referencia
buildindex(paste0(file_dir, "mm10_mainchrs"), paste0(file_dir, mm10_mchrs_fa), 
memory = 8000, indexSplit = TRUE)

# Alineamos nuestros datos filtrados al número fichero fasta con la secuencia del genoma de referencia y creamos un archivo BAM
myMapped <- align("mm10_mainchrs", "filtered_chipseq.fastq.gz", output_format = "BAM", 
output_file = "Myc_Mel_1.bam", type = "dna", phredOffset = 64, nthreads = 8)

# Ordenamos e indexamos archivo BAM
sortBam("Myc_Mel_1.bam", "SR_Myc_Mel_rep1")
indexBam("SR_Myc_Mel_rep1.bam")
# Representamos número de lecturas mapeadas con los índices del archivo BAM
mappedReads <- idxstatsBam("SR_Myc_Mel_rep1.bam")
TotalMapped <- sum(mappedReads[, "mapped"])

ggplot(mappedReads, aes(x = seqnames, y = mapped)) + geom_bar(stat = "identity") +
    coord_flip()

Por último, podemos crear un archivo BigWig, que contiene las coordenadas y el nivel de cobertura en el genoma, para poder visualizarlo en IGV.

# Representamos número de lecturas mapeadas con los índices del archivo BAM
forBigWig <- coverage("SR_Myc_Mel_rep1.bam")

export.bw(forBigWig, con = "SR_Myc_Mel_rep1.bw")

Bibliografía