#Codigo corrido en R 3.4.1
#Fijar el directorio de trabajo en la carpeta donde esta el archivo descargado y descomprimido
setwd("ruta de la carpeta descargada")

source("https://bioconductor.org/biocLite.R")
biocLite("clusterProfiler")
biocLite("pathview")
biocLite("GOSemSim")
biocLite("KEGG.db")
biocLite("STRINGdb")
biocLite("DOSE") #Warning message:package 'DOSE' was built under R version 3.2.4
biocLite("GO.db")
biocLite("org.Hs.eg.db")
biocLite("topGO")
biocLite("GSEABase")



library(org.Hs.eg.db)
keytypes(org.Hs.eg.db)

library(clusterProfiler)
#Traductor
x <- c("GPX3",  "GLRX",   "LBP",   "CRYAB", "DEFB1", "HCLS1",   "SOD2",   "HSPA2",
       "ORM1",  "IGFBP1", "PTHLH", "GPC3",  "IGFBP3","TOB1",    "MITF",   "NDRG1",
       "NR1H4", "FGFR3",  "PVR",   "IL6",   "PTPRM", "ERBB2",   "NID2",   "LAMB1",
       "COMP",  "PLS3",   "MCAM",  "SPP1",  "LAMC1", "COL4A2",  "COL4A1", "MYOC",
       "ANXA4", "TFPI2",  "CST6",  "SLPI",  "TIMP2", "CPM",     "GGT1",   "NNMT",
       "MAL",   "EEF1A2", "HGD",   "TCN2",  "CDA",   "PCCA",    "CRYM",   "PDXK",
       "STC1",  "WARS",  "HMOX1", "FXYD2", "RBP4",   "SLC6A12", "KDELR3", "ITM2B")
eg = bitr(x, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")
head(eg)

ids <- bitr(x, fromType="SYMBOL", toType=c("UNIPROT", "ENSEMBL"), OrgDb="org.Hs.eg.db")
head(ids)

#KEGG
data(gcSample)
hg <- gcSample[[1]]
head(hg)
eg2np <- bitr_kegg(hg, fromType='kegg', toType='ncbi-proteinid', organism='hsa')
head(eg2np)

#Clasificacion/agrupacion por GO
data(geneList, package="DOSE")
gene <- names(geneList)[abs(geneList) > 2]
head(gene)
ggo <- groupGO(gene     = gene,
               OrgDb    = org.Hs.eg.db,
               ont      = "CC",
               level    = 3,
               readable = TRUE)
head(ggo)
View(ggo)

#GO over-representation analysis
ego <- enrichGO(gene          = gene,
                universe      = names(geneList),
                OrgDb         = org.Hs.eg.db,
                ont           = "BP",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.01,
                qvalueCutoff  = 0.05,
                readable      = TRUE)
head(ego)


#KEGG 
#clusterProfiler provee funcion para ayudar en la busqueda de organismo soportados
search_kegg_organism('ece', by='kegg_code')
ecoli <- search_kegg_organism('Escherichia coli', by='scientific_name')
dim(ecoli)
head(ecoli)

#ORA usando KEGG
kk <- enrichKEGG(gene         = gene,
                 organism     = 'hsa',
                 pvalueCutoff = 0.05)
View(kk)


#Visualizaciones
barplot(ggo, drop=TRUE, showCategory=12)
barplot(ego, showCategory=12)
dotplot(ego)
enrichMap(ego)
## categorySize can be scaled by 'pvalue' or 'geneNum'
cnetplot(ego, categorySize="pvalue", foldChange=geneList)
plotGOgraph(ego)

browseKEGG(kk, 'hsa04110')

library("pathview")
hsa04110 <- pathview(gene.data  = geneList,
                     pathway.id = "hsa04110",
                     species    = "hsa",
                     limit      = list(gene=max(abs(geneList)), cpd=1))



#comparaciones por grupoo
data(gcSample)
lapply(gcSample, head)
ck <- compareCluster(geneCluster = gcSample, fun = "enrichKEGG",use_internal_data = TRUE)
View(as.data.frame(ck))
dotplot(ck)

# compareCluster usando formula personalizada
mydf <- data.frame(Entrez=names(geneList), FC=geneList)
mydf <- mydf[abs(mydf$FC) > 1,]
mydf$group <- "upregulated"
mydf$group[mydf$FC < 0] <- "downregulated"
mydf$othergroup <- "A"
mydf$othergroup[abs(mydf$FC) > 2] <- "B"
formula_res <- compareCluster(Entrez~group+othergroup, data=mydf, fun="enrichKEGG")
View(as.data.frame(formula_res))
dotplot(formula_res)
dotplot(formula_res, x=~group) + ggplot2::facet_grid(~othergroup)




################################
# STRING
library(STRINGdb)
especies<-get_STRING_species(version="10", species_name=NULL)
View(especies) #2031 especies!
string_db <- STRINGdb$new( version="10", species=9606,score_threshold=0, input_directory="" )
#Que podemos hacer
STRINGdb$methods()

data(diff_exp_example1)
head(diff_exp_example1)
#Mapear los genes desde los nombres HUGO
example1_mapped <- string_db$map( diff_exp_example1, "gene", removeUnmappedRows = TRUE )

#Extraemos 50 genes para producir la red.
hits <- example1_mapped$STRING_id[1:50]
#generamos la red
string_db$plot_network( hits )

# filtrar por p-value y agregar columna de color: verde sub, rojo sobreexpresado.
example1_mapped_pval05 <- string_db$add_diff_exp_color( subset(example1_mapped, pvalue<0.05),
                                                        logFcColStr="logFC" )
# lelvar la informacion al servidor
payload_id <- string_db$post_payload( example1_mapped_pval05$STRING_id,
                                      colors=example1_mapped_pval05$color )
# Muestra la red con un halo usando el filtro anterior
string_db$plot_network( hits, payload_id=payload_id )


#Enriquecimiento
## Grafica el enriquecimiento para los mejores 100 genes
string_db$plot_ppi_enrichment( example1_mapped$STRING_id[1:100], quiet=TRUE )
enrichmentGO <- string_db$get_enrichment( hits, category = "Process", methodMT = "fdr", iea = TRUE )
View(enrichmentGO)