RNA-seq analysis in R

Gene Set Testing for RNA-seq

The list of differentially expressed genes is sometimes so long that its interpretation becomes cumbersome and time consuming. A common downstream procedure is gene set testing. It aims to understand which pathways or gene networks the differentially expressed genes are implicated in.

Various ways exist to test for enrichment of biological pathways.

GSEA analysis

Gene Set Enrichment Analysis GSEA was tests whether a set of genes of interest, e.g. genes (Subramanian et al. 2005). The software is distributed by the Broad Institute and is freely available for use by academic and non-profit organisations.

In addition to the GSEA software the Broad also provide a number of very well curated gene sets for testing against your data - the Molecular Signatures Database (MSigDB). Unfortunately, these are collections of human genes, however, these lists have been translated to mouse equivalents by the Walter+Eliza Hall Institutes Bioinformatics service and made avaialble for download.

The analysis is performed by:

• 1. ranking all genes in the data set
• 2. identifying the rank positions of all members of the gene set in the ranked data set
• 3. calculating an enrichment score (ES) that represents the difference between the observed rankings and that which would be expected assuming a random rank distribution.

commentary on GSEA. The article describing the original software is available here.

fgsea

The fgsea package (Sergushichev 2016) implements the same algorithm in R vignette “fast preranked gene set enrichment analysis (GSEA)”.

library(fgsea)

load("Robjects/Annotated_Results_LvV.RData")

Ranking Data

We need to provide fgsea a vector containing numeric data by which it should rank the genes. To start with we will simply use a rank based on their fold change. We do not need to rank the genes ourselves, fgsea will do this based on the data we provide.

We need to exclude genes for which we do not have Entrez IDs. Also, we should use the shrunk LFC values.

gseaDat <- filter(shrinkLvV, !is.na(Entrez))
rankData <- gseaDat$logFC
names(rankData) <- gseaDat$Entrez
head(rankData)

    497097      19888      20671      27395      18777      21399 
 0.5851914 -1.0498457 -0.5893634  0.4313540  0.6478714 -0.1110589 

Load pathways

load("Robjects/mouse_H_v5.RData")
pathwaysH <- Mm.H

Conduct analysis

fgseaRes <- fgsea(pathwaysH, 
                  rankData, 
                  minSize=15, 
                  maxSize = 500, 
                  nperm=1000)

There are ties in the preranked stats (0.05% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.

The warning produced indicates that there are few genes that have the same fold change and so are ranked equally. fgsea with arbitrarily order determine which comes first in the ranked list. As long as this number is small it shouldn’t significantly effect the results. If the number is large something is suspicious about the fold change results.

Lets look at the top 10 results.

fgseaRes %>% 
    arrange(desc(abs(NES))) %>% 
    top_n(10, -padj)

Enrichment score plot

plotEnrichment(pathwaysH[["HALLMARK_ESTROGEN_RESPONSE_EARLY"]], rankData)

Remember to check the GSEA article for the complete explanation.

GSEA table plot

The function plotGseaTable allows us to plot a summary figue showing the results for multiple pathways.

topPathways <- fgseaRes %>% 
    top_n(20, wt=-padj) %>% 
    arrange(-NES) %>% 
    pull(pathway)
plotGseaTable(pathwaysH[topPathways], 
              rankData, 
              fgseaRes, 
              gseaParam = 0.5)

Challenge 1

Another common way to rank the genes is to order by pvalue, but also, sorting so that upregulated genes are at start and downregulated at the other - you can do this combining the sign of the fold change and the pvalue.
1. Rank the genes by statisical significance - you will need to create a new ranking value using -log10({p value}) * sign({Fold Change})
2. Load the “C2” pathways from the the data/mouse_c2_v5.RData file
3. Run fgsea using the new ranked genes and the C2 pathways
4. Run fgsea using the new ranked genes and the H pathways. How do these results differ from the ones we got when ranking by the fold change alone?

GO enrichment analysis

goseq

GOseq is a method to conduct Gene Ontology (GO) analysis suitable for RNA-seq data as it accounts for the gene length bias in detection of over-representation (Young et al. 2010).

From the GOseq vignette:

• GOseq first needs to quantify the length bias present in the dataset under consideration.
• This is done by calculating a Probability Weighting Function or PWF which can be thought of as a function which gives the probability that a gene will be differentially expressed (DE), based on its length alone.
• The PWF is calculated by fitting a monotonic spline to the binary data series of differential expression (1=DE, 0=Not DE) as a function of gene length.
• The PWF is used to weight the chance of selecting each gene when forming a null distribution for GO category membership.
• The fact that the PWF is calculated directly from the dataset under consideration makes this approach robust, only correcting for the length bias present in the data.

“GO analysis of RNA-seq data requires the use of random sampling in order to generate a suitable null distribution for GO category membership and calculate each category’s significance for over representation amongst DE genes. … In most cases, the Wallenius distribution can be used to approximate the true null distribution, without any significant loss in accuracy. The goseq package implements this approximation as its default option.”

library(goseq)
supportedOrganisms() %>% filter(str_detect(Genome, "mm"))

Create a list of differentially expressed genes

The input for goseq is a vector that indicates, for each gene, whether or not it is significantly differentially expressed. This should be a named vector, where the names are the gene ids and the values are 1 if the gene is significant and a 0 if it is not.

In this case we can use the Ensembl gene IDs.

sigData <- as.integer( shrinkLvV$FDR < 0.01 & !is.na(shrinkLvV$FDR) )
names(sigData) <- shrinkLvV$GeneID

Fit the Probability Weighting Function (PWF)

pwf <- nullp(sigData, "mm10", "ensGene", bias.data = shrinkLvV$medianTxLength)

initial point very close to some inequality constraints

Conduct GO enrichment analysis

goResults <- goseq(pwf, "mm10","ensGene", test.cats=c("GO:BP"))

package ‘AnnotationDbi’ was built under R version 3.5.1

Plot the top 10

goResults %>% 
    top_n(10, wt=-over_represented_pvalue) %>% 
    mutate(hitsPerc=numDEInCat*100/numInCat) %>% 
    ggplot(aes(x=hitsPerc, 
               y=term, 
               colour=over_represented_pvalue, 
               size=numDEInCat)) +
        geom_point() +
        expand_limits(x=0) +
        labs(x="Hits (%)", y="GO term", colour="p value", size="Count")

Get the GO information for the GO accessions

library(GO.db)
GOTERM[[goResults$category[1]]]

GOID: GO:0044281
Term: small molecule metabolic process
Ontology: BP
Definition: The chemical reactions and pathways involving small molecules,
    any low molecular weight, monomeric, non-encoded molecule.
Synonym: small molecule metabolism

Challenge 2

1. Create a vector showing genes that are statistically significant at FDR < 0.01 and that are up-regulated by at least 4x (logFC>2) in lactating mice
2. Run a goseq analysis on this gene list
3. Plot the results
4. How is this result different to the previous GO analysis?

KEGG pathway enrichment analysis

clusterProfiler

We can analyse for enrichment of KEGG pathways in much the same way as for GO terms. We could also use goseq for this, but this time we’re going to use clusterProfiler (Yu et al. 2012). clusterprofiler supports direct online access of the current KEGG database, rather than relying on R annotation packages, it also provides some nice visualisation options.

library(clusterProfiler)

package ‘clusterProfiler’ was built under R version 3.5.1

search_kegg_organism('mmu', by='kegg_code')

KEGG enrichment analysis

The input for the KEGG enrichment is list of gene IDs for significant genes.

This time we’ll only use gene with an absolute fold change greater than 2.

For this tool we need to use Entrez IDs, so once again we will have to eliminate the missing values.

sigGenes <- shrinkLvV %>% 
    filter(FDR < 0.05 & !is.na(FDR) & 
               abs(logFC) > 1 & 
               !is.na(Entrez)) %>% 
    pull(Entrez)
kk <- enrichKEGG(gene = sigGenes, organism = 'mmu')
head(kk, n=10)

Visualise a pathway

In a browser

clusterProfile has a function browseKegg that allows you to view the Kegg pathway with the genes colours in in your browser.

browseKEGG(kk, 'mmu03320')

As a file

The package pathview (Luo et al. 2013) can be used to generate figures of KEGG pathways.

One advantage over clusterProfiles browser method is that the genes are coloured according to their fold change levels in our data. To do this we need to pass pathview a named vector of fold change data (actually you could colour by any numeric vector, e.g. p-value).

The package plots the KEGG pathway to a png file in the working directory.

library(pathview)

package ‘pathview’ was built under R version 3.5.1Loading required package: org.Hs.eg.db

##############################################################################
Pathview is an open source software package distributed under GNU General
Public License version 3 (GPLv3). Details of GPLv3 is available at
http://www.gnu.org/licenses/gpl-3.0.html. Particullary, users are required to
formally cite the original Pathview paper (not just mention it) in publications
or products. For details, do citation("pathview") within R.

The pathview downloads and uses KEGG data. Non-academic uses may require a KEGG
license agreement (details at http://www.kegg.jp/kegg/legal.html).
##############################################################################

logFC <- annotLvV$logFC
names(logFC) <- annotLvV$Entrez
pathview(gene.data = logFC, 
         pathway.id = "mmu03320", 
         species = "mmu", 
         limit = list(gene=5, cpd=1))

Info: Downloading xml files for mmu03320, 1/1 pathways..
Info: Downloading png files for mmu03320, 1/1 pathways..
'select()' returned 1:1 mapping between keys and columns
Info: Working in directory /Users/sawle01/Documents/RNAseq_course/RNAseq_March_2019/Course_Materials
Info: Writing image file mmu03320.pathview.png

mmu03320.pathview.png:

mmu03320 - PPAR signaling pathway

Challenge 3

1. Use pathview to export a figure for “mmu04060”, but this time only use genes that are statistically significant at FDR < 0.01

---

References

Luo, Weijun, Brouwer, and Cory. 2013. “Pathview: An R/Bioconductor Package for Pathway-Based Data Integration and Visualization.” Bioinformatics 29 (14): 1830–1. doi:10.1093/bioinformatics/btt285.

Sergushichev, Alexey. 2016. “An Algorithm for Fast Preranked Gene Set Enrichment Analysis Using Cumulative Statistic Calculation.” bioRxiv. Cold Spring Harbor Labs Journals. doi:10.1101/060012.

Subramanian, Aravind, Pablo Tamayo, Vamsi K. Mootha, Sayan Mukherjee, Benjamin L. Ebert, Michael A. Gillette, Amanda Paulovich, et al. 2005. “Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-Wide Expression Profiles.” Proceedings of the National Academy of Sciences 102 (43). National Academy of Sciences: 15545–50. doi:10.1073/pnas.0506580102.

Young, Matthew D, Matthew J Wakefield, Gordon K Smyth, and Alicia Oshlack. 2010. “Gene Ontology Analysis for Rna-Seq: Accounting for Selection Bias.” Genome Biology 11: R14.

Yu, Guangchuang, Li-Gen Wang, Yanyan Han, and Qing-Yu He. 2012. “ClusterProfiler: An R Package for Comparing Biological Themes Among Gene Clusters.” OMICS: A Journal of Integrative Biology 16 (5): 284–87. doi:10.1089/omi.2011.0118.