text: user285: ## ------------------------------------------------------------------------
  library(limma)
  library(statmod)
  library(GEOquery)
  #library(GEOmetadb)
  require(Biobase)
  library(ggplot2)
  library(rio)
  library(affy)
  library(affyPLM)
  library(reshape)
  library(stringr)
  library(cowplot)
  #library(makecdfenv)
  #library(aroma.affymetrix)
  #library(RTN)
  library(gProfileR)
  library(corpcor)
  library(viper)
  library(rio)
  library(dplyr)
  library(igraph)
  library(RColorBrewer)
  library(reshape2)
  
  
  ## ------------------------------------------------------------------------
  gse89988 <- getGEO(filename = "GSE89988_family.soft.gz", GSEMatrix = TRUE)
  
  
  ## ------------------------------------------------------------------------
  # define samples
  samples_nt <- c('GSM2394292', 'GSM2394303')
  samples_pam <- c('GSM2394295', 'GSM2394306')
  samples_lps <- c('GSM2394296', 'GSM2394307')
  
  # make table with samples metadata
  samples <- c(samples_nt, samples_pam, samples_lps)
  metadata <- as.data.frame(samples)
  colnames(metadata) <- 'sample_id'
  
  # add group information both with the group acronym as a dummy variable -- which
  # will be useful for the DEA; also add color for differentiate groups in plots
  metadata$group <- c(rep(0, 2), rep(1, 2), rep(2, 2))
  metadata$group_name <- c(rep('NT', 2), rep('PAM', 2), rep('LPS', 2))
  metadata$group_color <- c(rep('gray', 2), rep('tomato', 2), rep('dodgerblue', 2))
  
  # get positions in GSE object
  match(samples, names(GSMList(gse89988)))
  
  # with the following I confirm that the samples correspond to the expected samples
  # and get that expression values are already normalized and log2-transformed
  GSMList(gse89988)[[49]]
  GSMList(gse89988)[[60]]
  GSMList(gse89988)[[52]]
  GSMList(gse89988)[[63]]
  GSMList(gse89988)[[53]]
  GSMList(gse89988)[[64]]
  
  
  ## ------------------------------------------------------------------------
  #First, we need to make sure that all of the GSMs are from the same platform
  gsmplatforms <- lapply(GSMList(gse89988),function(x) {Meta(x)$platform_id})
  head(gsmplatforms)
  # from the output of the command above we see that indeed all GSMs derive from
  # the same platform "GPL17543"
  
  # We can filter the original GSMList to include only those GSMs with the GPL96 platform and use this list for further processing
  gsmlist = Filter(function(gsm) {Meta(gsm)$platform_id=='GPL17543'},GSMList(gse89988))
  length(gsmlist)
  
  # get the probeset ordering
  probesets <- Table(GPLList(gse89988)[[1]])$ID
  
  # make the data matrix from the VALUE columns from each GSM
  # being careful to match the order of the probesets in the platform
  # with those in the GSMs
  data.matrix <- do.call('cbind', lapply(gsmlist,function(x) 
                                        {tab <- Table(x)
                                        mymatch <- match(probesets,tab$ID_REF)
                                        return(tab$VALUE[mymatch])
                                       }))
  data.matrix <- apply(data.matrix, 2, function(x) {as.numeric(as.character(x))})
  
  # go through the necessary steps to make a compliant ExpressionSet
  # and filter for the samples to be analyzed
  rownames(data.matrix) <- probesets
  colnames(data.matrix) <- names(gsmlist)
  data.matrix <- data.matrix[, samples]
  pdata <- data.frame(samples = samples)
  rownames(pdata) <- samples
  pheno <- as(pdata,"AnnotatedDataFrame")
  eset <- new('ExpressionSet', exprs = data.matrix, phenoData = pheno)
  
  
  ## ------------------------------------------------------------------------
  load_gene_annotation <- function() {
  
    ifile <- 'GPL17543-38936.txt'
    df <- read.delim(ifile, comment.char = '#')
  
    return(df)
  
  }
  
  gene_annotation <- load_gene_annotation()
  
  
  ## ------------------------------------------------------------------------
  apply_dea <- function(eset, fc_cutoff = 2, adj_p_cutoff = 0.05) {
  
    # output directory
    ODIR <- 'results/step1/'
    dir.create(ODIR, showWarnings = FALSE, recursive = TRUE)
  
    # make design matrix
    design <- model.matrix(~ 0 + factor(c(0, 0, 1, 1, 2, 2)))
    colnames(design) <- c("group0", "group1", "group2")
  
    # plot with the normalized array weights
    aWeights <- arrayWeights(eset, design = design)
    ofile <- paste0(ODIR, 'array_weights_gse89988_samples_used_in_dea.pdf')
    pdf(ofile, width = 10, height = 5)
    dotchart(aWeights, labels = samples,
             groups = as.factor(metadata$group_name),
             pch = 19,
             color = metadata$group_color,
             main = "Array Weights")  
    dev.off()
  
    # linear model fit
    fit <- lmFit(eset, design)
    contrast.matrix <- makeContrasts(group1-group0, group2-group0,
                                     levels = design)
    fit2 <- contrasts.fit(fit, contrast.matrix)
    fit2 <- eBayes(fit2)
  
  
    # save 
    ODIR='deliverables/'
    dir.create(ODIR, showWarnings = FALSE)
    ofile=paste0(ODIR, 'fit_for_figures_7_and_8.rds')
    saveRDS(fit2, ofile)
  
  
    # extract results, annotate and visualize DEA results
  
    extract_annotate_and_viz_dea_results <- function(coef_index, ref, treat) {
  
      # Extract DEA results for a given comparison
      df <- topTable(fit2, coef = coef_index, adjust = "BH", number = Inf)
      df$reference <- ref
      df$treatment <- treat
  
      # remove genes with missing values
      df <- df[complete.cases(df), ]
  
      # annotate with gene information -- gene annotation has many fields which
      # are not relevant here so I will make a reduced gene information version
      fields_included <- c('ID', 'Species', 'Source', 'Entrez_Gene_ID', 'Symbol')
      df <- merge(df, gene_annotation[fields_included], by.x = 0, by.y = 'ID',
                  sort = FALSE)
      colnames(df)[1] <- 'probe_id'
  
      # prepare data for visualization through volcano plot
      cond1 <- df$adj.P.Val < adj_p_cutoff
      cond2d <- df$logFC < -log2(fc_cutoff)
      cond2u <- df$logFC > log2(fc_cutoff)
      log10p <- -log10(df$adj.P.Val)
      down <- nrow(df[cond1 & cond2d, ])
      up <- nrow(df[cond1 & cond2u, ])
      df$category <- ifelse(cond1 & cond2d, 'down-regulated',
                            ifelse(cond1 & cond2u, 'up-regulated', 'non-regulated'))
      abs_max_logfc <- max(abs(df$logFC))
  
      # volcano plot
      g <- ggplot(df, aes(x = logFC, y = -log10(adj.P.Val), colour = category)) +
            geom_point(alpha = 0.4, size = 1.75) +
            scale_color_manual(values = c('dodgerblue', 'gray', 'tomato')) +
            xlim(-abs_max_logfc, abs_max_logfc) +
            ylim(c(0, max(-log10(df$adj.P.Val)))) +
            xlab("log2 fold change") +
            ylab("-log10(adjusted p-value)") +
            annotate("text", x = -4, y = 3, size =3, 
                     label = paste0(down, "\ndown-regulated DEGs"),
                     color = 'dodgerblue') +
            annotate("text", x = 4, y = 3, size = 3,
                     label = paste0(up, "\nup-regulated DEGs"),
                     color = 'tomato') +
            annotate("text", x = -abs_max_logfc, y = 0, size = 2,
                     label = paste0("DEGs defined with:\nAdj. p-value < ",
                                    adj_p_cutoff,
                                    '\nFold change > ',
                                    fc_cutoff),
                     color = 'black', hjust = 0) +
            ggtitle(paste0(ref, ' vs. ', treat)) +
            theme(plot.title = element_text(hjust = 0.5))
  
      # save plot to file
      ofile <- paste0(ODIR, 'volcano_plot_', tolower(ref), '_vs_', tolower(treat),
                      '.pdf')
      ggsave(ofile, width = 7, height = 5)
  
      return(df)
  
    }
  
    res <- rbind(extract_annotate_and_viz_dea_results(1, 'NT', 'PAM'),
                 extract_annotate_and_viz_dea_results(2, 'NT', 'LPS'))
  
  
    # export table with results
    ofile <- paste0(ODIR, 'differential_expression_analysis_results_combined.xlsx')
    export(res, ofile)
  
  }
  
  
  ## ------------------------------------------------------------------------
  apply_dea(eset)
  
  After all these i have to run this 
  gsca <- pathwayAnalysis(fit2, pvaluecutoff=0.05),where fit2 is made after running apply_dea(eset) 
  As I run the above code i get this 
  
  Error in pathwayAnalysis(fit2, pvaluecutoff = 0.05) :    could not
    find function "pathwayAnalysis"
  
  I would share the .rdata file for reproduciblty
  rdata
  The gse & the annotation files
  I would be really glad if i get to know what am i doing wrong.
  
  user57: It seems like you are missing function "pathwayAnalysis" presumably from package BLMA, i.e. library('BLMA') should resolve your error once you install that package. 
  -- edit 1 --
  Ha, stupid me. Package BLMA contains file "pathwayAnalysis.R", not the function. A package that contains a FUNCTION called pathwayAnalysis is IntClust. You should try that one.
  -- edit 2 --
  Alright, stupid me again, now I have not checked that IntClust has capitalized PathwayAnalysis, you need pathwayAnalysis. I found one more function matching the name form package metaX. However, it could be that it's function defined somewhere in their code. You can try to search the source code for example using grep or you try to use google to find packages that contain the function of this name (this is exactly what I was doing, many you should be a bit more careful than I was).
  
original_text: 285: <pre><code>## ------------------------------------------------------------------------
  library(limma)
  library(statmod)
  library(GEOquery)
  #library(GEOmetadb)
  require(Biobase)
  library(ggplot2)
  library(rio)
  library(affy)
  library(affyPLM)
  library(reshape)
  library(stringr)
  library(cowplot)
  #library(makecdfenv)
  #library(aroma.affymetrix)
  #library(RTN)
  library(gProfileR)
  library(corpcor)
  library(viper)
  library(rio)
  library(dplyr)
  library(igraph)
  library(RColorBrewer)
  library(reshape2)
  
  
  ## ------------------------------------------------------------------------
  gse89988 &lt;- getGEO(filename = "GSE89988_family.soft.gz", GSEMatrix = TRUE)
  
  
  ## ------------------------------------------------------------------------
  # define samples
  samples_nt &lt;- c('GSM2394292', 'GSM2394303')
  samples_pam &lt;- c('GSM2394295', 'GSM2394306')
  samples_lps &lt;- c('GSM2394296', 'GSM2394307')
  
  # make table with samples metadata
  samples &lt;- c(samples_nt, samples_pam, samples_lps)
  metadata &lt;- as.data.frame(samples)
  colnames(metadata) &lt;- 'sample_id'
  
  # add group information both with the group acronym as a dummy variable -- which
  # will be useful for the DEA; also add color for differentiate groups in plots
  metadata<span class="math-container">$group &lt;- c(rep(0, 2), rep(1, 2), rep(2, 2))
  metadata$</span>group_name &lt;- c(rep('NT', 2), rep('PAM', 2), rep('LPS', 2))
  metadata$group_color &lt;- c(rep('gray', 2), rep('tomato', 2), rep('dodgerblue', 2))
  
  # get positions in GSE object
  match(samples, names(GSMList(gse89988)))
  
  # with the following I confirm that the samples correspond to the expected samples
  # and get that expression values are already normalized and log2-transformed
  GSMList(gse89988)[[49]]
  GSMList(gse89988)[[60]]
  GSMList(gse89988)[[52]]
  GSMList(gse89988)[[63]]
  GSMList(gse89988)[[53]]
  GSMList(gse89988)[[64]]
  
  
  ## ------------------------------------------------------------------------
  #First, we need to make sure that all of the GSMs are from the same platform
  gsmplatforms &lt;- lapply(GSMList(gse89988),function(x) {Meta(x)$platform_id})
  head(gsmplatforms)
  # from the output of the command above we see that indeed all GSMs derive from
  # the same platform "GPL17543"
  
  # We can filter the original GSMList to include only those GSMs with the GPL96 platform and use this list for further processing
  gsmlist = Filter(function(gsm) {Meta(gsm)$platform_id=='GPL17543'},GSMList(gse89988))
  length(gsmlist)
  
  # get the probeset ordering
  probesets &lt;- Table(GPLList(gse89988)[[1]])$ID
  
  # make the data matrix from the VALUE columns from each GSM
  # being careful to match the order of the probesets in the platform
  # with those in the GSMs
  data.matrix &lt;- do.call('cbind', lapply(gsmlist,function(x) 
                                        {tab &lt;- Table(x)
                                        mymatch &lt;- match(probesets,tab<span class="math-container">$ID_REF)
                                        return(tab$</span>VALUE[mymatch])
                                       }))
  data.matrix &lt;- apply(data.matrix, 2, function(x) {as.numeric(as.character(x))})
  
  # go through the necessary steps to make a compliant ExpressionSet
  # and filter for the samples to be analyzed
  rownames(data.matrix) &lt;- probesets
  colnames(data.matrix) &lt;- names(gsmlist)
  data.matrix &lt;- data.matrix[, samples]
  pdata &lt;- data.frame(samples = samples)
  rownames(pdata) &lt;- samples
  pheno &lt;- as(pdata,"AnnotatedDataFrame")
  eset &lt;- new('ExpressionSet', exprs = data.matrix, phenoData = pheno)
  
  
  ## ------------------------------------------------------------------------
  load_gene_annotation &lt;- function() {
  
    ifile &lt;- 'GPL17543-38936.txt'
    df &lt;- read.delim(ifile, comment.char = '#')
  
    return(df)
  
  }
  
  gene_annotation &lt;- load_gene_annotation()
  
  
  ## ------------------------------------------------------------------------
  apply_dea &lt;- function(eset, fc_cutoff = 2, adj_p_cutoff = 0.05) {
  
    # output directory
    ODIR &lt;- 'results/step1/'
    dir.create(ODIR, showWarnings = FALSE, recursive = TRUE)
  
    # make design matrix
    design &lt;- model.matrix(~ 0 + factor(c(0, 0, 1, 1, 2, 2)))
    colnames(design) &lt;- c("group0", "group1", "group2")
  
    # plot with the normalized array weights
    aWeights &lt;- arrayWeights(eset, design = design)
    ofile &lt;- paste0(ODIR, 'array_weights_gse89988_samples_used_in_dea.pdf')
    pdf(ofile, width = 10, height = 5)
    dotchart(aWeights, labels = samples,
             groups = as.factor(metadata<span class="math-container">$group_name),
             pch = 19,
             color = metadata$</span>group_color,
             main = "Array Weights")  
    dev.off()
  
    # linear model fit
    fit &lt;- lmFit(eset, design)
    contrast.matrix &lt;- makeContrasts(group1-group0, group2-group0,
                                     levels = design)
    fit2 &lt;- contrasts.fit(fit, contrast.matrix)
    fit2 &lt;- eBayes(fit2)
  
  
    # save 
    ODIR='deliverables/'
    dir.create(ODIR, showWarnings = FALSE)
    ofile=paste0(ODIR, 'fit_for_figures_7_and_8.rds')
    saveRDS(fit2, ofile)
  
  
    # extract results, annotate and visualize DEA results
  
    extract_annotate_and_viz_dea_results &lt;- function(coef_index, ref, treat) {
  
      # Extract DEA results for a given comparison
      df &lt;- topTable(fit2, coef = coef_index, adjust = "BH", number = Inf)
      df<span class="math-container">$reference &lt;- ref
      df$</span>treatment &lt;- treat
  
      # remove genes with missing values
      df &lt;- df[complete.cases(df), ]
  
      # annotate with gene information -- gene annotation has many fields which
      # are not relevant here so I will make a reduced gene information version
      fields_included &lt;- c('ID', 'Species', 'Source', 'Entrez_Gene_ID', 'Symbol')
      df &lt;- merge(df, gene_annotation[fields_included], by.x = 0, by.y = 'ID',
                  sort = FALSE)
      colnames(df)[1] &lt;- 'probe_id'
  
      # prepare data for visualization through volcano plot
      cond1 &lt;- df<span class="math-container">$adj.P.Val &lt; adj_p_cutoff
      cond2d &lt;- df$</span>logFC &lt; -log2(fc_cutoff)
      cond2u &lt;- df<span class="math-container">$logFC &gt; log2(fc_cutoff)
      log10p &lt;- -log10(df$</span>adj.P.Val)
      down &lt;- nrow(df[cond1 &amp; cond2d, ])
      up &lt;- nrow(df[cond1 &amp; cond2u, ])
      df<span class="math-container">$category &lt;- ifelse(cond1 &amp; cond2d, 'down-regulated',
                            ifelse(cond1 &amp; cond2u, 'up-regulated', 'non-regulated'))
      abs_max_logfc &lt;- max(abs(df$</span>logFC))
  
      # volcano plot
      g &lt;- ggplot(df, aes(x = logFC, y = -log10(adj.P.Val), colour = category)) +
            geom_point(alpha = 0.4, size = 1.75) +
            scale_color_manual(values = c('dodgerblue', 'gray', 'tomato')) +
            xlim(-abs_max_logfc, abs_max_logfc) +
            ylim(c(0, max(-log10(df$adj.P.Val)))) +
            xlab("log2 fold change") +
            ylab("-log10(adjusted p-value)") +
            annotate("text", x = -4, y = 3, size =3, 
                     label = paste0(down, "\ndown-regulated DEGs"),
                     color = 'dodgerblue') +
            annotate("text", x = 4, y = 3, size = 3,
                     label = paste0(up, "\nup-regulated DEGs"),
                     color = 'tomato') +
            annotate("text", x = -abs_max_logfc, y = 0, size = 2,
                     label = paste0("DEGs defined with:\nAdj. p-value &lt; ",
                                    adj_p_cutoff,
                                    '\nFold change &gt; ',
                                    fc_cutoff),
                     color = 'black', hjust = 0) +
            ggtitle(paste0(ref, ' vs. ', treat)) +
            theme(plot.title = element_text(hjust = 0.5))
  
      # save plot to file
      ofile &lt;- paste0(ODIR, 'volcano_plot_', tolower(ref), '_vs_', tolower(treat),
                      '.pdf')
      ggsave(ofile, width = 7, height = 5)
  
      return(df)
  
    }
  
    res &lt;- rbind(extract_annotate_and_viz_dea_results(1, 'NT', 'PAM'),
                 extract_annotate_and_viz_dea_results(2, 'NT', 'LPS'))
  
  
    # export table with results
    ofile &lt;- paste0(ODIR, 'differential_expression_analysis_results_combined.xlsx')
    export(res, ofile)
  
  }
  
  
  ## ------------------------------------------------------------------------
  apply_dea(eset)
  </code></pre>
  
  <p>After all these i have to run this </p>
  
  <p><code>gsca &lt;- pathwayAnalysis(fit2, pvaluecutoff=0.05)</code>,where fit2 is made after running <code>apply_dea(eset)</code> </p>
  
  <p>As I run the above code i get this </p>
  
  <blockquote>
    <p>Error in pathwayAnalysis(fit2, pvaluecutoff = 0.05) :    could not
    find function "pathwayAnalysis"</p>
  </blockquote>
  
  <p>I would share the .rdata file for reproduciblty<br>
  <a href="https://drive.google.com/file/d/1TGHa0TKV9J6glq7Q5JhoFMUM94IXa0eN/view?usp=sharing" rel="nofollow noreferrer">rdata</a></p>
  
  <p>The gse &amp; the annotation <a href="https://drive.google.com/open?id=1QkxCaIlzNwa4BoBZ_KEhaVYocUrJNbuD" rel="nofollow noreferrer">files</a></p>
  
  <p>I would be really glad if i get to know what am i doing wrong.</p>
   57: <p>It seems like you are missing function "<a href="https://rdrr.io/bioc/BLMA/src/R/PathwayAnalysis.R" rel="nofollow noreferrer">pathwayAnalysis</a>" presumably from package BLMA, i.e. <code>library('BLMA')</code> should resolve your error once you install that package. </p>
  
  <p>-- edit 1 --</p>
  
  <p>Ha, stupid me. Package <code>BLMA</code> contains file "pathwayAnalysis.R", not the function. A package that contains a FUNCTION called <code>pathwayAnalysis</code> is <a href="https://www.rdocumentation.org/packages/IntClust/versions/0.1.0/topics/PathwayAnalysis" rel="nofollow noreferrer">IntClust</a>. You should try that one.</p>
  
  <p>-- edit 2 --</p>
  
  <p>Alright, stupid me again, now I have not checked that IntClust has capitalized <code>PathwayAnalysis</code>, you need <code>pathwayAnalysis</code>. I found <a href="https://www.rdocumentation.org/packages/metaX/versions/1.4.0/topics/pathwayAnalysis" rel="nofollow noreferrer">one more function</a> matching the name form package <code>metaX</code>. However, it could be that it's function defined somewhere in their code. You can try to <a href="https://readwrite.com/2010/11/10/how-to-search-your-source-with/" rel="nofollow noreferrer">search the source code for example using grep</a> or you try to use google to find packages that contain the function of this name (this is exactly what I was doing, many you should be a bit more careful than I was).</p>
  