<jupyter_start><jupyter_text>Gene Ontology enrichment analysisThis is a Jupyter notebook.To run all cells in the notebook use `Cell --> Run All`.To run cells one at a time click into the first code cell and key `Shift-Enter` in each cell in sequence.More information on Jupyter notebooks can be found[here](http://jupyter-notebook.readthedocs.io/en/stable/examples/Notebook/Running%20Code.html).<jupyter_code>## Set the plot window sizes within the notebook (in inches)
options(repr.plot.width=6, repr.plot.height=8)<jupyter_output><empty_output><jupyter_text>Please comment or uncomment the following declarations as appropriate for your run of this notebook:<jupyter_code># Set parameters
opt <- list()
opt$universe   <- './user_data/CMT_PEPs.csv'# File with all genes & their scores. 'interesting' genes chosen based on opt$cutoff
#opt$mapping    <- 'YOURMAPPINGFILE' # Comment this out to have topGO retrieve the latest mapping
opt$outdir     <- './user_data/results'
opt$cutoff     <- 0.05  # Minimum qvalue to be included in a PEP (default 0.05)
opt$num        <- 10 # Minimum number of genes in a geneset for the geneset to be included in enrichment analysis
opt$two        <- FALSE # Set to TRUE to run two-way histology enrichment
opt$workingdir <- './user_data'<jupyter_output><empty_output><jupyter_text>Setting up the notebook libraries and functions<jupyter_code>## Load libraries
print('Loading packages...')
library("topGO")
library("goseq")
library("dplyr")

## If the output directory doesn't exist, create it
if(!dir.exists(opt$outdir)) {
  print(paste('Creating output directory',opt$outdir))
  system(paste('mkdir -p',opt$outdir))
}
### Function to run topGO
run_topGO <- function(peps, pep.name, alg='weight01', ns=opt$num) {

  ## Whole gene list
  geneUniverse <- peps[, pep.name]
  names(geneUniverse) <- peps$HumanSymbol

  ## Load the GO annotations to use
  if( !is.null(opt$mapping) ) {
    print('mapping null')
    geneID2GO <- try(readMappings(file=opt$mapping))
    if(inherits(geneID2GO, "try-error")) { print('ERROR: Unable to read mapping file.'); quit(save='no',status=1) }
  } else {
    geneID2GO = getgo(peps$HumanSymbol, genome='hg19', 'geneSymbol')
  }

  ## How we select the genes of interest
  topDiffGenes <- function(allScore) { return(allScore < opt$cutoff) }

  ## Set the ontology (in this case we're always using BP)
  desiredOntology <- 'BP'

  ## Run topGO
  topGOdata <- new('topGOdata',
    description = paste(pep.name,desiredOntology),
    ontology = desiredOntology,
    allGenes = geneUniverse,
    geneSel = topDiffGenes,
    annot = annFUN.gene2GO,
    gene2GO = geneID2GO,
    nodeSize = ns)

  ## Calculate significance of results, adjust for multiple hypotheses
  resultKS     <- runTest(topGOdata, algorithm = alg, statistic = 'ks')
  dfresult     <- GenTable(topGOdata,pvalue=resultKS,topNodes=length(resultKS@score))
  dfresult$FDR <- p.adjust(dfresult$pvalue,method='BH')
  dfresult     <- dfresult[dfresult$pvalue<0.05,]

  ## Don't return results that are significantly under-represented
  dfresult <- dfresult[ dfresult$Significant > dfresult$Expected,]

  ## Add list of genes in each set
  dfresult$Genes <- sapply(dfresult$GO.ID, function(x) { paste(unlist(genesInTerm(topGOdata, x)), collapse=',') } )

  return( dfresult )
} # End run_topGO<jupyter_output><empty_output><jupyter_text>Run gene ontology enrichment on each PEP list<jupyter_code>### Run topGO on each PEP
peps <- read.csv(opt$universe)

print(paste('Successfully loaded gene universe and scores:',opt$universe))

pep.names <- c('Tumor_Expression_Pattern', 'Carcinoma_Expression_Pattern', 'Adenoma_Expression_Pattern')
go.res <- lapply(pep.names, function(x){ run_topGO(peps,pep.name=x) })
names(go.res) <- pep.names<jupyter_output><empty_output><jupyter_text>Run gene ontology enrichment on BRCA<jupyter_code>print('Running the BRCA comparisons'); flush.console()

## Scores generated in DESeq3.R and also provided in the repository
de.brca <- read.table('./data/BRCA_sig_genes.csv', sep=',', header=TRUE, row.names=1, check.names=FALSE)
geneUniverse <- de.brca[,'AbsLog10FoldChange' ]
names(geneUniverse) <- rownames(de.brca)

geneID2GO = getgo(rownames(de.brca), genome='hg19', 'geneSymbol')
topDiffGenes <- function(allScore) { return(allScore < 0.5) }
desiredOntology <- 'BP'

topGOdata <- new('topGOdata',
  description = paste('BRCA',desiredOntology),
  ontology = desiredOntology,
  allGenes = geneUniverse,
  geneSel = topDiffGenes,
  annot = annFUN.gene2GO,
  gene2GO = geneID2GO,
  nodeSize = opt$num)

resultKS     <- runTest(topGOdata, algorithm = 'weight01', statistic = 'ks')
brca.go.res     <-GenTable(topGOdata,pvalue=resultKS,topNodes=length(resultKS@score))
brca.go.res$FDR <-p.adjust(brca.go.res$pvalue,method='BH')
brca.go.res     <- brca.go.res[brca.go.res$pvalue<0.05,]
brca.go.res <- brca.go.res[ brca.go.res$Significant > brca.go.res$Expected,]
brca.go.res$Genes <- sapply(brca.go.res$GO.ID, function(x) { paste(unlist(genesInTerm(topGOdata, x)), collapse=',') } )

write.table(brca.go.res, file=paste(opt$outdir,'BRCA_GO_res.csv',sep='/'), sep=',', col.names=TRUE, row.names=TRUE, quote=TRUE) # TODO temp print statement<jupyter_output><empty_output><jupyter_text>If the option is set, run gene ontology enrichment on each pairwise histology comparison<jupyter_code>### Run topGO using each histology contrast m_n from LRTtidied
if( opt$two ) {
  print('Running the Histology comparisons.')

  ## Check to make sure the required rda files exist, if yes then load them & run analysis
  ## Make sure all of the required files exist - quit if any are missing
  for( rda in c('LRTtidied.rda','humanmapping.rda') ) {
    rda <- paste(opt$workingdir,rda, sep='/')
    if(!file.exists(rda)) { print(paste('ERROR: Unable to locate',rda)); quit(save='no',status=1) }
    load(rda)
  }

  ## Add human mappings
  LRTtidied$HumanSymbol <- NA #Not all map to human
  LRTtidied$HumanSymbol[ LRTtidied$gene %in% Map_CanEns2HumSymb_unique$Can_Ens ] <- Map_CanEns2HumSymb_unique$Hum_Symb
  LRTtidied <- LRTtidied[! is.na(LRTtidied$HumanSymbol),]

  # Malignant vs Normal - only print significantly enriched terms (not significantly under-enriched)
  print('Processing: Malignant vs Normal')
  dfresult <- run_topGO( as.data.frame(LRTtidied[ LRTtidied$contrast=='m_n',]),  pep.name='logFC')
  res.maligVSnormal <- dfresult # Save for later compiling
  write.table(dfresult, file=paste(opt$outdir,'2WAY_MvsN.csv',sep='/'), sep=',', col.names=TRUE, row.names=FALSE, quote=TRUE)

  # Malignant vs Benign
  print('Processing: Malignant vs Benign')
  dfresult <- run_topGO( as.data.frame(LRTtidied[ LRTtidied$contrast=='m_b',]),  pep.name='logFC')
  write.table(dfresult, file=paste(opt$outdir,'2WAY_MvsB.csv',sep='/'), sep=',', col.names=TRUE, row.names=FALSE, quote=TRUE)


  # Benign vs Normal
  print('Processing: Benign vs Normal')
  dfresult <- run_topGO( as.data.frame(LRTtidied[ LRTtidied$contrast=='b_n',]),  pep.name='logFC')
  write.table(dfresult, file=paste(opt$outdir,'2WAY_BvsN.csv',sep='/'), sep=',', col.names=TRUE, row.names=FALSE, quote=TRUE)
}<jupyter_output><empty_output><jupyter_text>Create the supplemental table with the results<jupyter_code>## Combine the results into 1 table, save to file
# P-val normal-carcinoma, Pvalue Tumor PEP, Pvalue Carcinoma PEP, Pvalue Adenoma PEP, SigGenes in Tumor PEP, SigGenes in Carcinoma PEP, SigGenes in Adenoma PEP
print('Generating the supplemental table with results.')

## Drop to columns we're interested in
go.res$Adenoma_Expression_Pattern   <- go.res$Adenoma_Expression_Pattern[,c('GO.ID','Term','pvalue','Genes')]
go.res$Carcinoma_Expression_Pattern <- go.res$Carcinoma_Expression_Pattern[,c('GO.ID','Term','pvalue','Genes')]
go.res$Tumor_Expression_Pattern     <- go.res$Tumor_Expression_Pattern[,c('GO.ID','Term','pvalue','Genes')]

## Change to unique column names
colnames(go.res$Adenoma_Expression_Pattern)   <- c('GO.ID','Term','PValueEnrich_in_Adenoma_Expression_Pattern','SigGenes_Adenoma_PEP')
colnames(go.res$Carcinoma_Expression_Pattern) <- c('GO.ID','Term','PValueEnrich_in_Carcinoma_Expression_Pattern','SigGenes_Carcinoma_PEP')
colnames(go.res$Tumor_Expression_Pattern)     <- c('GO.ID','Term','PValueEnrich_in_Tumor_Expression_Pattern','SigGenes_Tumor_PEP')

## Combine the PEP results tables into 1
res <- full_join(go.res$Tumor_Expression_Pattern, go.res$Carcinoma_Expression_Pattern, by=c('GO.ID','Term'))
res <- full_join(go.res$Adenoma_Expression_Pattern, res, by=c('GO.ID','Term'))

## Add Normal vs Carcinoma comparison if option is set
## Either way, reorder columns based on which comparisons are included
if( opt$two ) {
  res.maligVSnormal <- res.maligVSnormal[,c('GO.ID','Term','pvalue')]
  colnames(res.maligVSnormal) <- c('GO.ID','Term','PValueEnrich_in_Normal_vs_Carcinoma')
  res <- full_join(res.maligVSnormal, res, by=c('GO.ID','Term'))
  res <- res[, c('GO.ID','Term','PValueEnrich_in_Normal_vs_Carcinoma','PValueEnrich_in_Tumor_Expression_Pattern','PValueEnrich_in_Carcinoma_Expression_Pattern','PValueEnrich_in_Adenoma_Expression_Pattern','SigGenes_Tumor_PEP','SigGenes_Carcinoma_PEP','SigGenes_Adenoma_PEP')]
} else {
  res <- res[, c('GO.ID','Term','PValueEnrich_in_Tumor_Expression_Pattern','PValueEnrich_in_Carcinoma_Expression_Pattern','PValueEnrich_in_Adenoma_Expression_Pattern','SigGenes_Tumor_PEP','SigGenes_Carcinoma_PEP','SigGenes_Adenoma_PEP')]
}

## NA gene lists to empty string
res$SigGenes_Tumor_PEP[is.na(res$SigGenes_Tumor_PEP)] <- ''
res$SigGenes_Carcinoma_PEP[is.na(res$SigGenes_Carcinoma_PEP)] <- ''
res$SigGenes_Adenoma_PEP[is.na(res$SigGenes_Adenoma_PEP)] <- ''

## NA p-values to 1
res$PValueEnrich_in_Tumor_Expression_Pattern[is.na(res$PValueEnrich_in_Tumor_Expression_Pattern)] <- 1
res$PValueEnrich_in_Carcinoma_Expression_Pattern[is.na(res$PValueEnrich_in_Carcinoma_Expression_Pattern)] <- 1
res$PValueEnrich_in_Adenoma_Expression_Pattern[is.na(res$PValueEnrich_in_Adenoma_Expression_Pattern)] <- 1
if( opt$two ) { res$PValueEnrich_in_Normal_vs_Carcinoma[is.na(res$PValueEnrich_in_Normal_vs_Carcinoma)] <- 1 }

## Sort by Tumor P-value
res <- res[with(res, order(PValueEnrich_in_Tumor_Expression_Pattern)),]

## Save to file
write.table(res, file=paste(opt$outdir, 'GO_Enrichment.csv', sep='/'), sep=',', col.names=T, row.names=F, quote=T)

print(paste('Finished, printing files to:',opt$outdir))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13930013458950202, 'fraction_numerical': 0.005383580080753701, 'mean_word_length': 6.2343532684283725, 'pattern_counts': {'":': 0, '<': 99, '<?xml version=': 0, '>': 32, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '373548', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 3919, 'n_tokens_neox': 3657, 'n_words': 742, 'path': 'DataAnalysis/topGO Notebook.ipynb', 'repo_name': 'flatironinstitute/FREYA', 'size': '10402'}