WGCNA: Weighted Correlation Network Analysis

Documented in GOenrichmentAnalysis

# Function that returns GO terms with best enrichment and highest number of genes.

# So that I don't forget: information about GO categories is contained in the GO.db package

GOenrichmentAnalysis = function(labels, entrezCodes, 
               yeastORFs = NULL,
               organism = "human",
               ontologies = c("BP", "CC", "MF"),
               evidence = "all",
               includeOffspring = TRUE,
               backgroundType = "givenInGO",
               removeDuplicates = TRUE,
               leaveOutLabel = NULL,
               nBestP = 10, pCut = NULL, nBiggest = 0,
               getTermDetails = TRUE,
               verbose = 2, indent = 0 )
{
   warning(immediate. = TRUE,
     spaste("This function is deprecated and will be removed in the near future. \n",
            "We suggest using the replacement function enrichmentAnalysis \n",
            "in R package anRichment, available from the following URL:\n",
            "https://labs.genetics.ucla.edu/horvath/htdocs/CoexpressionNetwork/GeneAnnotation/"));

   sAF = options("stringsAsFactors")
   options(stringsAsFactors = FALSE);
   on.exit(options(stringsAsFactors = sAF[[1]]), TRUE)

   organisms = c("human", "mouse", "rat", "malaria", "yeast", "fly", "bovine", "worm", "canine",
                 "zebrafish", "chicken");
   allEvidence =  c("EXP", "IDA", "IPI", "IMP", "IGI", "IEP", "ISS", "ISO", "ISA",
                    "ISM", "IGC", "IBA", "IBD", "IKR", "IRD", "RCA", "TAS", "NAS", "IC", "ND", "IEA",
                    "NR");
   allOntologies = c("BP", "CC", "MF");

   backgroundTypes = c("allGiven", "allInGO", "givenInGO");

   spaces = indentSpaces(indent);
   orgInd = pmatch(organism, organisms);
   if (is.na(orgInd))
     stop(paste("Unrecognized 'organism' given. Recognized values are ", 
                paste(organisms, collapse = ", ")));

   if (length(evidence)==0)
     stop("At least one valid evidence code must be given in 'evidence'.");
   if (length(ontologies)==0)
     stop("At least one valid ontology code must be given in 'ontology'.");

   if (evidence=="all")
      evidence = allEvidence;

   evidInd = pmatch(evidence, allEvidence);
   if (sum(is.na(evidInd))!=0)
     stop(paste("Unrecognized 'evidence' given. Recognized values are ", 
                paste(allEvidence, collapse = ", ")));

   ontoInd = pmatch(ontologies, allOntologies);
   if (sum(is.na(ontoInd))!=0)
     stop(paste("Unrecognized 'ontologies' given. Recognized values are ", 
                paste(allEvidence, collapse = ", ")));

   backT = pmatch(backgroundType, backgroundTypes);
   if (is.na(backT))
     stop(paste("Unrecognized 'backgroundType' given. Recognized values are ", 
           paste(backgroundTypes,  collapse = ", ")));

   orgCodes = c("Hs", "Mm", "Rn", "Pf", "Sc", "Dm", "Bt", "Ce", "Cf", "Dr", "Gg");
   orgExtensions = c(rep(".eg", 4), ".sgd", rep(".eg", 6));
   reverseMap = c(rep(".egGO2EG", 4), ".sgdGO2ORF", rep(".egGO2EG", 6))

   missingPacks = NULL;
   packageName = paste("org.", orgCodes[orgInd], orgExtensions[orgInd], ".db", sep="");
   if (!eval(parse(text = "require(packageName, character.only = TRUE)")))
     missingPacks = c(missingPacks, packageName);

   if (!eval(parse(text="require(GO.db)")))
     missingPacks = c(missingPacks, "GO.db");

   if (!is.null(missingPacks)) 
     stop(paste("Could not load the requisite package(s)", 
           paste(missingPacks, collapse = ", "), ". Please install the package(s)."))

   if (verbose > 0)
   {
     printFlush(paste(spaces, "GOenrichmentAnalysis: loading annotation data..."));
   }

   if (orgInd==5)
   {
      # Yeast needs special care.
      if (!is.null(yeastORFs))
      {
        entrezCodes = yeastORFs
      } else {
        # Map the entrez IDs to yeast ORFs
        x = eval(parse(text = "org.Sc.sgd:::org.Sc.sgdENTREZID"))
        # x = org.Sc.sgd:::org.Sc.sgdENTREZID
        xx = as.list(x[mapped_genes])
        allORFs = names(xx);
        mappedECs = as.character(sapply(xx, as.character))
        entrez2orf = match(entrezCodes, mappedECs);
        fin = is.finite(entrez2orf);
        newCodes = paste("InvalidCode", c(1:length(entrezCodes)), sep = ".");
        newCodes[fin] = allORFs[entrez2orf[fin]];
        entrezCodes = newCodes;
     }
   } 

   labels = as.matrix(labels);

   nSets = ncol(labels);
   nGivenRaw = nrow(labels);

   if (removeDuplicates)
   {
     # Restrict given entrezCodes such that each code is unique
     ECtab = table(entrezCodes);
     uniqueEC = names(ECtab);
     keepEC = match(uniqueEC, entrezCodes);
     entrezCodes = entrezCodes[keepEC]
     labels = labels[keepEC, , drop = FALSE];
   } else
     keepEC = c(1:nGivenRaw);

   egGO = eval(parse(text = paste(packageName, ":::org.", orgCodes[orgInd], orgExtensions[orgInd], 
                                  "GO", sep = "")));

   if (orgInd==5)
   {
      mapped_genes = as.character(do.call(match.fun("mappedkeys"), list(egGO)));
      encodes2mapped = match(as.character(entrezCodes), mapped_genes);
   } else {
      mapped_genes = as.numeric(as.character(do.call(match.fun("mappedkeys"), list(egGO))));
      encodes2mapped = match(as.numeric(entrezCodes), mapped_genes);
   }
   encMapped = is.finite(encodes2mapped);
   nAllIDsInGO = sum(encMapped);

   mapECodes = entrezCodes[encMapped];
   mapLabels = labels[encMapped, , drop = FALSE];
   nMappedGenes = nrow(mapLabels);

   if (nMappedGenes==0)
     stop(paste("None of the supplied gene identifiers map to the GO database.\n",
                "Please make sure you have specified the correct organism (default is human)."))

   Go2eg = eval(parse(text = paste("AnnotationDbi::as.list(", packageName, ":::org.", orgCodes[orgInd], 
                                           reverseMap[orgInd],")", sep = "")));
   nTerms = length(Go2eg);

   goInfo = as.list(GO.db::GOTERM);
   if (length(goInfo) > 0)
   {
      orgGoNames = names(Go2eg);
      dbGoNames = as.character(sapply(goInfo, GOID));
      dbGoOntologies = as.character(sapply(goInfo, Ontology));
   } else {
      dbGoNames = "";
   }
   goOffSpr = list();
   if (includeOffspring)
   {
     goOffSpr[[1]] = as.list(GOBPOFFSPRING);
     goOffSpr[[2]] = as.list(GOCCOFFSPRING);
     goOffSpr[[3]] = as.list(GOMFOFFSPRING);
   }
   term2info = match(names(Go2eg), names(goInfo));
   termOntologies = dbGoOntologies[term2info];

   if (backT==1)
   {
     nBackgroundGenes = sum(!is.na(entrezCodes));
   } else
   {
     if (backT==2)
     {
       nBackgroundGenes = length(mapped_genes)
     } else
       nBackgroundGenes = nMappedGenes;
   }

   termCodes = vector(mode="list", length = nTerms);
   collectGarbage();
   nExpandLength = 0;
   blockSize = 3000; # For a more efficient concatenating of offspring genes
   nAllInTerm = rep(0, nTerms);
   if (verbose > 0)
   {
      printFlush(paste(spaces, " ..of the", length(entrezCodes), 
                      " Entrez identifiers submitted,", sum(encMapped), 
                      "are mapped in current GO categories."));
      printFlush(paste(spaces, " ..will use", nBackgroundGenes, 
                       "background genes for enrichment calculations."));
      cat(paste(spaces, " ..preparing term lists (this may take a while).."));
      pind = initProgInd();
   }

   for (c in 1:nTerms) if (!is.na(Go2eg[[c]][[1]]))
   {
      te = as.character(names(Go2eg[[c]])); # Term evidence codes
      tc = Go2eg[[c]];
      if (includeOffspring) 
      {
        termOffspring = NULL;
        for (ont in 1:length(goOffSpr))
        {
          term2off = match(names(Go2eg)[c], names(goOffSpr[[ont]]))
          if (!is.na(term2off))
            termOffspring = c(termOffspring, goOffSpr[[ont]][[term2off]]);
        }
        if (length(termOffspring)>0) 
        {
           maxLen = blockSize;
           tex = rep("", maxLen);
           tcx = rep("", maxLen);
           ind = 1;
           len = length(te);
           tex[ ind:len ] = te;
           tcx[ ind:len ] = tc;
           ind = len + 1;
           o2go = match(termOffspring, as.character(names(Go2eg)));
           o2go = o2go[is.finite(o2go)]
           if (length(o2go)>0) for (o in 1:length(o2go)) if (!is.na(Go2eg[[o2go[o]]][[1]]))
           {
             #printFlush(paste("Have offspring for term", c, ": ", names(Go2eg)[c], 
             #           Term(goInfo[[term2info[c]]])));
             newc = Go2eg[[o2go[o]]];
             newe = names(newc);
             newl = length(newe);
             if ((len + newl) > maxLen)
             {
               nExpand = ceiling( (len + newl - maxLen)/blockSize);
               maxLen = maxLen + blockSize * nExpand;
               tex = c(tex, rep("", maxLen - length(tex)));
               tcx = c(tcx, rep("", maxLen - length(tex)));
               nExpandLength = nExpandLength + 1;
             }
             tex[ind:(len + newl)] = newe;
             tcx[ind:(len + newl)] = newc;
             ind = ind + newl;
             len = len + newl;
           }
           te = tex[1:len];
           tc = tcx[1:len];
        }
      }
      use = is.finite(match(te, evidence));
      if (orgInd==5)
      {
         if (backT==2)
         {
            termCodes[[c]] = unique(as.character(tc[use]));
         } else
            termCodes[[c]] = as.character(intersect(tc[use], mapECodes));
      } else {
         if (backT==2)
         {
            termCodes[[c]] = unique(as.character(tc[use]));
         } else
            termCodes[[c]] = as.numeric(as.character(intersect(tc[use], mapECodes)));
      }
      nAllInTerm[c] = length(termCodes[[c]]); 
      if ( (c %%50 ==0) & (verbose > 0)) pind = updateProgInd(c/nTerms, pind);
   }
   if (verbose > 0) 
   {
      pind = updateProgInd(1, pind);
      printFlush("");
   }
   if ((verbose > 5) & (includeOffspring))
      printFlush(paste(spaces, " ..diagnostic for the developer: offspring buffer was expanded",
                       nExpandLength, "times."));

   ftp = function(...) { fisher.test(...)$p.value }

   setResults = list();

   for (set in 1:nSets)
   {
      if (verbose > 0)
        printFlush(paste(spaces, " ..working on label set", set, ".."));
      labelLevels = levels(factor(labels[, set]));
      if (!is.null(leaveOutLabel))
      {
        keep = !(labelLevels %in% as.character(leaveOutLabel));
        if (sum(keep)==0)
          stop("No labels were kept after removing labels that are supposed to be ignored.");
        labelLevels = labelLevels[keep]
      }
      nLabelLevels = length(labelLevels);

      modCodes = list();
      nModCodes = rep(0, nLabelLevels);
      if (backT==1)
      {
        for (ll in 1:nLabelLevels)
        {
           modCodes[[ll]] = entrezCodes[labels[, set]==labelLevels[ll]];
           nModCodes[ll] = length(modCodes[[ll]]);
        }
      } else 
      {
        for (ll in 1:nLabelLevels)
        {
           modCodes[[ll]] = mapECodes[mapLabels[, set]==labelLevels[ll]];
           nModCodes[ll] = length(modCodes[[ll]]);
        }
      }

      countsInTerm = matrix(0, nLabelLevels, nTerms);
      enrichment = matrix(1, nLabelLevels, nTerms);

      for (ll in 1:nLabelLevels)
        countsInTerm[ll, ] = sapply(lapply(termCodes, intersect, modCodes[[ll]]), length)

      nAllInTermMat = matrix(nAllInTerm, nLabelLevels, nTerms, byrow = TRUE);
      nModCodesMat = matrix(nModCodes, nLabelLevels, nTerms);
      tabArr = array(c(countsInTerm, nAllInTermMat - countsInTerm, 
                       nModCodesMat - countsInTerm,
                       nBackgroundGenes - nModCodesMat - nAllInTermMat + countsInTerm),
                     dim = c(nLabelLevels * nTerms, 2, 2));

      if (verbose > 0)
         printFlush(paste(spaces, "   ..calculating enrichments (this may also take a while).."));
      calculate = c(countsInTerm) > 0;
      enrichment[calculate] = apply(tabArr[calculate, , ], 1, ftp, alternative = "g");
    
      dimnames(enrichment) = list (labelLevels, names(Go2eg));
      dimnames(countsInTerm) = list (labelLevels, names(Go2eg));

      bestPTerms = list();
      modSizes = table(labels[ !(labels[, set] %in% leaveOutLabel), set]);

      if (!is.null(pCut) || nBestP > 0)
      {
         printFlush(paste(spaces, "   ..putting together terms with highest enrichment significance.."));
         nConsideredOntologies = length(ontologies)+1;
         for (ont in 1:nConsideredOntologies)
         {
            if (ont==nConsideredOntologies)
            {
               ontTerms = is.finite(match(termOntologies, ontologies))
               bestPTerms[[ont]] = list(ontology = ontologies);
               names(bestPTerms)[ont] = paste(ontologies, collapse = ", ");
            } else {
               ontTerms = termOntologies==ontologies[ont];
               bestPTerms[[ont]] = list(ontology = ontologies[ont]);
               names(bestPTerms)[ont] = ontologies[ont];
            }
            bestPTerms[[ont]]$enrichment = NULL;
            bestPTerms[[ont]]$termInfo = list();
            nOntTerms = sum(ontTerms)
            ontEnr = enrichment[, ontTerms, drop = FALSE];
            order = apply(ontEnr, 1, order);
            for (ll in 1:nLabelLevels)
            {
              if (!is.null(pCut))
              {
                 reportTerms = c(1:nTerms)[ontTerms][ontEnr[ll, ] < pCut];
                 reportTerms = reportTerms[order(ontEnr[ll, ][reportTerms])];
              } else
                 reportTerms =  c(1:nTerms)[ontTerms][order[1:nBestP, ll]];
              nRepTerms = length(reportTerms);
              enrTab = data.frame(module = rep(labelLevels[ll], nRepTerms), 
                                  modSize = rep(modSizes[ll], nRepTerms),
                                  bkgrModSize = rep(nModCodes[ll], nRepTerms), 
                                  rank = seq(length.out = nRepTerms),
                                  enrichmentP = enrichment[ll, reportTerms],
                                  BonferoniP = pmin(rep(1, nRepTerms), 
                                                    enrichment[ll, reportTerms] * nOntTerms),
                                  nModGenesInTerm = countsInTerm[ll, reportTerms],
                                  fracOfBkgrModSize = countsInTerm[ll, reportTerms]/nModCodes[ll],
                                  fracOfBkgrTermSize = 
                                       countsInTerm[ll, reportTerms]/nAllInTerm[reportTerms],
                                  bkgrTermSize = nAllInTerm[reportTerms],
                                  termID = names(Go2eg)[reportTerms],
                                  termOntology =  rep("NA", nRepTerms),
                                  termName = rep("NA", nRepTerms),
                                  termDefinition =  rep("NA", nRepTerms))
              bestPTerms[[ont]]$forModule[[ll]] = list();
              for (rci in seq(length.out = nRepTerms))
              {
                 term = reportTerms[rci];
                 termID = names(Go2eg)[term];
                 dbind = match(termID, dbGoNames);
                 if (is.finite(dbind))
                 {
                   enrTab$termName[rci] = eval(parse(text = "AnnotationDbi:::Term(goInfo[[dbind]])"));
                   enrTab$termDefinition[rci] = 
                         eval(parse(text = "AnnotationDbi:::Definition(goInfo[[dbind]])"));
                   enrTab$termOntology[rci] = 
                         eval(parse(text = "AnnotationDbi:::Ontology(goInfo[[dbind]])"));
                 } 
                 if (getTermDetails)
                 {
                   geneCodes = intersect(modCodes[[ll]], termCodes[[term]])
                   bestPTerms[[ont]]$forModule[[ll]][[rci]] = list(termID = termID,
                           termName = enrTab$termName[rci],
                           enrichmentP = enrTab$enrichmentP[rci],
                           termDefinition = enrTab$termDefinition[rci],
                           termOntology = enrTab$termOntology[rci],
                           geneCodes = geneCodes,
                           genePositions = keepEC[match(geneCodes, entrezCodes)]);
                }
              }
              if (ll==1) 
              {
                 bestPTerms[[ont]]$enrichment = enrTab;
              } else 
                 bestPTerms[[ont]]$enrichment = rbind(bestPTerms[[ont]]$enrichment, enrTab);
           }
         }
      }
           
      biggestTerms = list();
   
      if (nBiggest > 0)
      {
         printFlush(paste(spaces, "   ..putting together terms with largest number of genes in modules.."));
         nConsideredOntologies = length(ontologies)+1;
         for (ont in 1:nConsideredOntologies)
         {
            if (ont==nConsideredOntologies)
            {
               ontTerms = is.finite(match(termOntologies, ontologies))
               biggestTerms[[ont]] = list(ontology = ontologies);
               names(biggestTerms)[ont] = paste(ontologies, collapse = ", ");
            } else {
               ontTerms = termOntologies==ontologies[ont];
               biggestTerms[[ont]] = list(ontology = ontologies[ont]);
               names(biggestTerms)[ont] = ontologies[ont];
            }
            biggestTerms[[ont]]$enrichment = NULL;
            biggestTerms[[ont]]$termInfo = list();
            nOntTerms = sum(ontTerms)
            ontNGenes = countsInTerm[, ontTerms, drop = FALSE];
            order = apply(-ontNGenes, 1, order);
            for (ll in 1:nLabelLevels)
            {
              reportTerms = c(1:nTerms)[ontTerms][order[1:nBiggest, ll]];
              nRepTerms = length(reportTerms);
              enrTab = data.frame(module = rep(labelLevels[ll], nRepTerms),
                                  modSize = rep(modSizes[ll], nRepTerms),
                                  bkgrModSize = rep(nModCodes[ll], nRepTerms), 
                                  rank = seq(length.out = nRepTerms),
                                  enrichmentP = enrichment[ll, reportTerms],
                                  BonferoniP = pmin(rep(1, nRepTerms), 
                                                    enrichment[ll, reportTerms] * nOntTerms),
                                  nModGenesInTerm = countsInTerm[ll, reportTerms],
                                  fracOfModSize = countsInTerm[ll, reportTerms]/nModCodes[ll],
                                  fracOfBkgrTermSize = 
                                       countsInTerm[ll, reportTerms]/nAllInTerm[reportTerms],
                                  bkgrTermSize = nAllInTerm[reportTerms],
                                  termID = names(Go2eg)[reportTerms],
                                  termOntology =  rep("NA", nRepTerms),
                                  termName = rep("NA", nRepTerms),
                                  termDefinition =  rep("NA", nRepTerms))
              biggestTerms[[ont]]$forModule[[ll]] = list();
              for (rci in seq(length.out = nRepTerms))
              {
                 term = reportTerms[rci];
                 termID = names(Go2eg)[term];
                 dbind = match(termID, dbGoNames);
                 if (is.finite(dbind))
                 {
                   enrTab$termName[rci] = eval(parse(text="AnnotationDbi:::Term(goInfo[[dbind]])"));
                   enrTab$termDefinition[rci] =
                       eval(parse(text="AnnotationDbi:::Definition(goInfo[[dbind]])"));
                   enrTab$termOntology[rci] = eval(parse(text="AnnotationDbi:::Ontology(goInfo[[dbind]])"));
                 }
                 if (getTermDetails)
                 {
                   geneCodes = intersect(modCodes[[ll]], termCodes[[term]])
                   biggestTerms[[ont]]$forModule[[ll]][[rci]] = list(termID = termID,
                           termName = enrTab$termName[rci],
                           enrichmentP = enrTab$enrichmentP[rci],
                           termDefinition = enrTab$termDefinition[rci],
                           termOntology = enrTab$termOntology[rci],
                           geneCodes = geneCodes,
                           genePositions = keepEC[match(geneCodes, entrezCodes)]);
                }
              }
              if (ll==1)
              {
                 biggestTerms[[ont]]$enrichment = enrTab;
              } else
                 biggestTerms[[ont]]$enrichment = rbind(biggestTerms[[ont]]$enrichment, enrTab);
           }
         }
      }
      setResults[[set]] = list(countsInTerm = countsInTerm, 
                               enrichmentP = enrichment,
                               bestPTerms = bestPTerms,
                               biggestTerms = biggestTerms);
   }
   
   inGO = rep(FALSE, nGivenRaw);
   inGO[keepEC] = encMapped;
   kept = rep(FALSE, nGivenRaw);
   kept[keepEC] = TRUE;
   if (nSets==1)
   {
      list(keptForAnalysis = kept,
           inGO = inGO,
           countsInTerms = setResults[[1]]$countsInTerm,
           enrichmentP = setResults[[1]]$enrichmentP,
           bestPTerms = setResults[[1]]$bestPTerms,
           biggestTerms = setResults[[1]]$biggestTerms);
   } else {
      list(keptForAnalysis = kept,
           inGO = inGO,
           setResults = setResults);
   }
}