SCRIPT; Encontrar anotacion de secuencias de interes

+---
+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()  
+```
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+