# 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)