hypergeoDistMatrix <- matrix(c(10,490,90,9410),2,2)
hypergeoDistMatrix [,1] [,2]
[1,] 10 90
[2,] 490 9410Functional analysis
Functional analysis is the process of taking a list of genes which we cannot meaningfully interact with and distilling information that we as researchers can then interpret. Generally, this involves grouping genes based on functions or processes, and determining which processes or functions appear to be different between your experimental groups.
There are different approaches to functional analysis e.g.,:
Over-representation analysis: are certain processes, pathways or biological functions over-represented in a list of differentially expressed genes relative to the background genome (often called Gene Ontology analysis).
Functional class scoring (Gene Set Enrichment Analysis (GSEA)): Rather than taking a list of differentially expressed genes (which are themselves based on thresholds and cutoffs defined by a user), GSEA takes the ranked position of all genes to see how pathways and processes are perturbed.
Pathway topology: Pathway topology analysis takes a higher level approach, compiling information on up- and down-regulation at the level of the pathway rather than the level of the gene.
Today we will cover over-representation analysis and functional class scoring. Links to pathway topology workflows are provided.
Annotation
The basis of this stage of the analysis is to take individual genes, associate them with a function, and then to identify whether there are common themes of function in genes which have undergone changes in gene expression.
This first requires that our genes are annotated - we must be able to associate some type of function or process with a given gene. There are now many databases that maintain information and annotation on genes, including proposed and confirmed functions, pathways, processes etc.,. Annotation and databases is an entire area in itself, and not something we can cover today. For now, we will retrieve annotation information from different sources and attach the annotation to genes.
The annotations that we are going to be using today are called Gene Ontologies or GO terms. There are three major classes of GO – Molecular Function (MF), Cellular Component (CC), and Biological Process (BP). GO terms are hierarchical, and each child term is more specialised than the parent term. Check out the documentation here on gene ontologies.
Over-representation Analysis
The hypergeometric distribution
Let’s say you have a list of 100 differentially expressed genes and, after annotation, you notice that 10 of these genes are involved in apoptosis. Is this a significant finding? Initially, you might think that yes it is significant, since 10% of your genes of interest share a role. However, to understand whether this is significant we need to consider the number of apoptosis-related genes in the genome. If, for example, we find that 12% of all genes in the genome are annotated as apoptosis-related, then seeing 10% of your differentially expressed genes with this annotation isn’t significant - it’s approximately what you would have expected. However, if you’d found 25% (or 0%!) THEN you might have something significant.
Over-representation analysis is the use of a hypergeometric distribution to determine whether or not a sample group is undergoing coordinated gene expression changes. Statistically, it is asking what is the probability of getting 10 apoptosis genes in my list of 100 differentially expressed genes, given the number of apoptosis genes in the whole genome.
A word of warning about defining your background gene list
Here we use the whole genome to generate our background gene list. This approach may not always be appropriate, and it is up to you to decide what is the most biologically appropriate gene list to use as the background. For example, if you are comparing gene expression across samples from a tissue type that is known to be naturally high in apoptotic procesess, then using the whole genome as a background gene list may ‘drown out’ the signal. Instead, you may want to use all of the genes that are generally known to be expressed in that tissue as your background (e.g., this could be by using a set of control samples, or by using all genes detected as differentially expressed in all your samples together).
We can use the Fisher’s Exact test and the hypergeometric distribution to ask whether being involved in apoptosis is independent of being significantly differentially expressed.
What does this look like in practice? We can use a 2 x 2 table, where the four data points are (from top left to bottom right): the number of genes that are both in the category and are differentially expressed, the number of differentially expressed genes not in the category, the number of category genes that are not differentially expressed, and the number of genes that are neither in the category nor differentially expressed.
Let’s attach some real numbers to those categories and visualise them. We have 10 apoptosis genes in our 100 differentially expressed genes, there are 500 genes annotated as apoptosis, and a total of 10,000 genes in the genome.
Quickly, where did these numbers come from?
Top left: 10 differentially expressed genes which are apoptosis related.
Top right: 90 genes that are differentially expressed but are not apoptosis related.
Bottom left: 490 genes that are apoptosis related but not differentially expressed.
Bottom right: 9410 genes that are neither differentially expressed or apoptosis related.
Row 1 sums to 100 (differentially expressed genes)
Row 2 sums to 9,900 (non-differentially expressed genes)
Col 1 sums to 500 (apoptosis genes)
Col 2 sums to 9,500 (non-apoptosis genes)
Now we can test for an association - or, more accurately, we can test the null hypothesis that being in the apoptosis category is independent of being in the differentially expressed category.
fisher.test(hypergeoDistMatrix)
Fisher's Exact Test for Count Data
data: hypergeoDistMatrix
p-value = 0.03328
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.9832142 4.1416491
sample estimates:
odds ratio
2.133664 The p-value of 0.033 let’s us reject the null hypothesis and we conclude there is an association. In our hypothetical experiment, apoptosis-related genes are over-represented (enriched) in our differentially expressed gene list.
We don’t want to perform this test manually for every possible functional category (adjusting for multiple testing as we go), so we will look to existing software packages to do this for us. But first, let’s look at some caveats and considerations.
Caveats
The hypergeometric distribution and Fisher’s test makes the assumption that all genes are equal. As an exercise, take 30 seconds to think about the ways that genes might not be equal, about RNA-sequencing, and what impact this could have on our test.
One major difference when looking across genes is gene length. RNA-seq tends to produce greater counts for longer genes, which is itself worth remembering, and this gives these genes a greater chance of being statistically differentially expressed. An additional problem is that some gene sets tend to be made up of primarily long genes (because certain processes or functions may require large, complex proteins). Because long genes have a greater chance of being differentially expressed, gene sets which are made up of long genes will have a greater chance of being enriched or over-represented.
We must carry out some type of correction for gene length when conducting over-representation analyses.
GOseq
We will begin with a method called GOseq, which was developed as part of one of the publications initially identifying the gene length bias. An illustration from the full paper is below, demonstrating the relationship between differential gene expression and gene length and between differential expression and number of reads. 
There is a strong positive correlation between differential gene expression and both total gene length and number of reads. Unless corrected for, gene length will bias your over-representation analysis.
Load the goseq package, the dplyr package (if you haven’t already), and use the supportedOrganisms() function to check whether your organism of interest is supported. The code for Saccharomyces cerevisiae is sacCer.
# BiocManager::install("goseq")
library(dplyr)
library(goseq)
supportedOrganisms() %>% View()We need to create an object that will pass GOseq a list of all our genes, with the differentially expressed genes highlighted. To do this, we will use an ifelse() statement which will check every gene in our experiment and ask test if it has an adjusted p-value of less than 0.05. If adj-p-value is < 0.05, return a “1” and if greater than (or equal to) 0.05, return a “0”. After that, add the gene names to this vector of 1s and 0s.
Note: the tt object was created when we used limma to identify differentially expressed genes. tt is 7,127 rows, and has genes stored as row names, includes the logFC, and adj.P.Val columns.
load("tt.RData")genes <- ifelse(tt$adj.P.Val <= 0.05 & tt$logFC > 1, 1, 0)
names(genes) <- rownames(tt)
head(genes) YOR161C YAL038W YML128C YMR105C YHL021C YMR122W-A
0 0 0 0 0 0 table(genes)genes
0 1
6019 1108 Note: because of how we selected our thresholds for differentially expressed or not, genes with a 1 (i.e., a 1 indicates yes) are scattered throughout our gene vector. Using table() is a good way to see how many genes fall into each category.
This is a useful reminder about the conditions of this experiment: of the 7127 genes in the tt object, 1108 are upregulated. This is 15.5 % of all genes. Whether or not this level of differential expression is normal for your samples and experimental conditions is up to you to investigate and make an informed decision on.
Methodology
What is GOSeq doing? How does it correct for gene length?
Remember that we will still use the hypergeometric distribution and Fisher’s Exact test with the 2 x 2 table. Because we know that longer genes are more likely to in the differentially expressed category, we can think of this category as having more weight than it should. GOSeq will calculate a value for each gene which will offset this artificial weight. In other words, if our 2 x 2 table would have 10 genes in the top left square, and some of those genes are very long, GOSeq will treat the number as something slightly less than 10. By treating the value as less than 10, we have taken into account the fact that there shouldn’t really have been 10 genes in there in the first place if not for gene length bias.
Calculate the weighting that should be assigned to each gene with the Probability Weighting Function (the nullp() function).
- Specify the object “gene” which is our list of all genes and whether they are differentially expressed or not, as a binary named vector
- Specify the genome (“sacCer1” for our yeast genome)
- Specify the gene ID type (in this case, we are using ensembl Gene IDs)
This will create a plot in which genes are placed into “bins” based on length, and then gene length vs proportion of differentially expressed is plotted. We can then inspect the pwf object.
Note: you will frequently get a Warning in pcls(G): initial point very close to some inequality constraints message when running
nullp(). You can ignore this.
pwf <- nullp(genes, "sacCer1", "ensGene")Loading sacCer1 length data...
pwf %>% head() DEgenes bias.data pwf
YOR161C 0 1621 0.07861976
YAL038W 0 1504 0.08455472
YML128C 0 1543 0.08231682
YMR105C 0 1711 0.07549338
YHL021C 0 1399 0.09203172
YMR122W-A 0 256 0.26171180pwf %>% tail() DEgenes bias.data pwf
YOR122C 0 383 0.24247801
YJR139C 0 1081 0.12913044
YEL046C 0 1165 0.11738449
YDL084W 0 1342 0.09706676
YOR326W 0 4726 0.06160229
YOL120C 0 563 0.21406738Here we can see that genes have been given a different pwf value based on the bias.data column. The pwf value will be less than 1, and indicates how much a gene should be ‘counted for’ if it is in the differentially expressed category in the 2 x 2 table.
We can predict that pwf and bias (length) are opposing values, and we can visualise that data:
par(mfrow=c(1,2))
hist(pwf$bias.data,30)
hist(pwf$pwf,30)
We have a large number of genes which have low/no bias, and a correspondingly high number of genes with the max pwf weighting. A smaller number of genes have a greater bias and are given lower weighting.
We can now carry out our Fishers Exact test using the pwf value instead of the raw counts. We will use the goseq() function to perform this test, and output the over-representation data into an object.
# BiocManager::install("org.Sc.sgd.db")
# BiocManager::install("AnnotationDbi")
library(org.Sc.sgd.db)
GO.wall <- goseq(pwf, "sacCer1", "ensGene")
GO.wall %>% head() category over_represented_pvalue under_represented_pvalue numDEInCat
1238 GO:0005730 2.056091e-66 1 164
4691 GO:0042254 9.898069e-54 1 187
3463 GO:0022613 4.530403e-51 1 202
3626 GO:0030684 5.381016e-51 1 106
1542 GO:0006364 2.209670e-50 1 150
2879 GO:0016072 5.476986e-49 1 164
numInCat term ontology
1238 300 nucleolus CC
4691 406 ribosome biogenesis BP
3463 485 ribonucleoprotein complex biogenesis BP
3626 164 preribosome CC
1542 293 rRNA processing BP
2879 355 rRNA metabolic process BPWe will want to filter for only those categories with an adjusted p-value < 0.05, and save that information to a new object for easy browsing.
You might have noticed that the GO.wall object doesn’t automatically perform an adjustment for multiple testing, so we will use the p.adjust() function to generate corrected p-values in a new columns for both under- and over- represented P-values.
We will then use the adjusted p-values as a filtering criteria, and include the columns that correspond to category, term, and ontology (using the colnames() function, we can see that these are columns 1, 6 and 7).
GO.wall.padj <- GO.wall |>
mutate(padj.over = p.adjust(over_represented_pvalue, method="fdr")) |>
mutate(padj.under = p.adjust(under_represented_pvalue, method="fdr"))
GO.wall.sig <- GO.wall.padj[GO.wall.padj$padj.over < 0.05 , c(1,6,7)]
GO.wall.sig %>% dim()[1] 105 3GO.wall.sig %>% head() category term ontology
1238 GO:0005730 nucleolus CC
4691 GO:0042254 ribosome biogenesis BP
3463 GO:0022613 ribonucleoprotein complex biogenesis BP
3626 GO:0030684 preribosome CC
1542 GO:0006364 rRNA processing BP
2879 GO:0016072 rRNA metabolic process BPAn optional filtering step can be applied here, which is to remove categories that have a large number of genes. Categories which have a large number of genes tend to be very broad terms, which are not very informative e.g., the categories “organelle”, “biological_process”, and “protein localization”.
If you choose to apply this filter, re-make the GO.wall.sig object. Here we are filtering with a requirement that the category contains fewer than 500 genes.
GO.wall.sig <- GO.wall.padj[GO.wall.padj$padj.over < 0.05 & GO.wall$numInCat < 500, c(1,6,7)]
GO.wall.sig %>% dim()[1] 78 3GO.wall.sig %>% head() category term ontology
1238 GO:0005730 nucleolus CC
4691 GO:0042254 ribosome biogenesis BP
3463 GO:0022613 ribonucleoprotein complex biogenesis BP
3626 GO:0030684 preribosome CC
1542 GO:0006364 rRNA processing BP
2879 GO:0016072 rRNA metabolic process BPThis is a very useful way of getting a broad overview of what processes our differentially expressed genes are involved in.
Plot
We’re going to need a new column Gene Ratio for plotting. The Gene Ratio is the number of DE genes in the GO term divided by the total number of genes in that GO term. This is the proportion of DE genes with that GO term e.g., a Gene ratio of 0.5 means that 50% of the genes in that GO term are DE.
GO.wall.padj <- GO.wall.padj |>
mutate(GeneRatio = numDEInCat / numInCat)Now let’s plot the top 10 over-represented GO terms for upregulated genes in MT. We will use the ggplot2 package to create a bubble plot, where the x-axis is the Gene Ratio, the y-axis is the GO term, the size of the bubble is the number of DE genes in that GO term, and the colour of the bubble is the -log10 of the over-represented p-value.
library(stringr)
library(ggplot2)
GO.wall.padj |>
filter(!is.na(ontology)) |>
slice_min(padj.over, n = 10) |>
ggplot(aes(x = GeneRatio,
y = reorder(str_wrap(term, width = 35), GeneRatio),
size = numDEInCat,
colour = -log10(over_represented_pvalue))) +
geom_point() +
scale_colour_gradient(low = "steelblue", high = "firebrick",
name = "-log10(p-value)") +
scale_size_continuous(name = "DE gene count", range = c(2, 10)) +
labs(title = str_wrap("Over-represented GO terms for upregulated genes in MT", width = 45),
x = "Gene Ratio", y = NULL) +
theme_bw()
Retrieving more information on GO terms
We can use the GO.db package to retrieve more detailed information about the categories we have identified. GO.db will use the unique identifier in the category column and provide more information.
library(GO.db)
GOTERM[[GO.wall$category[1]]]GOID: GO:0005730
Term: nucleolus
Ontology: CC
Definition: A small, dense body one or more of which are present in the
nucleus of eukaryotic cells. It is rich in RNA and protein, is not
bounded by a limiting membrane, and is not seen during mitosis. Its
prime function is the transcription of the nucleolar DNA into 45S
ribosomal-precursor RNA, the processing of this RNA into 5.8S, 18S,
and 28S components of ribosomal RNA, and the association of these
components with 5S RNA and proteins synthesized outside the
nucleolus. This association results in the formation of
ribonucleoprotein precursors; these pass into the cytoplasm and
mature into the 40S and 60S subunits of the ribosome.Alternative approaches to functional analysis and next steps
Further processing and filtering
Our output above contains both over-represented and under-represented GO terms for our upregulated genes from our MT vs WT pairwise comparison, which need to be considered separately.
For each pairwise comparison, we should end up with 4 objects or files:
- Upregulated genes - over-represented GO terms
- Upregulated genes - under-represented GO terms
- Downregulated genes - over-represented GO terms
- Downregulated genes - under-represented GO terms
As you can imagine, this can end up being a lot of files if you have many pairwise comparisons!
You need to consider how you want to further analyse and plot your GO with this in mind.
Investigate the under-represented GO terms! Underrepresented GO terms can be just as interesting as overrepresented GO terms. It can give you a good idea of what functions are depleted in your samples, compared to your background.
Gene set enrichment analysis (GSEA)
GSEA takes an alternative approach to the previous section in that it does not take as input a list of differentially expressed genes. GSEA moves away from the binary “differentially expressed or not” approach and considers the state, or position, of all genes in the list.
This approach has the advantage of not relying on our somewhat arbitrary cutoff thresholds (e.g., a gene with a p-value 0.0501 probably should still contribute to what we think is happening to a pathway).
GO.nobias <- goseq(pwf, "sacCer1", "ensGene", method="Hypergeometric")ClusterProfiler for plotting
ClusterProfiler is a widely used package with the support to create excellent figures (check these out).
Check out the note below on how to create an enrichment plot using semantic similarity.
GO enrichment plot
# BiocManager::install(c("clusterProfiler", "enrichplot", "GO.db"))
library(clusterProfiler)
library(enrichplot)
library(org.Sc.sgd.db)
library(GO.db)
library(dplyr)
library(ggplot2)
# Step 1: define the DE genes and background
# These come from the same objects we used for goseq:
# - background = every gene in the tt object
# - DE_genes_UP = the genes flagged as 1 in our named binary vector
background <- names(genes)
DE_genes_UP <- names(genes)[genes == 1]
length(background)[1] 7127length(DE_genes_UP)[1] 1108# Step 2: build TERM2GENE from the yeast annotation package
# enricher() needs a two-column dataframe: GO term first, gene ID second.
# We use GOALL (rather than GO) so each gene is also annotated to the parent
# terms of its direct GO terms, which is what over-representation needs.
# We keep only Biological Process terms and genes in our background.
term2gene <- AnnotationDbi::select(org.Sc.sgd.db,
keys = background,
keytype = "ORF",
columns = c("GOALL", "ONTOLOGYALL")) |>
filter(ONTOLOGYALL == "BP", !is.na(GOALL)) |>
dplyr::select(GOALL, ORF) |>
distinct()'select()' returned 1:many mapping between keys and columnshead(term2gene) GOALL ORF
1 GO:0006810 YOR161C
2 GO:0008150 YOR161C
3 GO:0009987 YOR161C
4 GO:0051179 YOR161C
5 GO:0051234 YOR161C
6 GO:0055085 YOR161C# Step 3: generate term2name
# We need the GO ID and its name (e.g., GO:0000001 and 'mitochondrion inheritance')
# so the plot has readable labels. We only keep terms present in term2gene.
term2name <- AnnotationDbi::select(GO.db,
keys = unique(term2gene$GOALL),
columns = c("GOID", "TERM")) |>
dplyr::select(GOID, TERM)'select()' returned 1:1 mapping between keys and columnshead(term2name) GOID TERM
1 GO:0006810 transport
2 GO:0008150 biological_process
3 GO:0009987 cellular process
4 GO:0051179 localization
5 GO:0051234 establishment of localization
6 GO:0055085 transmembrane transport# Step 4: run enricher()
# universe = the background gene list (same as goseq)
# gene = the DE gene list (a character vector of IDs, not the named binary vector)
# The gene set size limits match the <500 genes filter used earlier.
cp_result <- enricher(
gene = DE_genes_UP,
universe = background,
TERM2GENE = term2gene,
TERM2NAME = term2name,
pvalueCutoff = 0.05,
qvalueCutoff = 0.05,
pAdjustMethod = "BH",
minGSSize = 10,
maxGSSize = 500
)# Step 5: calculate semantic similarity
# emapplot requires the pairwise similarity between terms to be calculated first.
# This creates a similarity matrix.
cp_result_sim <- pairwise_termsim(cp_result)
# Step 6: visualisation
emapplot(cp_result_sim,
showCategory = 10, # change how many GOs on plot
size_category = 1, # change size of bubbles
node_label_size = 2) + # change GO font size
labs(title = "GO term enrichment map: upregulated genes in MT vs WT")
Additional resources
- The Carpentries also provide a guide in an incubator (in development) lesson on Gene Set Enrichment Analysis.