Skip to content
Snippets Groups Projects
Commit 68a73286 authored by Holger Brandl's avatar Holger Brandl
Browse files

fixed package loading for R3.3

parent 05130a0f
No related branches found
No related tags found
No related merge requests found
......@@ -20,14 +20,14 @@ Options:
opts <- docopt(doc, commandArgs(TRUE)) ## does not work when spining
# opts <- docopt(doc, "--overlay_expr_data ../plot_score_matrix.txt ../degs_by_contrast.txt contrast" )
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.25/R/core_commons.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.25/R/ggplot_commons.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.25/R/bio/diffex_commons.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.26/R/core_commons.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.26/R/ggplot_commons.R")
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.26/R/bio/diffex_commons.R")
loadpack(knitr)
loadpack(DT)
loadpack(clusterProfiler)
#loadpack(clusterProfiler)
#devtools::session_info() # nice!
......@@ -63,7 +63,7 @@ resultsBaseName <- if(str_length(opts$project)>0) paste0(opts$project, ".") else
qCutoff <- as.numeric(opts$qcutoff)
reload_dplyr()
#reload_dplyr()
########################################################################################################################
#' ## Enrichment Analysis
......@@ -139,7 +139,7 @@ annoDb <- guess_anno_db(geneLists$ensembl_gene_id) # e.g. "org.Hs.eg.db"
## convert to entrez gene ids
glMappedRaw <- clusterProfiler::bitr(geneLists$ensembl_gene_id, fromType="ENSEMBL", toType="ENTREZID", annoDb=annoDb) %>%
glMappedRaw <- clusterProfiler::bitr(geneLists$ensembl_gene_id, fromType="ENSEMBL", toType="ENTREZID", OrgDb=annoDb) %>%
set_names("ensembl_gene_id", "entrez_gene_id") %>% right_join(geneLists)
......@@ -147,6 +147,11 @@ glMappedRaw <- clusterProfiler::bitr(geneLists$ensembl_gene_id, fromType="ENSEM
#' Check how many failed to map
count(glMappedRaw, is.na(entrez_gene_id))
unloadNamespace('clusterProfiler')
#loadpack(clusterProfiler)
reload_dplyr()
glMapped <- glMappedRaw %>% filter(!is.na(entrez_gene_id)) %>% select(-ensembl_gene_id) %>%
## regroup the data
# group_by_(cluster)
......@@ -158,11 +163,11 @@ glMapped <- glMappedRaw %>% filter(!is.na(entrez_gene_id)) %>% select(-ensembl_g
#count(glMappedSub, cluster)
## sync reactome pacakge to node
if(Sys.getenv("LSF_SERVERDIR")!="" && dir.exists("/tmp/local_r_packages")){
system("if [ ! -d '/tmp/local_r_packages/reactome.db/' ]; then scp -r falcon:/tmp/local_r_packages /tmp ; fi")
}
#if(Sys.getenv("LSF_SERVERDIR")!="" && dir.exists("/tmp/local_r_packages")){
# system("if [ ! -d '/tmp/local_r_packages/reactome.db/' ]; then scp -r falcon:/tmp/local_r_packages /tmp ; fi")
#}
require_auto(ReactomePA)
loadpack(ReactomePA)
cp_test <- function(geneIds){
# DEBUG geneIds <- glMapped %>% filter(cluster %in% c("cluster_9")) %$% entrez_gene_id %>% as.integer
......@@ -179,11 +184,11 @@ cp_test <- function(geneIds){
# browser()
# pantherResults <- enricher(gene = geneIds, organism = cpSpecies, qvalueCutoff = qCutoff, readable = TRUE, TERM2GENE = PANTHER10_ontology) %>% summary()
keggResults <- enrichKEGG(gene = geneIds, organism = cpSpecies, qvalueCutoff = qCutoff, readable = TRUE, use_internal_data=T) %>% summary()
reactomeResults <- enrichPathway(gene = geneIds, organism = cpSpecies, qvalueCutoff = qCutoff, readable = TRUE) %>% summary()
goResultsCC <- enrichGO(gene = geneIds, organism = cpSpecies, qvalueCutoff = qCutoff, readable = TRUE, ont = "CC") %>% summary()
goResultsMF <- enrichGO(gene = geneIds, organism = cpSpecies, qvalueCutoff = qCutoff, readable = TRUE, ont = "MF") %>% summary()
goResultsBP <- enrichGO(gene = geneIds, organism = cpSpecies, qvalueCutoff = qCutoff, readable = TRUE, ont = "BP") %>% summary()
keggResults <- clusterProfiler::enrichKEGG(gene = geneIds, organism = cpSpecies, qvalueCutoff = qCutoff, use_internal_data=T) %>% summary()
reactomeResults <- enrichPathway(gene = geneIds, organism = cpSpecies, qvalueCutoff = qCutoff) %>% summary()
goResultsCC <- clusterProfiler::enrichGO(gene = geneIds, OrgDb = annoDb, qvalueCutoff = qCutoff, ont = "CC") %>% summary()
goResultsMF <- clusterProfiler::enrichGO(gene = geneIds, OrgDb = annoDb, qvalueCutoff = qCutoff, ont = "MF") %>% summary()
goResultsBP <- clusterProfiler::enrichGO(gene = geneIds, OrgDb = annoDb, qvalueCutoff = qCutoff, ont = "BP") %>% summary()
#cp-bug: if no pathways are enriched odd strucuture is retured ##todo file issue
if(!("data.frame" %in% class(keggResults))) keggResults <- filter(goResultsBP, Description="foobar")
......@@ -211,7 +216,7 @@ enrResults <- quote(glMapped %>% do(cp_test(.$entrez_gene_id))) %>% cache_it(pa
#}, .progress="text", .parallel=F) ##%>% cache_it(paste0("enrdata_", digest(glMapped)))
#rbind_all(enrResults)
reload_dplyr()
#reload_dplyr()
#enrResults %<>% rename(Category=Description)
#enrResults %<>% rename(ontology=Category)
......@@ -223,7 +228,7 @@ write_tsv(enrResults, path=paste0(resultsBaseName, "enrResults.txt"))
# enrResults <- read.delim(paste0(resultsBaseName, "enrResults.txt"))
#' [Enrichment Results](`r paste0(resultsBaseName, "enrResults.txt")`)
require_auto(DT)
loadpack(DT)
datatable(enrResults)
#enrResults %>% ggplot(aes(pvalue)) + geom_histogram() + scale_x_log10()
......@@ -376,10 +381,14 @@ keggPathways <- enrResults %>% filter(ontology=="kegg") %$% ac(ID) %>% unique()
# cat("No enriched pathways found")
#}
loadpack(pathview)
loadpack(png)
loadpack(naturalsort)
#reload_dplyr()
#unloadNamespace('pathview')
# keggPathways <- keggPathways[1]
#if(nrow(keggPathways)==0){
......@@ -413,12 +422,34 @@ dir.create(pwPlotDir)
# set_names("entrez_gene_id", "ensembl_gene_id" ) %>%
# inner_join( sliceData %>% add_rownames("ensembl_gene_id") )
# install.packages("session")
# session::save.session(".pview_debug.dat");
# session::restore.session(".pview_debug.dat");
## from http://stackoverflow.com/questions/7505547/detach-all-packages-while-working-in-r
## todo remove this once we use v1.27 of core_commons
detachAllPackages <- function() {
basic.packages <- c("package:stats","package:graphics","package:grDevices","package:utils","package:datasets","package:methods","package:base")
package.list <- search()[ifelse(unlist(gregexpr("package:",search()))==1,TRUE,FALSE)]
package.list <- setdiff(package.list,basic.packages)
if (length(package.list)>0) for (package in package.list) detach(package, character.only=TRUE)
}
detachAllPackages()
devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.26/R/core_commons.R")
#loadpack(plyr)
#loadpack(dplyr)
#loadpack(dplyr)
#loadpack(stringr)
loadpack(pathview)
loadpack(png)
plot_pathway <- function(pathwayID, sliceData){
# pathwayID="mmu04015"; sliceData <- sliceData
# pathwayID="mmu05145"; sliceData <- sliceData
# pathwayID="mmu05216"; sliceData <- sliceData
# browser()
# echo("processing pathway", pathwayID)
pv.out <- pathview(
......@@ -480,6 +511,8 @@ pathwayPlots <- pathwayPlots[lapply(pathwayPlots, is.character) == 0]
stopifnot(llply(pathwayPlots, function(x) file.exists(x$plotfile)) %>% unlist %>% all)
.instpack("biomaRt")
## prepare tooltips with expression scores
ens2entrez <- quote({
# mart <- biomaRt::useDataset(guess_mart(geneLists$ensembl_gene_id), mart = biomaRt::useMart("ensembl"))
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment