RNA-seq analysis in R

Gene Set Testing for RNA-seq

Original Authors: Belinda Phipson, Anna Trigos, Matt Ritchie, Maria Doyle, Harriet Dashnow, Charity Law, Stephane Ballereau, Oscar Rueda, Ashley Sawle Based on the course RNAseq analysis in R delivered on May 11/12th 2016 and modified by Cancer Research Uk Cambridge Centre for the Functional Genomics Autumn School 2017

Resources and data files

This material has been created using the following resources:
http://www.statsci.org/smyth/pubs/QLedgeRPreprint.pdf [@Lun2016]
http://monashbioinformaticsplatform.github.io/RNAseq-DE-analysis-with-R/99-RNAseq_DE_analysis_with_R.html
http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html https://bioconductor.github.io/BiocWorkshops/rna-seq-data-analysis-with-deseq2.html

This section also uses code from Stephen Turner’s guide to fgsea https://stephenturner.github.io/deseq-to-fgsea/

Gene Set Testing

In the early days of microarray analysis, people were happy if they got a handful of differentially-expressed genes that they could validate or follow-up. However, with later technologies (and depending on the experimental setup) we might have thousands of statistically-significant results, which no-one has the time to follow-up. Also, we might be interested in pathways / mechanisms that are altered and not just individual genes.

In this section we move towards discovering if our results are biologically significant. Are the genes that we have picked statistical flukes, or are there some commonalities.

There are two different approaches one might use, and we will cover the theory behind both

• There is also a bunch of websites for doing the tests
  • we will show how they are done in Bioconductor so the theory is clear
• We will assume we have done a differential-expression analysis, but the same techniques can be used for other situations when we have a gene list
  • ChIP-seq
  • RNA-seq

fgsea analysis

The fgsea package is a free implementation of the Broad’s GSEA software and is described in more detail in the package vignette “fast preranked gene set enrichment analysis (GSEA)”:

The GSEA analysis is performed by:

• 1. ranking all genes in the data set based on their correlation to the chosen phenotype,
• 2. identifying the rank positions of all members of the gene set, and
• 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.

“After establishing the ES for each gene set across the phenotype, GSEA reiteratively randomizes the sample labels and retests for enrichment across the random classes. By performing repeated class label randomizations, the ES for each gene set across the true classes can be compared to the ES distribution from the random classes. Those gene sets that significantly outperform iterative random class permutations are considered significant.”

The article describing the original software is available here and there is also a commentary on GSEA.

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. These gene sets use Entrez ID as their identifier.

library(fgsea)

Loading required package: Rcpp

load("Robjects/LvsV_annotated.rdata")
head(results.annotated)

             ENSEMBL   baseMean log2FoldChange     lfcSE     stat        pvalue
1 ENSMUSG00000000381  398943.44       8.722325 0.3815879 22.94608 1.612341e-116
2 ENSMUSG00000061937 1134857.22       7.351405 0.3761569 19.72683  1.268851e-86
3 ENSMUSG00000026417   11320.43       6.078220 0.3332017 18.13304  1.748099e-73
4 ENSMUSG00000022491  514850.04       7.765956 0.4539439 17.96995  3.350319e-72
5 ENSMUSG00000061388   30555.67       8.282907 0.4740455 17.36398  1.546169e-67
6 ENSMUSG00000024903   60039.52       7.312158 0.4522371 16.27868  1.398616e-59
           padj  SYMBOL                                         GENENAME ENTREZID
1 2.785480e-112     Wap                              whey acidic protein    22373
2  1.096033e-82 Csn1s2a                           casein alpha s2-like A    12993
3  1.006672e-69    Pigr                polymeric immunoglobulin receptor    18703
4  1.447003e-68 Glycam1 glycosylation dependent cell adhesion molecule 1    14663
5  5.342324e-64 Csn1s2b                           casein alpha s2-like B    12992
6  4.027083e-56    Lao1                           L-amino acid oxidase 1   100470

An appealing feature of the method is that it does not require us to impose arbitrary cut-offs on the dataset to decide what is differentially-expressed or not. The steps in producing the input required for GSEA are i) retrieving the ranked statistics ii) naming each one according to Entrez ID.

library(dplyr)
gseaInput <- filter(results.annotated, !is.na(ENTREZID)) %>% 
  arrange(stat)
ranks <- pull(gseaInput,stat)
names(ranks) <- gseaInput$ENTREZID

The Walter+Eliza Hall Institutes Bioinformatics service have made mouse versions of the MSigDB datasets available for download. This should already be available in the Robjects folder.

load("Robjects/mouse_H_v5p2.rdata")
pathways <- Mm.H

The analysis is one call to the fgsea function. We can automatically exclude any pathways with too many or too few genes.

library(fgsea)
fgseaRes <- fgsea(pathways, ranks, minSize=15, maxSize = 500, nperm=1000)

There are duplicate gene names, fgsea may produce unexpected results

dim(fgseaRes)

[1] 50  8

#head(fgseaRes)

The results table gives the names of each pathway that was tested and the stats from doing the test. We can make this into a “tidy” object with the following code.

fgseaResTidy <- fgseaRes %>%
  as_tibble() %>%
  arrange(desc(NES))
# Show in a nice table:
fgseaResTidy 

library(ggplot2)
ggplot(fgseaResTidy, aes(reorder(pathway, NES), NES)) +
  geom_col(aes(fill=padj<0.05)) +
  coord_flip() +
  labs(x="Pathway", y="Normalized Enrichment Score",
       title="Hallmark pathways NES from GSEA")

The enrichment plot will show where the genes belonging to a particular gene set are towards the top or the bottom of the genelist, and how the enrichment score is calculated across the dataset.

Here we show the enrichment plot for the pathway with the most positive enrichment score.

plotEnrichment(pathways[["HALLMARK_OXIDATIVE_PHOSPHORYLATION"]],
               ranks)

Challenge 1

1. What pathway has the most extreme negative enrichment score? How can you identify this pathway from the results table?
2. Make the enrichment plot for this pathway

Visualising pathways using a heatmap

The names of genes involved in particular pathways can be extracted from the pathways object. We can then use these genes to make a heatmap, or other visualisation, to show how these genes partition the dataset.

### Load the dds object if not present
load("Robjects/preprocessing.Rdata")

Loading required package: DESeq2
Loading required package: S4Vectors
Loading required package: stats4
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: ‘BiocGenerics’

The following objects are masked from ‘package:parallel’:

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ, clusterExport,
    clusterMap, parApply, parCapply, parLapply, parLapplyLB, parRapply, parSapply,
    parSapplyLB

The following objects are masked from ‘package:dplyr’:

    combine, intersect, setdiff, union

The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs

The following objects are masked from ‘package:base’:

    anyDuplicated, append, as.data.frame, basename, cbind, colnames, dirname,
    do.call, duplicated, eval, evalq, Filter, Find, get, grep, grepl, intersect,
    is.unsorted, lapply, Map, mapply, match, mget, order, paste, pmax, pmax.int,
    pmin, pmin.int, Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
    table, tapply, union, unique, unsplit, which, which.max, which.min


Attaching package: ‘S4Vectors’

The following objects are masked from ‘package:dplyr’:

    first, rename

The following object is masked from ‘package:base’:

    expand.grid

Loading required package: IRanges

Attaching package: ‘IRanges’

The following objects are masked from ‘package:dplyr’:

    collapse, desc, slice

Loading required package: GenomicRanges
Loading required package: GenomeInfoDb
Loading required package: SummarizedExperiment
Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with 'browseVignettes()'. To cite
    Bioconductor, see 'citation("Biobase")', and for packages 'citation("pkgname")'.

Loading required package: DelayedArray
Loading required package: matrixStats

Attaching package: ‘matrixStats’

The following objects are masked from ‘package:Biobase’:

    anyMissing, rowMedians

The following object is masked from ‘package:dplyr’:

    count

Loading required package: BiocParallel

Attaching package: ‘DelayedArray’

The following objects are masked from ‘package:matrixStats’:

    colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges

The following objects are masked from ‘package:base’:

    aperm, apply, rowsum

library(pheatmap)
my_genes <- filter(results.annotated, ENTREZID %in% pathways[["HALLMARK_MYC_TARGETS_V2"]]) %>% 
  pull(ENSEMBL)
vsd <- vst(dds)

using 'avgTxLength' from assays(dds), correcting for library size

mat <- assay(vsd)[my_genes,]
mat <- mat - rowMeans(mat)
dim(mat)

[1] 65 12

rownames(sampleinfo) <- sampleinfo$run
pheatmap(mat,
         annotation_col = sampleinfo[,c("Status","CellType")])

Gene Set Testing - competitive gene set tests

GOseq analysis

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 (GOseq article)

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.”

First, create list of DEGs. We don’t have to be too strict about our criteria for a gene to be differentially-expressed, as too few genes will not give us many enriched pathways.

genes <- results.annotated$padj < 0.05 & !is.na(results.annotated$padj)
names(genes) <- results.annotated$ENSEMBL

Fit the Probability Weighting Function (PWF):

library(goseq)
if(!require(TxDb.Mmusculus.UCSC.mm10.knownGene)) BiocManager::install("TxDb.Mmusculus.UCSC.mm10.knownGene")
#print(supportedGeneIDs())
#print(supportedGenomes())
pwf <- nullp(genes, "mm10","ensGene")

libraries ‘/usr/local/lib/R/site-library’, ‘/usr/lib/R/site-library’ contain no packagesargument 'pattern' has length > 1 and only the first element will be usedcollapsing to unique 'x' valuescollapsing to unique 'x' valuescollapsing to unique 'x' valuescollapsing to unique 'x' valuescollapsing to unique 'x' valuescollapsing to unique 'x' valuescollapsing to unique 'x' valuescollapsing to unique 'x' valuescollapsing to unique 'x' valuescollapsing to unique 'x' valuescollapsing to unique 'x' valuescollapsing to unique 'x' values

Conduct gene set enrichment analysis:

rr #?goseq goseq_res <- goseq(pwf, 10,,test.cats=:BP) r goseq_res

?goseq

Analysis with clusterProfiler

clusterProfiler is another Bioconductor package for over-representation analysis. It’s main advantage is that it provides some nice visualisation methods.

Firstly, we can identify over-represented GO terms and visualise these as a network.

library(clusterProfiler)
universe <- results.annotated %>% pull(ENTREZID)
sigGenes <- results.annotated %>% 
  filter(padj < 0.05, !is.na(ENTREZID)) %>% pull(ENTREZID)

enrich_go <- enrichGO(
  gene= sigGenes,
  OrgDb = org.Mm.eg.db,
  keyType = "ENTREZID",
  ont = "BP",
  universe = universe,
  qvalueCutoff = 0.05,
  readable=TRUE
)

rr dotplot(enrich_go)

A dot plot can show us the most enriched pathways, and the size of each.

dotplot(enrich_go)

emapplot(enrich_go)

We also perform enrichment of KEGG pathways with the clusterProfiler package. We could do this with goseq, but again clusterProfiler has some nice visualisations and uses online resources rather than relying on Bioconductor annotation.

To run the analysis we can use the same list of genes as before, and need to know the kegg code for our given organism.

search_kegg_organism('mouse', by='common_name')
keg_res <- enrichKEGG(gene=sigGenes, organism="mmu")

head(keg_res,n=10)

If we want to view the network for any pathway, we can use the browseKEGG function which will take us to the relevant page on the KEGG website.

We will choose an example of a smaller pathway to make things a bit clearer.

browseKEGG(keg_res, 'mmu03320')

The pathview package will produce a version of this plot (png) in your working directory.

library(pathview)
logFC <- results.annotated$log2FoldChange
names(logFC) <- results.annotated$ENTREZID

pathview(gene.data = logFC, 
         pathway.id = "mmu03320", 
         species = "mmu", 
         limit = list(gene=5, cpd=1))

Creating Gene lists to use with an online tool

There are also many online tools that one could use to perform a gene set or ontology analysis.

• DAVID
• GeneTrail
• GORilla

The tools generally require your input genes lists to be uploaded as a simple text file. In this final challenge, we will create some files that you might use in one of these tools.

A file containing names of background genes

This file has one column which lists all the gene names present in the analysis. Gene Symbols are commonly used, although a tool may accept Ensembl or Refseq names

A file containing names of significant genes

This file has one column which lists the genes that passed the threshold for statistical significance (e.g. p-value less than 0.05) in your analysis. Gene Symbols are commonly used, although a tool may accept Ensembl or Refseq names

Challenge

Create two text files that can be imported into online tools for further analysis 1. A list of background genes 2. A list of differentially expressed genes 3. Load these files into GOrilla for analysis HINT: the write.table function is able to write a data frame to a txt file in R. You will need to set the appropriate arguments to make sure that a text file with only one column is created.