permuted_gene_sets_tcga_padog.R

# scriptDir=/wynton/group/gladstone/biocore/projects/pfocr_pathway_enrichment_evaluation/PFOCRInPathwayAnalyses/
# containerDir=/wynton/group/gladstone/biocore/containers
# dataDir=/wynton/group/gladstone/biocore/projects/PFOCR/
# dataDir1=/wynton/group/gladstone/biocore/projects/pfocr_pathway_enrichment_evaluation/permuted_geneset_databases_tcga_gsea_results_july_2021/
# export SINGULARITY_BINDPATH="$containerDir,$scriptDir,$dataDir,$dataDir1"
# singularity exec $containerDir/pfocr_path_analyses_tools_latest.sif R

##################################required libraries
require(rSEA)
require(GSEABenchmarkeR)
require("clusterProfiler")
require("org.Hs.eg.db")
require(dplyr)
require(GO.db)
require(dplyr)
require(magrittr)
require(BiocParallel)

register(BatchtoolsParam(workers = 8), default = TRUE)
bp.par <- BiocParallel::registered()[[1]]

args <- commandArgs(trailingOnly=TRUE)

GSE_index <- as.integer(args[1])
print(GSE_index)


tcga <- loadEData("tcga", nr.datasets=40, cache = FALSE)
tcga <- tcga[GSE_index]

##run the DE analyses

tcga %<>% runDE(., de.method="DESeq2", padj.method="BH")


##load gene set databases
outdir <- "/wynton/group/gladstone/biocore/projects/pfocr_pathway_enrichment_evaluation/permuted_geneset_databases_tcga_padog_results_july_2021/"

min_set_size <- 3
max_set_size <- 500 
dataDir <- "/wynton/group/gladstone/biocore/projects/PFOCR/PFOCRInPathwayAnalyses/"
database_lists1 <- load("/wynton/group/gladstone/biocore/projects/PFOCR/PFOCRInPathwayAnalyses/databases_pfocr_3intersect_v2.RData")#pfocr_3sets
database_lists2 <- load("/wynton/group/gladstone/biocore/projects/PFOCR/PFOCRInPathwayAnalyses/databases.RData")#has wp, pfocr, go
database_lists3 <- load("/wynton/group/gladstone/biocore/projects/PFOCR/PFOCRInPathwayAnalyses/pfocr_miny_April2021.RData")#has wp, pfocr, go

database_lists <- append(database_lists1, database_lists2)
database_lists <- append(database_lists, database_lists3)

database_annotation <- unname(unlist(sapply(database_lists, grep, pattern="annotation", value = T, perl = T)))
database_lists <- unname(unlist(sapply(database_lists, grep, pattern="list$", value = T, perl = T)))
for (db in database_lists) {
  eval(call("<-", as.name(db),  Filter(Negate(is.null), lapply(get(db), function(x){
    if(length(x) < min_set_size | length(x) > max_set_size)
      NULL
    else
      x
  }))
  ))
}

##need to filter annotation data as well
wp_annotation %<>% filter(set_id %in% names(wp_list))
go_annotation %<>% filter(set_id %in% names(go_list))
pfocr_annotation %<>% filter(set_id %in% names(pfocr_list))
pfocr_annotation_3sets %<>% filter(set_id %in% names(pfocr_3sets_list))
pfocr_miny_annotation %<>% filter(set_id %in% names(pfocr_miny_list))


gene_entrez=readRDS("/wynton/group/gladstone/biocore/projects/PFOCR/PFOCRInPathwayAnalyses/gene_entrez.rds")

set.seed(1234)

Nperm <- 100
qLevels <- c(0.5, 0.75, 0.9, 0.95, 0.99)
names_qLevels <- paste("quantile", round(100*qLevels), sep="_")


choose_random_genes_in_path <- function(path, path_annotation, gene_prob) {
  
  random_genes <- sample(names(gene_prob), size=length(path), replace = FALSE, prob = gene_prob)
  
  return(random_genes)
}

PermuteDatabase <- function(p, GSE_index, path_list, path_annotation,tcga, gene_entrez,path_gene_prob) {
  #print(p)
  temp_annotation <- path_annotation
  temp_list <- lapply(path_list, choose_random_genes_in_path, path_annotation, path_gene_prob)
  for(i in 1:length(temp_list)) {
    temp_annotation$gene[temp_annotation$set_id == names(temp_list)[i]] <- temp_list[[i]]
  }
  
  gsea.res <- runEA(tcga[[1]], method="padog", gs=temp_list, perm=1000, parallel = bp.par)
  res <- gsea.res$padog[[1]]$ranking
  
  Nsig <- sum(res$PVAL < 0.01, na.rm = TRUE)
  TempRes <- data.frame(p.value=res$PVAL, aNES=res$PADOG0)
  TempRes %<>% filter(p.value < 0.01)
  TDP_bound_90 <- TempRes %>%
    .$aNES %>%
    quantile(.,c(0.5, 0.75, 0.9, 0.95, 0.99), na.rm=TRUE)
  Res <- append(Nsig, TDP_bound_90)
  names(Res) <- append("Nsig", names_qLevels)
  print(Res)
  return(Res)
}


##WP
gsea.res.wp <- runEA(tcga[[1]], method="padog", gs=wp_list, perm=1000, parallel = bp.par)
res.wp <- gsea.res.wp$padog[[1]]$ranking

TempRes <-  data.frame(p.value=res.wp$PVAL, aNES=res.wp$PADOG0)
TempRes %<>% filter(p.value < 0.01)
TDP_bound_90 <- TempRes %>%
  .$aNES %>%
  quantile(.,c(0.5, 0.75, 0.9, 0.95, 0.99), na.rm=TRUE)
oTDP_bound_90_wp <- TDP_bound_90
oNsig_wp <- sum(res.wp$PVAL < 0.01, na.rm = TRUE)
oTDP_bound_90_wp <- TDP_bound_90
print(oNsig_wp)

unique_genes <- wp_annotation %>%
  .$gene %>%
  unique()

wp_gene_prob <- wp_annotation %>%
  .$gene %>%
  table() %>%
  divide_by(length(unique_genes))

rNsig_wp <- t(sapply(1:Nperm, PermuteDatabase, GSE_index, wp_list, wp_annotation,tcga, gene_entrez, wp_gene_prob))

##GO

gsea.res.go <- runEA(tcga[[1]], method="padog", gs=go_list, perm=1000, parallel = bp.par)
res.go <- gsea.res.go$padog[[1]]$ranking

TempRes <-  data.frame(p.value=res.go$PVAL, aNES=res.go$PADOG0)
TempRes %<>% filter(p.value < 0.01)

TDP_bound_90 <- TempRes %>%
  .$aNES %>%
  quantile(.,c(0.5, 0.75, 0.9, 0.95, 0.99), na.rm=TRUE)
oTDP_bound_90_go <- TDP_bound_90

oNsig_go <- sum(res.go$PVAL < 0.01, na.rm = TRUE)
print(oNsig_go)

unique_genes <- go_annotation %>%
  .$gene %>%
  unique()

go_gene_prob <- go_annotation %>%
  .$gene %>%
  table() %>%
  divide_by(length(unique_genes))

rNsig_go <- t(sapply(1:Nperm, PermuteDatabase, GSE_index, go_list, go_annotation,tcga, gene_entrez, go_gene_prob))

##PFCOR-3sets
gsea.res.pfocr <- runEA(tcga[[1]], method="padog", gs=pfocr_3sets_list, perm=1000, parallel = bp.par)
res.pfocr <- gsea.res.pfocr$padog[[1]]$ranking

TempRes <-  data.frame(p.value=res.pfocr$PVAL, aNES=res.pfocr$PADOG0)
TempRes %<>% filter(p.value < 0.01)

TDP_bound_90 <- TempRes %>%
  .$aNES %>%
  quantile(.,c(0.5, 0.75, 0.9, 0.95, 0.99), na.rm=TRUE)
oTDP_bound_90_pfocr <- TDP_bound_90

unique_genes <- pfocr_annotation_3sets %>%
  .$gene %>%
  unique()

pfocr_3sets_gene_prob <- pfocr_annotation_3sets %>%
  .$gene %>%
  table() %>%
  divide_by(length(unique_genes))

oNsig_pfocr <- sum(res.pfocr$PVAL < 0.01, na.rm = TRUE)
print(oNsig_pfocr)
rNsig_pfcor <- t(sapply(1:Nperm, PermuteDatabase, GSE_index, pfocr_3sets_list, pfocr_annotation_3sets,tcga, gene_entrez, pfocr_3sets_gene_prob))

##PFCOR-miny
gsea.res.pfocr.miny <- runEA(tcga[[1]], method="padog", gs=pfocr_miny_list, perm=1000, parallel = bp.par)
res.pfocr.miny <- gsea.res.pfocr.miny$padog[[1]]$ranking

TempRes <-  data.frame(p.value=res.pfocr.miny$PVAL, aNES=res.pfocr.miny$PADOG0)
TempRes %<>% filter(p.value < 0.01)

TDP_bound_90 <- TempRes %>%
  .$aNES %>%
  quantile(.,c(0.5, 0.75, 0.9, 0.95, 0.99), na.rm=TRUE)
oTDP_bound_90_pfocr_miny <- TDP_bound_90

oNsig_pfocr_miny <- sum(res.pfocr.miny$PVAL < 0.01, na.rm = TRUE)
print(oNsig_pfocr_miny)

unique_genes <- pfocr_miny_annotation %>%
  .$gene %>%
  unique()

pfocr_miny_gene_prob <- pfocr_miny_annotation %>%
  .$gene %>%
  table() %>%
  divide_by(length(unique_genes))

rNsig_pfcor_miny <- t(sapply(1:Nperm, PermuteDatabase, GSE_index, pfocr_miny_list, pfocr_miny_annotation,tcga, gene_entrez, pfocr_miny_gene_prob))

GSE_index_2_100_perm_res <- list(oNsig_wp=oNsig_wp, rNsig_wp=rNsig_wp, oNsig_go=oNsig_go, rNsig_go=rNsig_go, oNsig_pfocr=oNsig_pfocr, rNsig_pfcor=rNsig_pfcor, oNsig_pfocr_miny=oNsig_pfocr_miny, rNsig_pfcor_miny=rNsig_pfcor_miny, oTDP_bound_90_wp=oTDP_bound_90_wp, oTDP_bound_90_go=oTDP_bound_90_go, oTDP_bound_90_pfocr=oTDP_bound_90_pfocr, oTDP_bound_90_pfocr_miny=oTDP_bound_90_pfocr_miny)
saveRDS(GSE_index_2_100_perm_res, file=paste0(outdir, "GSE_index_",GSE_index,"_",Nperm,"_padog_tcga_perm_result.rds"))