2.1 Workflow package installation from Bioconductor

You can install the package via the BiocManager.

if (!require("BiocManager"))
    install.packages("BiocManager")
BiocManager::install("maEndToEnd", version = "devel")

Currently, the workflow is available in the development version of Bioconductor (3.8), which will become the release version in October 2018.

For details on how to use this version of Bioconductor see:

http://bioconductor.org/developers/how-to/useDevel/

2.2 Workflow package installation from Github

#In order to download the package from GitHub, we need the "install_github" 
#function from the "remotes" package. We download the latest developer 
#version of "remotes" from GitHub with the devtool::install_github
#function; note that this is necessary as the current "remotes" version on
#CRAN doesn't allow us to correctly download the "maEndToEnd" package: 

install.packages("devtools")
library(devtools)

devtools::install_github("r-lib/remotes")
library(remotes)
packageVersion("remotes") # has to be 1.1.1.9000 or later

remotes::install_github("b-klaus/maEndToEnd", ref="master")

2.3 Workflow package import

Once the workflow package has been successfully installed, we can use a call to library() in order to load it. This will also load all the other packages neccessary to run the workflow.

suppressPackageStartupMessages({library("maEndToEnd")})

   
   groupGOTerms:    GOBPTerm, GOMFTerm, GOCCTerm environments built.

2.4 List of packages required for the workflow

Below, you find a list of packages that are required by the workflow. Some Helper/Styling packages have been commented here, as they are not strictly neccesary to execute the workflow.

#General Bioconductor packages
    library(Biobase)
    library(oligoClasses)
     
#Annotation and data import packages
    library(ArrayExpress)
    library(pd.hugene.1.0.st.v1)
    library(hugene10sttranscriptcluster.db)
     
#Quality control and pre-processing packages
    library(oligo)
    library(arrayQualityMetrics)
     
#Analysis and statistics packages
    library(limma)
    library(topGO)
    library(ReactomePA)
    library(clusterProfiler)
     
#Plotting and color options packages
    library(gplots)
    library(ggplot2)
    library(geneplotter)
    library(RColorBrewer)
    library(pheatmap)
     
#Formatting/documentation packages
   #library(rmarkdown)
   #library(BiocStyle)
    library(dplyr)
    library(tidyr)

#Helpers:
    library(stringr)
    library(matrixStats)
    library(genefilter)
    library(openxlsx)
   #library(devtools)

4.1 Information stored in ArrayExpress

Each dataset at ArrayExpress is stored according to the MAGE-TAB (MicroArray Gene Expression Tabular) specifications as a collection of tables bundled with the raw data. The MAGE-TAB format specifies up to five different types of files:

• Investigation Description Format (IDF)
• Array Design Format (ADF)
• Sample and Data Relationship Format (SDRF)
• raw data files
• processed data files

Other than the raw data files, the IDF and the SDRF file are important for us. The IDF file contains top level information about the experiment including title, description, submitter contact details and protocols. The SDRF file contains essential information on the experimental samples, e.g. the experimental group(s) they belong to.

Before we move on to the actual raw data import, we will briefly introduce the ExpressionSet class contained in the Biobase package. It is commonly used to store microarray data in Bioconductor.

4.2 Bioconductor ExpressionSets

Genomic data can be very complex, usually consisting of a number of different components, e.g. information on the experimental samples, annotation of genomic features measured as well as the experimental data itself. In Bioconductor, the approach is taken that these components should be stored in a single structure to easily manage the data.

The package Biobase contains standardized data structures to represent genomic data. The ExpressionSet class is designed to combine several different sources of information (i.e. as contained in the various MAGE-TAB files) into a single convenient structure. An ExpressionSet can be manipulated (e.g., subsetted, copied), and is the input to or output of many Bioconductor functions.

The data in an ExpressionSet consist of:

• assayData: Expression data from microarray experiments with microarray probes in rows and sample identifiers in columns
• metaData
  • phenoData: A description of the samples in the experiment with sample identifiers in rows and description elements in columns; holds the content of the SDRF file
  • featureData: metadata about the features on the chip
    or technology used for the experiment with same rows as assayData by default and freely assignable columns
  • further annotations for the features, for example gene annotations from biomedical databases (annotation).
• experimentData: A flexible structure to describe the experiment.

The ExpressionSet class coordinates all of these data, so that one does not have to worry about the details. However, one should keep in mind that the rownames of the phenoData have to match the column names of the assay data, while the row names of the assay data have to match the row names of the featureData. This is illustrated in Figure 1.

Figure 1: Structure of Bioconductor’s ExpressionSet class

You can use the functions pData and fData to extract the sample and feature annotation, respectively, from an ExpressionSet. The function exprs will return the expression data itself as a matrix.

7.1 Background adjustment

After the initial import and quality assessment, the next step in processing of microarray data is background adjustment. This is essential because a proportion of the measured probe intensities are due to non-specific hybridization and the noise in the optical detection system. Therefore, observed intensities need to be adjusted to give accurate measurements of specific hybridization.

7.2 Across-array normalization (calibration)

Normalization across arrays is needed in order to be able to compare measurements from different array hybridizations due to many obscuring sources of variation. These include different efficiencies of reverse transcription, labeling or hybridization reactions, physical problems with the arrays, reagent batch effects, and laboratory conditions.

7.3 Summarization

After normalization, summarization is necessary to be done because on the Affymetrix platform, transcripts are represented by multiple probes, that is multiple locations on the array. For each gene, the background-adjusted and normalized intensities of all probes need to be summarized into one quantity that estimates an amount proportional to the amount of RNA transcript.

After the summarization step, the summarized data can be annotated with various information, e.g. gene symbols and ENSEMBL gene identifiers. There is an annotation database available from Bioconductor for our platform, namely the package hugene10sttranscriptcluster.db.

You can view its content like this:

head(ls("package:hugene10sttranscriptcluster.db"))

   [1] "hugene10sttranscriptcluster"           
   [2] "hugene10sttranscriptcluster.db"        
   [3] "hugene10sttranscriptclusterACCNUM"     
   [4] "hugene10sttranscriptclusterALIAS2PROBE"
   [5] "hugene10sttranscriptclusterCHR"        
   [6] "hugene10sttranscriptclusterCHRLENGTHS"

Additional information is available from the reference manual of the package. Essentially, the package provides a mapping from the transcript cluster identifiers to the various annotation data.

7.4 Old and new “probesets” of Affymetrix microarrays

Traditionally, Affymetrix arrays (the so-called 3’ IVT arrays) were probeset based: a certain fixed group of probes were part of a probeset which represented a certain gene or transcript (note however, that a gene can be represented by multiple probesets).

The more recent “Gene” and “Exon” Affymetrix arrays are exon based and hence there are two levels of summarization to get to the gene level. The “probeset” summarization leads to the exon level. The gene / transcript level is given by “transcript clusters”. Hence, the appropriate annotation package for our chip type is called hugene10sttranscriptcluster.db.

“Gene” arrays were created as affordable versions of the “Exon” arrays, by only taking the “good” probes from the Exon array. Initially on the Exon array, at least four probes were part of one “Exon”. With the thinned out “Gene” array, many probesets were made up of three or fewer probes. This is visualized in Figure 4: Single probesets are indicated by single colours; probesets representing one gene are indicated by a colour shade: e.g., all yellow probes belong to one Exon, and all yellow, orange and red probesets belong to one gene:

Figure 4: Visualization of the difference between “Exon” type array (left) and “Gene” type array (right)

On the left side, we see plenty of probes for each Exon / probeset (i.e. each colour): therefore, a summarization on the probeset / exon level makes sense. In the gene type array, however, only a small proportion of the original probes per probeset is included. Thus, a summarization on the probeset / exon level is not recommended for “Gene” arrays but nonetheless possible by using the hugene10stprobeset.db annotation package.

Note that furthermore, there are also no longer designated match/mismatch probes present on “Gene” and “Exon” type chips. The mismatch probe was initially intended as base-level for background correction, but hasn’t prevailed due to more elaborate background correction techniques that do not require a mismatch probe.

7.5 One-step preprocessing in oligo

The package oligo allows us to perform background correction, normalization and summarization in one single step using a deconvolution method for background correction, quantile normalization and the RMA (robust multichip average) algorithm for summarization.

This series of steps as a whole is commonly referred to as RMA algorithm, although strictly speaking RMA is merely a summarization method (6–8).

8.1 Computing the RLE

The RLE is performed by calculating the median log2 intensity of every transcript across all arrays.

We do this by calculating the row medians of exprs(palmieri_eset), as the transcripts are represented by the rows and the single microarrays by the columns.

Note that we do not have to apply the log2 manually, as the output data of the RMA function is in log2 scale by default.

We then substract this transcript median intensity from every transcript intensity via the sweep function.

8.2 Plotting the RLE

We finally reshape the data into a format in which we can use to create a boxplot for each array, as before:

row_medians_assayData <- 
  Biobase::rowMedians(as.matrix(Biobase::exprs(palmieri_eset)))

RLE_data <- sweep(Biobase::exprs(palmieri_eset), 1, row_medians_assayData)

RLE_data <- as.data.frame(RLE_data)
RLE_data_gathered <- 
  tidyr::gather(RLE_data, patient_array, log2_expression_deviation)

ggplot2::ggplot(RLE_data_gathered, aes(patient_array,
                                       log2_expression_deviation)) + 
  geom_boxplot(outlier.shape = NA) + 
  ylim(c(-2, 2)) + 
  theme(axis.text.x = element_text(colour = "aquamarine4", 
                                  angle = 60, size = 6.5, hjust = 1 ,
                                  face = "bold"))

Figure 5: Boxplot for the RLE values

Note that the y-axis now displays for each microarray the deviation of expression intensity from the median expression of the respective single transcripts across arrays.

Boxes with a larger extension therefore indicate an unusually high deviation from the median in a lot of transcripts, suggesting that these arrays are different from most of the others in some way.

Boxes that are shifted in y-direction indicate a systematically higher or lower expression of the majority of transcripts in comparison to most of the other arrays. This could be caused by quality issues or batch effects.

Therefore, if shape and median of a given box varies too much from the bulk, they should be inspected and potentially removed.

By inspecting the boxplot in Figure 5, five arrays could be considered as outliers: 2826_I, 2826_II, 3262_II, 3302_II and 3332_II are negatively y-shifted.

We will keep these samples in mind for heatmap cluster analysis later on in the workflow. Arrays that are confirmed to be outliers by heatmap analysis could be removed for subsequent analysis.

9.1 Some mathematical background on normalization (calibration) and background correction

A generic model for the value of the intensity \(Y\) of a single probe on a microarray is given by

\[ Y = B + \alpha \cdot S \]

where B is a random quantity due to background noise, usually composed of optical effects and non-specific binding, \(\alpha\) is a gain factor, and \(S\) is the amount of measured specific binding. The signal \(S\) is considered a random variable as well and accounts for measurement error and probe effects. The measurement error is typically assumed to be multiplicative so we can write:

\[ \log(S) = \theta + \varphi + \varepsilon \]

Here \(\theta\) represents the logarithm of the true abundance, \(\varphi\) is a probe-specific effect, and \(\varepsilon\) accounts for the nonspecific error. This is the additive-multiplicative-error model for microarray data used by RMA and also the vsn algorithm (10). The algorithms differ in the way that \(B\) is removed and an estimate of \(\theta\) is obtained.

9.2 Quality assessment of the calibrated data

We now produce a clustering heatmap and another PCA plot using the calibrated data.

PCA analysis

First, we perform a PCA analysis of the calibrated data analogously to the one with the raw data:

exp_palmieri <- Biobase::exprs(palmieri_eset_norm)
PCA <- prcomp(t(exp_palmieri), scale = FALSE)

percentVar <- round(100*PCA$sdev^2/sum(PCA$sdev^2),1)
sd_ratio <- sqrt(percentVar[2] / percentVar[1])

dataGG <- data.frame(PC1 = PCA$x[,1], PC2 = PCA$x[,2],
                    Disease = 
                     Biobase::pData(palmieri_eset_norm)$Factor.Value.disease.,
                    Phenotype = 
                     Biobase::pData(palmieri_eset_norm)$Factor.Value.phenotype.)


ggplot(dataGG, aes(PC1, PC2)) +
      geom_point(aes(shape = Disease, colour = Phenotype)) +
  ggtitle("PCA plot of the calibrated, summarized data") +
  xlab(paste0("PC1, VarExp: ", percentVar[1], "%")) +
  ylab(paste0("PC2, VarExp: ", percentVar[2], "%")) +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = sd_ratio) +
  scale_shape_manual(values = c(4,15)) + 
  scale_color_manual(values = c("darkorange2", "dodgerblue4"))

Figure 6: PCA plot of the calibrated, summarized data

In comparison to the first PCA analysis before RMA (Figure 2), we see that now the first principal component separates between the tissues types (Figure 6). This indicates that now differential expression between the tissue types is the dominant source of variation. Note that the second principal component separates the diseases.

Heatmap clustering analysis

We want to plot a heatmap with the sample-to-sample distances with the sample names as row-names. Also, we want to see how well the samples cluster for phenotype (inflamed and non-inflamed tissue) and disease (UC and CD), respectively. We use row annotation for that: it means that these features get a colour code and will be displayed to the left of each row.

phenotype_names <- ifelse(str_detect(pData
                                    (palmieri_eset_norm)$Factor.Value.phenotype.,
                             "non"), "non_infl.", "infl.")

disease_names <- ifelse(str_detect(pData
                                    (palmieri_eset_norm)$Factor.Value.disease.,
                             "Crohn"), "CD", "UC")

annotation_for_heatmap <- 
  data.frame(Phenotype = phenotype_names,  Disease = disease_names)

row.names(annotation_for_heatmap) <- row.names(pData(palmieri_eset_norm))

In order to map the sample-to-sample distances, we first compute the distances using the dist function. We need to transpose the expression values since the function computes the distances between the rows (i.e. genes in our case) by default. The default distance is the Euclidean one. However this can be changed and we choose the Manhattan distance here (it uses absolute distances along rectangular paths instead of squared distances of the direct path), as it is more robust. We set the diagonal of the distance matrix to NA in order to increase the contrast of the color coding. Those diagonal entries do not contain information since the distance of a sample to itself is always equal to zero.

dists <- as.matrix(dist(t(exp_palmieri), method = "manhattan"))

rownames(dists) <- row.names(pData(palmieri_eset_norm))
hmcol <- rev(colorRampPalette(RColorBrewer::brewer.pal(9, "YlOrRd"))(255))
colnames(dists) <- NULL
diag(dists) <- NA

ann_colors <- list(
  Phenotype = c(non_infl. = "chartreuse4", infl. = "burlywood3"),
  Disease = c(CD = "blue4", UC = "cadetblue2")
                   )
pheatmap(dists, col = (hmcol), 
         annotation_row = annotation_for_heatmap,
         annotation_colors = ann_colors,
         legend = TRUE, 
         treeheight_row = 0,
         legend_breaks = c(min(dists, na.rm = TRUE), 
                         max(dists, na.rm = TRUE)), 
         legend_labels = (c("small distance", "large distance")),
         main = "Clustering heatmap for the calibrated samples")

Figure 7: Heatmap of the sample-to-sample distances

On the heatmap plot (Figure 7) we also see that the samples do not cluster strongly by tissue, confirming the impression from the PCA plot (Figure 6) that the separation between the tissues is not perfect. The yellow stripes in the heatmap might correspond to outliers that could potentially be removed: the ones that could be flagged here are 2826_II, 3262_II, 3271_I, 2978_II and 3332_II. 2826_II, 3262_II and 3332_II were found to be outliers in both RLE and heatmap analysis and could therefore potentially be removed; however, in order to stay as close as possible to the original paper, we continue with the complete set of samples. Note again that more elaborate metrics to identify and remove outliers are provided by the arrayQualityMetrics package.

11.1 Removing multiple mappings

Many transcript-cluster identifiers will map to multiple gene symbols, i.e. they can’t be unambigously assigned.

We compute a summary table in the code below to see how many there are:

anno_grouped <- group_by(anno_palmieri, PROBEID)
anno_summarized <- 
  dplyr::summarize(anno_grouped, no_of_matches = n_distinct(SYMBOL))

head(anno_summarized)

   # A tibble: 6 x 2
     PROBEID no_of_matches
     <chr>           <int>
   1 7896742             3
   2 7896754             1
   3 7896759             1
   4 7896761             1
   5 7896779             1
   6 7896798             1

anno_filtered <- filter(anno_summarized, no_of_matches > 1)

head(anno_filtered)

   # A tibble: 6 x 2
     PROBEID no_of_matches
     <chr>           <int>
   1 7896742             3
   2 7896937             3
   3 7896961             3
   4 7897006             2
   5 7897632             3
   6 7897774             2

probe_stats <- anno_filtered 

nrow(probe_stats)

   [1] 1763

First, we grouped anno_palmieri by their PROBEID; that way, the subsequent operations are not carried through for each single row, but for each group, i.e. each PROBEID. We then summarized the groups and indicate the number of different genes assigned to a transcript cluster in the column no_of_matches. Finally, we filtered for PROBEIDs with multiple matches, i.e. no_of_matches > 1.

With dim(probe_stats), we could see how many probes have been mapped to multiple genes.

We have close to 2000 transcript clusters that map to multiple gene symbols. It is difficult to decide which mapping is “correct”. Therefore, we exclude these transcript clusters.

We want to remove those probe IDs that match the ones in probe_stats, as those are the probes with multiple mappings. We assign these IDs to the variable ids_to_exclude. Then, we generate palmieri_final, an expression set without the ids_to_exclude.

ids_to_exlude <- (featureNames(palmieri_manfiltered) %in% probe_stats$PROBEID)

table(ids_to_exlude)

   ids_to_exlude
   FALSE  TRUE 
   21041  1763

palmieri_final <- subset(palmieri_manfiltered, !ids_to_exlude)

validObject(palmieri_final)

   [1] TRUE

As we have just excluded probe IDs from the assay data, we now have to also exclude them from the feature data anno_palmieri:

head(anno_palmieri)

        PROBEID       SYMBOL
   2657 7896742  LINC00266-1
   2658 7896742       PCMTD2
   2659 7896742    LINC01881
   2664 7896754 LOC100287497
   2665 7896759    LINC01128
   2666 7896761       SAMD11
                                                                            GENENAME
   2657                                 long intergenic non-protein coding RNA 266-1
   2658 protein-L-isoaspartate (D-aspartate) O-methyltransferase domain containing 2
   2659                                  long intergenic non-protein coding RNA 1881
   2664                                                          septin 7 pseudogene
   2665                                  long intergenic non-protein coding RNA 1128
   2666                                     sterile alpha motif domain containing 11

Recall that fData enables us to access the feature data of an expression set. Until now, no feature data whatsoever is stored in the fData(palmieri_final). Only the row names are the row names of the assay data by default, which are the PROBEIDs of the transcripts.

Therefore, we generate a column PROBEID in fData(palmieri_final) and assign the row names of fData(palmieri_final) to it:

fData(palmieri_final)$PROBEID <- rownames(fData(palmieri_final))

Then, we left-join fData(palmieri_final)with anno_palmieri, which already contains the columns “SYMBOL” and “GENENAME”. A left-join keeps the rows and columns of the first argument and adds the corresponding column entries of the second argument:

fData(palmieri_final) <- left_join(fData(palmieri_final), anno_palmieri)

   Joining, by = "PROBEID"

# restore rownames after left_join
rownames(fData(palmieri_final)) <- fData(palmieri_final)$PROBEID 
    
validObject(palmieri_final)

   [1] TRUE

By left-joining with anno_palmieri, we thus add the “SYMBOL” and “GENENAME” columns from anno_palmieri for only the PROBEIDs that are in fData(palmieri_final) and thus get the feature Data for the filtered probes.

11.2 Building custom annotations

Alternatively, one can re-map the probes of the array to a current annotation. A workflow to do this for Illumina arrays is given in (13). Essentially, the individual probe sequences are re-aligned to an in-silico “exome” that consists of all annotated transcript exons.

In any case, the package pdInfoBuilder can be used to build customannotation packages for use with oligo. In order to do this, PGF / CLF files (called “Library files” on the Affymetrix website) as well as the probeset annotations are required. The probesets typically represent small stretches of the genome (such as a single exon) and multiple probesets are then used to form a transcript cluster.

The CLF file contains information about the location of individual probes on the array. The PGF file then contains the individual probe sequences and shows the probeset they belong to. Finally, the probeset annotation .csv then contains information about which probesets are used in which transcript cluster. Commonly, multiple probesets are used in one transcript cluster and some probesets are contained in multiple transcript clusters.

12.1 Linear models for microarrays

We now apply linear models to microarrays. Specifically, we discuss how to use the limma package for differential expression analysis. The package is designed to analyze complex experiments involving comparisons between many experimental groups simultaneously while remaining reasonably easy to use for simple experiments. The main idea is to fit a linear model to the expression data for each gene. Empirical Bayes and other methods are used to borrow information across genes for the residual variance estimation leading to “moderated” \(t\)-statistics, and stabilizing the analysis for experiments with just a small number of arrays (12). Conceptually, the final per gene variance is a mix of a prior variance and the per gene variance.

Typical experimental designs are disussed in chapter 9 of limma “User Guide”, which can be found on the Bioconductor landing page of limma.

In the following, we use appropriate design and contrast matrices for our linear models and fit a linear model to each gene separately.

12.2 A linear model for the data

For the subsequent linear modelling of the data, we introduce the abbreviations “UC” and “CD” for the disease types, and “non_infl.” and “infl.” for the phenotypes, respectively:

individual <- 
  as.character(Biobase::pData(palmieri_final)$Characteristics.individual.)

tissue <- str_replace_all(Biobase::pData(palmieri_final)$Factor.Value.phenotype.,
                  " ", "_")

tissue <- ifelse(tissue == "non-inflamed_colonic_mucosa",
                 "nI", "I")

disease <- 
  str_replace_all(Biobase::pData(palmieri_final)$Factor.Value.disease.,
                  " ", "_")
disease <- 
  ifelse(str_detect(Biobase::pData(palmieri_final)$Factor.Value.disease., 
                    "Crohn"), "CD", "UC")

The original paper is interested in changes in transcription that occur between inflamed and adjacent non-inflamed mucosal areas of the colon. This is studied in both inflammatory bowel disease types.

For building our linear model, we have to think about which experimental variables we want to consider. As we want to find differential expression between non-inflamed and inflamed tissue, in principle, those are the only two variables we would have to consider.

However, since we have two arrays per individual patient, we have a “Paired Samples” design (see section 9.4 of the limma user guide). This means that the samples might be biased by the person they come from. Whenever a feature in an experimental setup is expected to have a systematic influence on the result, blocking factors on these features should be introduced.

Thus, the first factor we need is a blocking factor for the individuals that will absorb differences in expression between them. Therefore, we block on patients, which means that the patient IDs become variables of the linear model.

Then we create factors that give us the grouping for the tissue types (non-inflamed and inflamed).

Finally, we create two design matrices, one for each of the two diseases as we will analyze them separately in order to follow the analysis strategy of the original paper closely (one could also fit a joint model to the complete data set; however, the two diseases might show very different characteristics so that a joint fit might not be appropriate).

i_CD <- individual[disease == "CD"]
design_palmieri_CD <- model.matrix(~ 0 + tissue[disease == "CD"] + i_CD)
colnames(design_palmieri_CD)[1:2] <- c("I", "nI")
rownames(design_palmieri_CD) <- i_CD 

i_UC <- individual[disease == "UC"]
design_palmieri_UC <- model.matrix(~ 0 + tissue[disease == "UC"] + i_UC )
colnames(design_palmieri_UC)[1:2] <- c("I", "nI")
rownames(design_palmieri_UC) <- i_UC 

We can inspect the design matrices:

head(design_palmieri_CD[, 1:6])

        I nI i_CD183 i_CD2114 i_CD2209 i_CD2255
   164  0  1       0        0        0        0
   164  1  0       0        0        0        0
   183  0  1       1        0        0        0
   183  1  0       1        0        0        0
   2114 0  1       0        1        0        0
   2114 1  0       0        1        0        0

head(design_palmieri_UC[, 1:6])

        I nI i_UC2424 i_UC2987 i_UC2992 i_UC2995
   2400 0  1        0        0        0        0
   2400 1  0        0        0        0        0
   2424 0  1        1        0        0        0
   2424 1  0        1        0        0        0
   2987 0  1        0        1        0        0
   2987 1  0        0        1        0        0

In the design matrix, the rows represent the patient array, and the columns are the variables we include in our linear model. The variables correspond to our blocking factors: there are two for non-inflamed and inflamed tissue, respectively, and one for each patient. “i_UC2424” for example is the blocking variable of patient 2424; UC stands for the disease the patient is suffering from. For example, the first two rows of the design matrix design_palmieri_CD correspond to the two arrays for individual “164”.

The design matrix entries are 0 or 1 and thus tell us which variables are “active” for which sample:

It can be seen as a switch that turns the single variables “on” (with a 1 at the corresponding position) or “off” (with a 0 at the corresponding position) for each row, i. e. each patient array. If we take a closer look at the single rows, we see that for each sample, there are two “ones” assigned: one to one of the variables “nI” or “I” corresponding to the tissue the array came from, and one to the corresponding patient-specific blocking variable.

Note that in the linear model, individual 164 serves as the baseline for all other individuals and thus isn’t included in the sample variables.

Note that instead of blocking on individuals, it would also be possible to use a “mixed model” approach with the duplicateCorrelation()function from the limma package. It has advantages over the “fixed patient effect model” presented here in terms of applicability to more complicated experimental designs, where we want to perform comparisons both within and between the patients (e.g. comparing between the two diseases; “split-plot-designs”).

More information on it can be found in the limma User’s Guide (section 17.3.6). However, the above explained is more intuitive and is therefore used here.

Before heading on to find all differentially expressed genes for both diseases, we will first have a look at how this works in principle for one gene. We will fit the linear model for one gene and run a t-test in order to see whether the gene is differentially expressed or not.

12.3 Analysis of differential expression based on a single gene

For linear model fitting and subsequent testing for differential expression by t-test, we will pick the gene with the PROBEID 8164535. It has the gene symbol CRAT and will be named as such in the following code.

Illustration of the fitted linear model on the CRAT gene

Before fitting the linear model, we have a look at the expression intensities of this gene for each patient in non-inflamed and inflamed tissue, respectively:

tissue_CD <- tissue[disease == "CD"]
crat_expr <- Biobase::exprs(palmieri_final)["8164535", disease == "CD"]
crat_data <- as.data.frame(crat_expr)
colnames(crat_data)[1] <- "org_value"
crat_data <- mutate(crat_data, individual = i_CD, tissue_CD)

crat_data$tissue_CD <- factor(crat_data$tissue_CD, levels = c("nI", "I"))

ggplot(data = crat_data, aes(x = tissue_CD, y = org_value, 
                             group = individual, color = individual)) +
      geom_line() +
      ggtitle("Expression changes for the CRAT gene")

Figure 10: Visualization of expression changes

We see that overall, this gene is expressed less in inflamed tissue (Figure 10). We also see that the absolute expression intensities vary greatly between patients. However, we have already taken care of this problem by introducing blocking factors based on the individuals, which allows us to compare the tissues within each individual as represented by the single lines.

If we had not blocked for individuals, the linear model would treat them interchangably and a graphical depiction would only include a single line. Since the individuals haver very different baseline expression levels, this would lead to a very high variance of the estimated fold changes.

We now compute the variable coefficients by fitting a linear model. We get a vector crat_coef with one entry for each variable.

crat_coef <- lmFit(palmieri_final[,disease == "CD"],
                design = design_palmieri_CD)$coefficients["8164535",]

crat_coef

          I       nI  i_CD183 i_CD2114 i_CD2209 i_CD2255  i_CD255 i_CD2826 
    6.76669  7.19173  0.12382 -0.22145  0.55759 -0.39905  0.29204 -1.07285 
   i_CD2853 i_CD2978  i_CD321 i_CD3262 i_CD3266 i_CD3271 i_CD3302 i_CD3332 
   -0.78285 -0.11633  0.01692 -0.62480 -0.46209 -0.61732 -0.30257 -0.09709

In order to now obtain the expression values fitted by the model, we multiply the design matrix with this vector of coefficients crat_coef:

crat_fitted <- design_palmieri_CD %*% crat_coef
rownames(crat_fitted) <- names(crat_expr)
colnames(crat_fitted) <- "fitted_value"

crat_fitted

               fitted_value
   164_I_.CEL         7.192
   164_II.CEL         6.767
   183_I.CEL          7.316
   183_II.CEL         6.891
   2114_I.CEL         6.970
   2114_II.CEL        6.545
   2209_A.CEL         7.749
   2209_B.CEL         7.324
   2255_I.CEL         6.793
   2255_II.CEL        6.368
   255_I.CEL          7.484
   255_II.CEL         7.059
   2826_I.CEL         6.119
   2826_II.CEL        5.694
   2853_I.CEL         6.409
   2853_II.CEL        5.984
   2978_I.CEL         7.075
   2978_II.CEL        6.650
   321_I.CEL          7.209
   321_II.CEL         6.784
   3262_I.CEL         6.567
   3262_II.CEL        6.142
   3266_I.CEL         6.730
   3266_II.CEL        6.305
   3271_I.CEL         6.574
   3271_II.CEL        6.149
   3302_I.CEL         6.889
   3302_II.CEL        6.464
   3332_I.CEL         7.095
   3332_II.CEL        6.670

Recall that for every row in the design matrix (i.e. every patient sample) only the variables with a 1 in the design matrix are taken into account for calculating the fitted expression value.

This means that as output of the multiplication, we get a vector crat_fitted whose entries are the sum of relevant variable coefficients for each sample, respectively.

For example, the fitted value for patient sample 2114_I.CEL is 6.9703: it is the sum of the respective activated variable coefficients “nI” (7.1917) and “i_CD2114” (-0.2215).

Let’s visualize the difference between non-inflamed and inflamed tissue again after fitting:

crat_data$fitted_value <- crat_fitted

ggplot(data = crat_data, aes(x = tissue_CD, y = fitted_value, 
                             group = individual, color = individual)) +
      geom_line() +
      ggtitle("Fitted expression changes for the CRAT gene")

Figure 11: Expression changes for the CRAT gene

Note that the difference of the fitted expression values between inflamed and non-inflamed samples of one patient is the same for all patients and is determined by the difference between the variable coefficients of I (6.7667) and nI (7.1917), which is -0.425 (Figure 11).

This is the case because the same blocking variable is activated by the design matrix for both samples from a single patient, leading to a comparison within patients only. These blocking variables correct the fitted tissue specific expression values towards the expression levels of the individual patients. Therefore the final estimate is like an average of all the within-individual differences.

The “difference” between non-inflamed and inflamed tissue of -0.425 is actually a log2 fold change, as our expression data is on the log2 scale. -0.425 therefore is our log2 fold change for the CRAT gene.

Differential expression analysis of the CRAT gene

In order to test whether the gene is differentially expressed or not, a \(t\)-test with the null hypothesis that there is no difference in the expression between non-inflamed and inflamed tissue is carried out. Our blocking design is conceptually similar to a paired t-test for which the statistic is given by:

\[ t = \frac{\bar{d} }{s/\sqrt{n}} \]

Where, \(\bar{d}\) is the mean difference in expression values between the individuals. The paired t-test computes the variance \(s^2\) from the paired differences. This is lower than the variance of a standard t-test and thus the paired t-test has higher power as long as the expression values for the same individual are correlated (see e.g. the article on Wikipedia).

We thus have improved the power of the ordinary \(t\)-test by reducing the variance via blocking on individuals.

We now conduct the \(t\)-test on the linear model in order to find out whether the difference between non-inflamed and inflamed tissue differs significantly from 0:

crat_noninflamed <- na.exclude(crat_data$org_value[tissue == "nI"])
crat_inflamed <- na.exclude(crat_data$org_value[tissue == "I"])
res_t <- t.test(crat_noninflamed ,crat_inflamed , paired = TRUE)
res_t

   
    Paired t-test
   
   data:  crat_noninflamed and crat_inflamed
   t = 6.8, df = 14, p-value = 8e-06
   alternative hypothesis: true difference in means is not equal to 0
   95 percent confidence interval:
    0.2919 0.5581
   sample estimates:
   mean of the differences 
                     0.425

We get a low p-value close to 0 and thus can conclude that the CRAT gene is differentially expressed between non-inflamed and inflamed tissue.

Note that the p-value isn’t exactly the same one as below when analyzing the differential expression of all genes. This is due to the variance moderation performed by limma .

12.4 Contrasts and hypotheses tests

We now fit the linear model for all genes and define appropriate contrasts to test hypotheses of interest.

We want to compare the inflamed to the non-inflamed tissue. Thus, we create a contrast matrix consisting of only one contrast “I-nI”: limma’s function makeContrasts creates this matrix from a symbolic description of the contrast of interest.

We now fit a linear model to our data and apply the contrasts.fit() function to it in order to find genes with significant differential expression between non-inflamed and inflamed tissue:

contrast_matrix_CD <- makeContrasts(I-nI, levels = design_palmieri_CD)

palmieri_fit_CD <- eBayes(contrasts.fit(lmFit(palmieri_final[,disease == "CD"],
                                design = design_palmieri_CD),
                                contrast_matrix_CD))

contrast_matrix_UC <- makeContrasts(I-nI, levels = design_palmieri_UC)

palmieri_fit_UC <- eBayes(contrasts.fit(lmFit(palmieri_final[,disease == "UC"],
                                design = design_palmieri_UC),
                                contrast_matrix_UC))

We applied the empirical Bayes variance moderation method to the model via the eBayes() function, which computes moderated \(t\)-statistics. In microarray analysis, the number of arrays often is quite small, and thus variance estimation is difficult. Using a combination of the per-gene-variance and a prior variance we can improve the variance estimate, hence the term “moderation”. “Empirical Bayes” means that the prior is estimated from the data.

The result of the eBayes() step is that the individual variances are shrunken towards the prior value.

12.5 Extracting results

Finally, we extract the number of differentially expressed genes. Results can be extracted by use of the topTable function. We extract the results for both Crohn’s disease and ulcerative colitis, and the results are sorted by their absolute \(t\)-statistics. As a diagnostic check, we also plot the p-value histogram: We expect a uniform distribution for the p-values that correspond to true null hypotheses, while a peak near zero shows an enrichment for low p-values corresponding to differentially expressed (DE) genes.

Note that if the p-value distribution for a dataset is very different from the ones in the histograms below, this might lead to quality loss in the subsequent analysis. Reasons for a divergent p-value-distribution might be batch effects or a lack of consideration of other blocking factors in the design matrix. Thus, if the p-value is not as expected, try to include possible blocking factors and batches and rerun the analysis. If this does not help, empirical Bayes / null estimation methods for multiple testing are useful.

A good starting point to learn about these methods is the article on false discovery rate estimation by Korbininan Strimmer (15) and chapter 1–6 of Efron’s book on Large-Scale Inference (16), as well as the blog-post on “How to interpret a p-value histogram” by David Robinson (17).

table_CD <- topTable(palmieri_fit_CD, number = Inf)
head(table_CD)

           PROBEID  SYMBOL                            GENENAME   logFC AveExpr
   7928695 7928695  PRXL2A               peroxiredoxin like 2A -0.5990   7.739
   8123695 8123695    ECI2         enoyl-CoA delta isomerase 2 -0.4855   6.876
   8164535 8164535    CRAT       carnitine O-acetyltransferase -0.4250   6.732
   8009746 8009746 SLC16A5   solute carrier family 16 member 5 -0.5182   5.562
   7952249 7952249    USP2      ubiquitin specific peptidase 2 -0.8484   5.606
   8105348 8105348    GPX8 glutathione peroxidase 8 (putative)  0.8312   5.301
                t   P.Value adj.P.Val     B
   7928695 -7.059 1.383e-06   0.02092 5.305
   8123695 -6.317 5.907e-06   0.02092 4.027
   8164535 -6.230 7.037e-06   0.02092 3.872
   8009746 -6.206 7.384e-06   0.02092 3.829
   7952249 -6.203 7.426e-06   0.02092 3.824
   8105348  6.074 9.653e-06   0.02092 3.590

hist(table_CD$P.Value, col = brewer.pal(3, name = "Set2")[1],
     main = "inflamed vs non-inflamed - Crohn's disease", xlab = "p-values")

Figure 12: Histogram of the p–values for Crohn’s disease

table_UC <- topTable(palmieri_fit_UC, number = Inf)
head(table_UC)

           PROBEID  SYMBOL                            GENENAME   logFC AveExpr
   8003875 8003875   SPNS2          sphingolipid transporter 2  0.7412   6.478
   8082012 8082012 SLC15A2   solute carrier family 15 member 2  0.8061   5.319
   7952290 7952290  TRIM29      tripartite motif containing 29  1.0140   5.855
   7961693 7961693    LDHB             lactate dehydrogenase B  0.3968   9.534
   8072015 8072015    GRK3 G protein-coupled receptor kinase 3  0.4713   5.584
   8096070 8096070    BMP3        bone morphogenetic protein 3 -1.6961   6.420
                t   P.Value adj.P.Val     B
   8003875  7.801 4.552e-07  0.003984 6.441
   8082012  7.744 5.032e-07  0.003984 6.352
   7952290  7.482 8.008e-07  0.003984 5.937
   7961693  7.401 9.263e-07  0.003984 5.806
   8072015  7.308 1.097e-06  0.003984 5.654
   8096070 -7.289 1.136e-06  0.003984 5.623

hist(table_UC$P.Value, col = brewer.pal(3, name = "Set2")[2],
     main = "inflamed vs non-inflamed - Ulcerative colitis", xlab = "p-values")

Figure 13: Histogram of the p–values for ulcerative colitis

12.6 Multiple testing FDR, and comparison with results from the original paper

In the original paper, a p-value of 0.001 was used as a significance cutoff. Using this we get 947 genes identified as differentially expressed for UC:

nrow(subset(table_UC, P.Value < 0.001))

   [1] 947

However, it is impossible to determine a precise bound on the number of false positive genes in this list. All that we can say using p-values is that we have at most 21041 (total number of tests) * 0.001 = 21.041 false positive genes in our list. Therefore, by choosing a p-value cutoff of 0.001, as much as 2.22% of our genes identified as differentially expressed might be false positives.

Thus, we can see that the “raw” p-values are very “liberal” when looking at many tests simultaneously. We therefore need error rates adapted to the multiple testing situation. By far the most popular one in molecular biology is the false discovery rate or FDR for short. It is the percentage of false positives among all positives. As we have seen, the FDR of our genes list using a simple p-value cutoff might be quite high.

On the other hand, we can see a clear peak in the p-value histogram (Figure 13), caused by the differentially expressed genes. There we expect the actual FDR of our list to be lower.

The FDR at a given cutoff is given by the “adjusted” p-value in the results table.

tail(subset(table_UC, P.Value < 0.001))

           PROBEID   SYMBOL
   7915640 7915640   EIF2B3
   7894577 7894577     <NA>
   7897877 7897877 TNFRSF1B
   8142671 8142671     WASL
   7941946 7941946   NDUFV1
   8140371 8140371 TMEM120A
                                                            GENENAME   logFC
   7915640 eukaryotic translation initiation factor 2B subunit gamma  0.2795
   7894577                                                      <NA>  0.4048
   7897877                        TNF receptor superfamily member 1B  0.3564
   8142671                             Wiskott-Aldrich syndrome like -0.4281
   7941946            NADH:ubiquinone oxidoreductase core subunit V1 -0.3832
   8140371                                transmembrane protein 120A -0.3657
           AveExpr      t   P.Value adj.P.Val       B
   7915640   6.271  3.957 0.0009867   0.02204 -0.6476
   7894577   5.598  3.956 0.0009891   0.02205 -0.6500
   7897877   6.656  3.956 0.0009898   0.02205 -0.6506
   8142671   7.907 -3.955 0.0009905   0.02205 -0.6513
   7941946   8.867 -3.952 0.0009975   0.02219 -0.6579
   8140371   6.519 -3.951 0.0009989   0.02220 -0.6592

The adjusted p-value for a raw p-value of 0.001 in the table is 0.0222, which is an order of magnitude lower than the FDR we can infer from p-values alone.

So although this is not recommended in general, we also use a p-value cutoff at 0.001 in the following in order to be able to compare our workflow results to the paper results.

The paper results can be downloaded as excel files from http://links.lww.com/IBD/A795 and should be saved as an .xlsx file named palmieri_DE_res.xlsx in your working directory.

Note that while the paper uses p-value cutoffs, it also reports the corresponding FDRs (just as we did for the UC data here).

For a p-value cutoff of 0.001, the corresponding FDRs are 0.05 in Crohn’s disease and 0.02 in ulcerative colitis. There are four tables in total, giving the list of up and downregulated genes in CD and UC, respectively. We calculate the overlap between our results and the ones from the paper as the ratio of the genes that were found in both analyses and the genes that were only found in the paper.

We also calculate the total number of diffentially expressed genes that we find in our workflow analysis.

fpath <- system.file("extdata", "palmieri_DE_res.xlsx", package = "maEndToEnd")
palmieri_DE_res <- sapply(1:4, function(i) read.xlsx(cols = 1, fpath, 
                                                     sheet = i, startRow = 4))

names(palmieri_DE_res) <- c("CD_UP", "CD_DOWN", "UC_UP", "UC_DOWN")
palmieri_DE_res <- lapply(palmieri_DE_res, as.character)
paper_DE_genes_CD <- Reduce("c", palmieri_DE_res[1:2])
paper_DE_genes_UC <- Reduce("c", palmieri_DE_res[3:4])

overlap_CD <- length(intersect(subset(table_CD, P.Value < 0.001)$SYMBOL,  
                               paper_DE_genes_CD)) / length(paper_DE_genes_CD)


overlap_UC <- length(intersect(subset(table_UC, P.Value < 0.001)$SYMBOL,
                               paper_DE_genes_UC)) / length(paper_DE_genes_UC)
overlap_CD

   [1] 0.6443

overlap_UC 

   [1] 0.6731

total_genenumber_CD <- length(subset(table_CD, P.Value < 0.001)$SYMBOL)
total_genenumber_UC <- length(subset(table_UC, P.Value < 0.001)$SYMBOL)

total_genenumber_CD

   [1] 576

total_genenumber_UC

   [1] 947

We find 576 (CD) and 947 (UC) differentially expressed genes (“DE-genes”).

In the paper, 298 (CD) and 520 (UC) DE-genes were found for the two diseases at the same cutoff. This higher number of DE-genes identified is probably due to the increased power of the blocking according to the individuals and the moderated variance estimation that limma performs.

We see that we get a moderate overlap of 0.6443 for CD and 0.6731 for UC, showing that both analyses lead to somewhat comparable results.

12.7 Visualization of DE analysis results - volcano plot

For a visualization of the differentially expressed genes, we create a volcano plot, which is commonly used to summarize the results of a differential expression analysis in a single figure.

For a better overview, we only show gene symbols of genes with a fold change greater than 1, which we define in the volcano_names object. The highlight option in the volcanoplot function is set to 100 and thus only labels the 100 genes with the lowest p-values.

volcano_names <- ifelse(abs(palmieri_fit_CD$coefficients)>=1, 
                        palmieri_fit_CD$genes$SYMBOL, NA)
             
             
volcanoplot(palmieri_fit_CD, coef = 1L, style = "p-value", highlight = 100, 
            names = volcano_names,
            xlab = "Log2 Fold Change", ylab = NULL, pch=16, cex=0.35)

Figure 14: Volcano plot of the DE-genes

We can now do a little research on the biological function of genes that show a high foldchange, for example the gene with the symbol S100A8 on the right side of the plot (Figure 14). If we search for this gene symbol on genecards.org, we find that it encodes for a protein that builds a pro-inflammatory complex in association with another protein.

13.1 Matching the background set of genes

The function genefinder from the genefilter package (21) will be used to find a background set of genes that are similar in expression to the differentially expressed genes. We then check whether the background has roughly the same distribution of average expression strength as the foreground.

We do this in order not to select a biased background since the gene set testing is performed by a simple Fisher test on a 2x2 table. Note that this approach is very similar to commonly used web tools like GOrilla (22).

For every differentially expressed gene, we try to find genes with similar expression with genefinder. The genefinder function returns a list with two elements for each gene: one with the indices of the background genes found and one with the distances to the DE-genes:

back_genes_idx <- genefilter::genefinder(palmieri_final, 
                                        as.character(DE_genes_CD), 
                                        method = "manhattan", scale = "none")

We have to extract the PROBEIDs, which correspond to the indices. We do that by using the sapply function, which gives us a single matrix with the DE-genes as column names and the PROBEIDs of the corresponding background genes in the cells below:

back_genes_idx <- sapply(back_genes_idx, function(x)x$indices)

We then create a vector back_genes containing all background gene PROBEIDs:

In order to eliminate foreground genes, i.e. DE-genes, from the back_genes set, we use the setdiff function. It returns all elements from the first argument (back_genes) that are not part of the second argument (DE_genes_CD). With the intersect function, we verify that we were successful: it should return 0, as there shouldn’t be any intersect anymore between back_genesand DE_genes_CD:

back_genes <- featureNames(palmieri_final)[back_genes_idx]
back_genes <- setdiff(back_genes, DE_genes_CD)

    
intersect(back_genes, DE_genes_CD)

   character(0)

length(back_genes)

   [1] 9756

We create a multidensity plot with mean expression on the x-axis and curves for all genes, foreground genes and background genes, respectively (Figure 15). We want to see whether the background genes show a plot similar to the foreground genes so that the background is not biased for the gene enrichment analysis:

multidensity(list(
        all = table_CD[,"AveExpr"] ,
        fore = table_CD[DE_genes_CD , "AveExpr"],
        back = table_CD[rownames(table_CD) %in% back_genes, "AveExpr"]),
        col = c("#e46981", "#ae7ee2", "#a7ad4a"),
     xlab = "mean expression",
   main = "DE genes for CD-background-matching")

Figure 15: Selecting a background set of genes for the gene ontology analysis

When comparing the “background gene” curve to the “foreground gene” curve, we see a similar curve shape, indicating a sensible background matching (Figure 15). Note that the right-shift of the “foreground-gene” curve in comparison to the “background-gene” curve indicates that DE-genes are generally very highly expressed, so that it wasn’t possible to find background-genes with exactly equal overall expression distribution.

The “all gene” curve has the leftmost curve maximum; this can be explained by a high number of lowly expressed genes in all samples and shows that a background matching is sensible in order to avoid biases.

For the actual testing of which GO gene sets are enriched in inflamed tissue, we use the topGO package which implements a nice interface to Fisher testing and also has additional algorithms taking the GO structure into account, by e.g. only reporting the most specific gene set in the hierarchy (23).

The GO has three top ontologies: Cellular component (CC), biological processes (BP), and molecular function (MF). For illustrative purposes we limit ourselves to the BP category here.

13.2 Running topGO

topGO requires a topGOdata object containing the necessary information for the analysis. We follow the steps described in the topGO vignettes: First, we will create a named vector all_genes with all genes to be analyzed, i.e. DE-genes and background genes:

gene_IDs <- rownames(table_CD)
in_universe <- gene_IDs %in% c(DE_genes_CD, back_genes)
in_selection <- gene_IDs %in% DE_genes_CD 

all_genes <- in_selection[in_universe]
all_genes <- factor(as.integer(in_selection[in_universe]))
names(all_genes) <- gene_IDs[in_universe] 

The following steps were carried through:

1. we created an in_universe vector by using the %in% matching function. We want to know which elements from gene_IDs are also contained in DE_genes_CD and back_genes, as the latter two are our gene universe we use for enrichment analysis. We got a vector in_universe with the length of gene_IDs that has the entry TRUE when the corresponding gene in gene_IDs could be also found in DE_genes_CD or back_genes, and FALSE otherwise.

2. We did the same for our DE-genes and call this vector in_selection.

3. We created the all_genes vector: a) First, we selected all the elements from in_selection that are TRUE in in_universe by applying all_genes <- in_selection[in_universe]. b) Then, we converted the elements in all_genes from TRUE and FALSE to 0 and 1 by converting the vector to an integer vector. This way, each element in the vector is a 0 if the corresponding gene is a background gene and a 1 if the corresponding gene is a DE-gene. Also, we converted the vector to a factor. c) We named the vector elements with the corresponding gene_IDs.

We now initialize the topGO data set, using the GO annotations contained in the annotation data base for the chip we are using. The nodeSize parameter specifies a minimum size of a GO category we want to use: i.e. here, categories with less than 10 genes are not included in the testing.

top_GO_data <- new("topGOdata", ontology = "BP", allGenes = all_genes,
 nodeSize = 10, annot = annFUN.db, affyLib = "hugene10sttranscriptcluster.db")

Now the tests can be run. topGO offers a wide range of options, for details see the paper (23) or the package vignette.

We run two common tests: an ordinary Fisher test for every GO category, and the “elim” algorithm, which tries to incorporate the hierarchical structure of the GO and tries to “decorrelate” it in order to report the most specific significant term in the hierarchy.

The algorithm starts processing the nodes / GO categories from the highest (bottommost) level and then iteratively moves to nodes from a lower level. If a node is scored as significant, all of its genes are marked as removed in all ancestor nodes. This way, the “elim” algorithm aims at finding the most specific node for every gene.

The test uses a 0.01 p-value cutoff by default.

result_top_GO_elim <- 
  runTest(top_GO_data, algorithm = "elim", statistic = "Fisher")
result_top_GO_classic <- 
  runTest(top_GO_data, algorithm = "classic", statistic = "Fisher")

We can now inspect the results. We look at the top 100 GO categories according to the “Fisher elim” algorithm. The function GenTable produces a table of significant GO categories, the function printGenes gives genes annotated to them; the significant ones are denoted with a “2” in the “raw p-value” column, the non-significant ones with a “1”. We therefore select raw p-value == 2.

Note that we do not get the actual p-values here because our all_genes vector doesn’t contain this information; it only tells us whether a gene is differentially expressed or not.

res_top_GO <- GenTable(top_GO_data, Fisher.elim = result_top_GO_elim,
        Fisher.classic = result_top_GO_classic,
        orderBy = "Fisher.elim" , topNodes = 100)

genes_top_GO <- printGenes(top_GO_data, whichTerms = res_top_GO$GO.ID,
    chip = "hugene10sttranscriptcluster.db", geneCutOff = 1000)

res_top_GO$sig_genes <- sapply(genes_top_GO, function(x){
                str_c(paste0(x[x$'raw p-value' == 2, "Symbol.id"],";"), 
                      collapse = "")
    })

head(res_top_GO[,1:8], 20)

           GO.ID                                        Term Annotated Significant
   1  GO:0030593                       neutrophil chemotaxis        61          39
   2  GO:0051897 positive regulation of protein kinase B ...       107          53
   3  GO:0006954                       inflammatory response       456         200
   4  GO:0032496              response to lipopolysaccharide       216          94
   5  GO:0006911                    phagocytosis, engulfment        48          27
   6  GO:0030335       positive regulation of cell migration       341         133
   7  GO:0030198           extracellular matrix organization       232          97
   8  GO:0050900                         leukocyte migration       291         138
   9  GO:0017144                      drug metabolic process       480         182
   10 GO:0070374 positive regulation of ERK1 and ERK2 cas...       149          58
   11 GO:0048661 positive regulation of smooth muscle cel...        66          35
   12 GO:0070098        chemokine-mediated signaling pathway        51          26
   13 GO:0016525         negative regulation of angiogenesis        72          33
   14 GO:0002675 positive regulation of acute inflammator...        20          14
   15 GO:0007186 G protein-coupled receptor signaling pat...       441         147
   16 GO:0045429 positive regulation of nitric oxide bios...        33          19
   17 GO:0030449         regulation of complement activation        44          23
   18 GO:0015718               monocarboxylic acid transport        79          39
   19 GO:0002548                         monocyte chemotaxis        27          19
   20 GO:0001774                  microglial cell activation        28          19
      Expected Rank in Fisher.classic Fisher.elim Fisher.classic
   1     13.40                     68     1.9e-10        2.0e-12
   2     23.51                    104     2.6e-10        2.6e-10
   3    100.18                      1     4.8e-08        8.1e-27
   4     47.45                     61     7.5e-08        7.1e-13
   5     10.55                    220     2.4e-07        2.4e-07
   6     74.91                     55     2.4e-07        3.3e-13
   7     50.97                     76     2.6e-07        5.5e-12
   8     63.93                      6     6.2e-07        1.5e-22
   9    105.45                     31     6.6e-07        2.9e-16
   10    32.73                    277     1.8e-06        1.8e-06
   11    14.50                    164     2.2e-06        3.1e-08
   12    11.20                    308     4.9e-06        4.9e-06
   13    15.82                    312     5.5e-06        5.5e-06
   14     4.39                    315     5.8e-06        5.8e-06
   15    96.88                    151     6.9e-06        1.2e-08
   16     7.25                    339     9.4e-06        9.4e-06
   17     9.67                    344     1.0e-05        1.0e-05
   18    17.36                    186     1.0e-05        6.8e-08
   19     5.93                    196     1.0e-05        1.0e-07
   20     6.15                    221     1.0e-05        2.5e-07

13.3 Visualization of the GO-analysis results

A graph of the results can also be produced. Here we visualize the three most significant nodes according to the Fisher elim algorithm in the context of the GO hierarchy.

showSigOfNodes(top_GO_data, score(result_top_GO_elim), firstSigNodes = 3,
               useInfo = 'def')

Figure 16: Significantly enriched GO nodes in the GO hierarchy

We can see that indeed GO categories related to inflammation, signalling and immune response come up as significant (Figure 16) Gene set enrichment analysis has been a field of very extensive research in bioinformatics. For additional approaches see the topGO vignette and the references therein and also in the GeneSetEnrichment view.

14.1 Visualizing the reactome based analysis results

The top pathways can be displayed as a bar chart that displays all categories with a p-value below the specified cutoff (Figure 17).

barplot(reactome_enrich)

Figure 17: Enriched Reactome pathways and their p–values as a bar chart

The “enrichment map” from the package enrichplot displays the results of the enrichment analysis as a graph, where the color represents the p-value of the pathway and the edge-thickness (that is the line connecting two pathways) is proportional to the number of overlapping genes between two pathways.

emapplot(reactome_enrich, showCategory = 10)

Figure 18: Enriched Reactome pathways enrichment results as a graph

Again, the graph in Figure 18 shows pathways related to signalling and immune response.

The package clusterProfiler (27) can also perform these analyses using downloaded KEGG data. Furthermore, the package EnrichmentBrowser (28) additionally offers network-based enrichment analysis of individual pathways. This allows the mapping of the expression data at hand to known regulatory interactions.