---
title: "GH"
author: "EQUIPO 5"
date: "22/4/2020"
output: html_document
---

Este script tiene la finalidad de obtener la anotacion en el genoma de Homo sapiens de aquellas coordenadas que corresponden a los genes de interes.

```{r}
#se mandan a llamar a todas las librerias a utilizar 
library(Biobase)
library(IRanges)
library(rtracklayer)
library(GenomicRanges)
library(Rsamtools)
library(ggplot2)
```

#Obtenemos la anotación de un genoma
Para tener con que comparar el genoma que vamos a analizar, es necesario tener las coordenadas y la cadena en la que se encuentra las regiones del genoma anotadas, para esto ensembl tiene estos archivos.       

  ```{r}
  #se descarga el genoma de referencia con sus anotaciones, en este caso fueron obtenidos por ensembl100
setwd("/home/aschafer/Documentos/Genomicas/Semestre_4/Genomica_humana/Enfermedades_monogenicas")
homoS = import("Homo_sapiens.GRCh38.100.gff3.gz")
head(homoS)
```

```{r}
#visualización de los datos que tiene el genoma de referencia. 
table(mcols(homoS)$type)
```

```{r}
#nos quedamos solo con las columnas que deseamos.
mcols(homoS) = mcols(homoS)[,c("source","type","ID","Name")]
```


#Importamos nuestros datos a analizar
Despues de haber trabajado nuestros datos con bowtie, utilizamos el archivo resultante .bam    

```{r}
#cargamos nuestros datos de las enfermedades monogenicas.
bamFile <- "/home/aschafer/Documentos/Genomicas/Semestre_4/Genomica_humana/Enfermedades_monogenicas/sequences_aligned_A_sort.bam"
informacion <- c("rname", "strand", "pos", "qwidth")
informacion <- ScanBamParam(what = informacion)
bam <- scanBam(bamFile, param = informacion)
lapply(bam, names)
```


```{r}
#Construimos un GRanges
mapGR <- GRanges(seqnames = bam[[1]]$rname, ranges = IRanges(start = bam[[1]]$pos, width = bam[[1]]$qwidth), strand = bam[[1]]$strand)
head(mapGR)
```


#Contamos mapeos dentro de las anotaciones 
```{r}
#buscamos en que regiones se traslapan 
traslapados <- countOverlaps(homoS, mapGR)
typeCounts <- aggregate(traslapados, by = list("type"=mcols(homoS)$type), sum)
head(typeCounts)
```

```{r}
#vamos a graficar todos los resultados obtenidos, pero para eso vamos a juntas todos aquellos que tengan menos de 1000 juntos. 
valMin <- 1000
Cuentas <- typeCounts[typeCounts$x > valMin,]
Cuentas <- Cuentas[order(Cuentas$x),]
Cuentas_totales <- rbind(data.frame("type" = "other", "x" = sum(typeCounts$x[typeCounts$x <= valMin])), Cuentas)
porcentaje <- round(100*Cuentas_totales$x/sum(Cuentas_totales$x), 1)

#png(filename="/home/aschafer/Documentos/Genomicas/Semestre_4/Genomica_humana/Enfermedades_monogenicas/Grafica1.png", width = 550, height = 400)
pie(Cuentas_totales$x, labels = paste(porcentaje, "%", sep = ""), col = rev(rainbow(nrow(Cuentas_totales))), main="Porcentaje de lecturas totales alineadas por tipo", cex=0.9)
legend("topleft", legend = Cuentas_totales$type, cex = 0.8, fill = rev(rainbow(nrow(Cuentas_totales))))
#dev.off()
```
Debido que vimos que los tipos: cromosoma, región biológica, exon, mRNA, gen, cinco prima UTR. No eran tan significativos, decidimos eliminarlos para los próximos análisis. 

```{r}
#graficamos ahora los nuevos datos, pero ahora en lugar de utilizar 1000 como el valor minimo, utilizamos 100.
nuevosDatos <- typeCounts[typeCounts$type != "chromosome" & typeCounts$type != "biological_region" & typeCounts$type != "exon" & typeCounts$type != "mRNA" & typeCounts$type != "gene" & typeCounts$type != "five_prime_UTR",]
colores <- c("violetred1", "dodgerblue4", "green2", "yellow", "red",  "steelblue1")
valMin2 <- 100
nuevas_cuentas <- nuevosDatos[nuevosDatos$x > valMin2,]
nuevas_cuentas <- nuevas_cuentas[order(nuevas_cuentas$x),]
nuevas_cuentas_totales <- rbind(data.frame("type"="other", "x" = sum(nuevosDatos$x[nuevosDatos$x <= valMin2])), nuevas_cuentas)
porcentaje2 <- round(100*nuevas_cuentas_totales$x/sum(nuevas_cuentas_totales$x), 1)
  
#png(filename="/home/aschafer/Documentos/Genomicas/Semestre_4/Genomica_humana/Enfermedades_monogenicas/Grafica2.png", width = 550, height = 400)
pie(nuevas_cuentas_totales$x, labels = paste(porcentaje2, "%", sep = ""), cex = 1, main = "Porcentaje de lecturas alineadas por tipo seleccionados", col = colores)
legend("topleft", legend = nuevas_cuentas_totales$type, cex = 0.9, fill = colores)
#dev.off()
```

```{r}
# Quitamos las categorías de las antocaciones que son redundantes a la clasificación codificante-no codificante y graficamos ahora los nuevos datos, pero ahora en lugar de utilizar 1000 como el valor minimo, utilizamos 100.
nuevosDatos <- typeCounts[-c(1,2,6,9,11,13,16,24),]
colores <- c("violetred1", "orangered1", "dodgerblue4", "green2", "yellow", "purple", "red",  "steelblue1")
valMin2 <- 100
nuevas_cuentas <- nuevosDatos[nuevosDatos$x > valMin2,]
nuevas_cuentas <- nuevas_cuentas[order(nuevas_cuentas$x),]
nuevas_cuentas_totales <- rbind(data.frame("type"="other", "x" = sum(nuevosDatos$x[nuevosDatos$x <= valMin2])), nuevas_cuentas)
porcentaje2 <- round(100*nuevas_cuentas_totales$x/sum(nuevas_cuentas_totales$x), 1)
total <- sum(nuevosDatos$x) #sacamos el total de lecturas
#establecer el número de lecturas que pertenecen a categorías codificantes
codificante <- nuevosDatos[c(3,6,7),] 
codificante <- sum(codificante$x)
#ahora para las no codificantes
no_codificante <- total - codificante
datos_finales<-data.frame("type" = c("codificante", "no codificante"), "x" = c(codificante, no_codificante))
colores2<-c("slateblue3", "green2")
porcentaje3 <- round(100*datos_finales$x/sum(datos_finales$x), 1)
  
#png(filename="/home/aschafer/Documentos/Genomicas/Semestre_4/Genomica_humana/Enfermedades_monogenicas/Grafica2.png", width = 550, height = 400)
#par(mfrow = c(1,2))
pie(nuevas_cuentas_totales$x, labels = paste(porcentaje2, "%", sep = ""), cex = 0.9, main = "Porcentaje de cada clasificación", col = colores, cex.main = 1)
legend("bottomleft", legend = nuevas_cuentas_totales$type, cex = 1, fill = colores )
#dev.off()
```




Datos recuperados por OMIM
```{r}
no_coding <- sum(20, 15, 9, 7, 4, 3, 2, 2, 2, 2, 1, 1, 1, 1, 1, 228, 144)
protein_coding <- 14886
OMIM<- data.frame("type" = c("Regiones codificantes", "Regiones no codificantes"), "x" = c(protein_coding, no_coding))
```

```{r}
#png(filename="/home/aschafer/Documentos/Genomicas/Semestre_4/Genomica_humana/Enfermedades_monogenicas/Grafica4.png", width = 550, height = 400)
OMIM$x
porcentaje4 <- round(100*as.numeric(OMIM$x)/sum(as.numeric(OMIM$x)), 1)
pie(as.numeric(OMIM$x), labels = paste(porcentaje4, "%", sep = ""), cex = 0.9, main = "Porcentaje de regiones codificantes y no codificantes de acuerdo a OMIM y a DISEASE", col = colores2, cex.main = 0.95)
legend("bottomleft", legend = OMIM$type, cex = 0.8, fill = colores2)
#dev.off()  
```