Script; obtener anotacion de zonas de interes

+# ENFERMEDADES MONOGENICAS
+# EQUIPO 5
+
+# 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.
+
+
+#se mandan a llamar a todas las librerias a utilizar 
+library(Biobase)
+library(IRanges)
+library(rtracklayer)
+library(GenomicRanges)
+library(Rsamtools)
+
+
+#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.       
+
+#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)
+
+#visualización de los datos que tiene el genoma de referencia. 
+table(mcols(homoS)$type)
+
+#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    
+
+#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)
+
+#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 
+#buscamos en que regiones se traslapan 
+traslapados <- countOverlaps(homoS, mapGR)
+typeCounts <- aggregate(traslapados, by = list("type"=mcols(homoS)$type), sum)
+head(typeCounts)
+
+#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)
+
+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))))
+
+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. 
+
+
+#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)
+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)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+