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 Based on the course RNAseq analysis in R delivered on May 11/12th 2016 ## Resources and data files

This material has been created using the following resources:
http://www.statsci.org/smyth/pubs/QLedgeRPreprint.pdf (Lun, Chen, and Smyth 2016)
https://github.com/MonashBioinformaticsPlatform/RNAseq-DE-analysis-with-R

Data files:
https://figshare.com/s/1d788fd384d33e913a2a : mouse_c2_v5.rdata mouse_H_v5.rdata

Gene Set Testing

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. Two types of test answer two different questions: competitive and self-contained gene set tests.

Competitive gene set tests, such as those implemented in GOseq and camera, ask the question whether the differentially expressed genes tend to be over-represented in the gene set, compared to all the other genes in the experiment.

Self-contained tests, which include the ROAST procedure, ask the question “Are the genes in the set/pathway differentially expressed as a whole?”

load("Robjects/DE.Rdata")

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

Create list of DEGs:

# Retrieve list of all genes tested:
results <- as.data.frame(topTags(lrt_BvsL, n = Inf))
print(head(results))

# Derive list of DEGs by filtering on FDR:
genes <- results$FDR < 0.01
# Add gene names to that list:
names(genes) <- rownames(results)
print(head(genes))

110308  50916  12293  56069  24117  12818 
  TRUE   TRUE   TRUE   TRUE   TRUE   TRUE 

library(goseq)
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
#print(supportedGeneIDs())
#print(supportedGenomes())

Fit the Probability Weighting Function (PWF):

pwf <- nullp(genes, "mm10","knownGene")

Can't find mm10/knownGene length data in genLenDataBase...
Found the annotaion package, TxDb.Mmusculus.UCSC.mm10.knownGene
Trying to get the gene lengths from it.
initial point very close to some inequality constraints

Conduct gene set enrichment analysis:

#?goseq
go.results <- goseq(pwf, "mm10","knownGene")

Fetching GO annotations...

For 1098 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

go.results

fgsea analysis

From the fgsea vignette “fast preranked gene set enrichment analysis (GSEA)”:

This 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.” commentary on GSEA. The article describing the original software is available here.

library(fgsea)

Loading required package: Rcpp

Create ranks

Rank all genes based on their fold change.

results.ord <- results[ order(-results[,"logFC"]), ]
head(results.ord)

ranks <- results.ord$logFC
names(ranks) <- rownames(results.ord)
head(ranks)

  791407   382036   192200    21990   246133    12310 
8.197169 7.947515 7.894963 7.794611 7.708981 7.660325 

Plot the ranked fold changes.

#plot(ranks)
barplot(ranks)

Load pathways

load("data/mouse_H_v5.rdata")
pathways <- Mm.H

Conduct analysis

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

[1] "data.table" "data.frame"

dim(fgseaRes)

[1] 50  8

#head(fgseaRes)

Table of results

head(fgseaRes[order(padj), ])

Plot outcome for the ‘HALLMARK_MYOGENESIS’ pathway

First find rank of the ‘HALLMARK_MYOGENESIS’ pathway genes in the sorted genes. Then plot as before.

# We will create a barplot of logFC for the sorted genes and add one vertical red bar for each gene in the 'HALLMARK_MYOGENESIS' pathway
tmpInd <- match(pathways[["HALLMARK_MYOGENESIS"]],names(ranks))
tmpInd <- tmpInd[!is.na(tmpInd)]
ranks2 <- rep(0,length(ranks))
ranks2[tmpInd] <- ranks[tmpInd]
barplot(ranks2)

Enrichment score plot

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

Remember to check the GSEA article for the complete explanation.

GSEA table plot

Select top pathways and plot outcome for all these

pathwaysUp <- fgseaRes[ES > 0]
topPathwaysUp <- head( pathwaysUp$pathway[ order(pathwaysUp$pval) ], n=10 )
pathwaysDown <- fgseaRes[ES < 0]
topPathwaysDown <- head( pathwaysDown$pathway[ order( pathwaysDown$pval ) ], n=10 )
topPathways <- c(topPathwaysUp, rev(topPathwaysDown))
plotGseaTable(pathways[topPathways], ranks, fgseaRes, gseaParam = 0.5)

?plotGseaTable

Gene Set Testing - self-contained gene set tests

ROAST gene set testing

ROAST is an example of a self-contained gene set test (Wu et al. 2010). It asks the question, “Do the genes in my set tend to be differentially expressed between my conditions of interest?”. ROAST does not use information on the other genes in the experiment, unlike camera (see below). ROAST is a good option for when you’re interested in a specific set, or a few sets. It is not really used to test thousands of sets at one time.

Here we will be using the C2 gene sets for mouse, available as .rdata files from the WEHI bioinformatics page http://bioinf.wehi.edu.au/software/MSigDB/index.html. The C2 gene sets contain 4725 curated gene sets collected from a variety of places: BioCarta, KEGG, Pathway Interaction Database, Reactome as well as some published studies. It doesn’t include GO terms.

# Load in the mouse c2 gene sets
# The R object is called Mm.c2
load("data/mouse_c2_v5.rdata")
# Have a look at the first few gene sets
names(Mm.c2)[1:5]

[1] "KEGG_GLYCOLYSIS_GLUCONEOGENESIS"               "KEGG_CITRATE_CYCLE_TCA_CYCLE"                 
[3] "KEGG_PENTOSE_PHOSPHATE_PATHWAY"                "KEGG_PENTOSE_AND_GLUCURONATE_INTERCONVERSIONS"
[5] "KEGG_FRUCTOSE_AND_MANNOSE_METABOLISM"         

# Number of gene sets in C2
length(Mm.c2)

[1] 4725

The gene identifiers are Entrez Gene ID, as are the rownames of our DGEList object ‘dgeObj’. We need to map the Entrez gene ids between the list of gene sets and our DGEList object. We can do this using the ids2indices function.

c2.ind <- ids2indices(Mm.c2, rownames(dgeObj$counts))

Let’s see if there are any MYC signalling pathways in MsigDB C2 collection. We can do this with the grep command on the names of the gene sets.

grep("MYC_", names(c2.ind))

 [1]  496  544  621 1591 1592 1783 1784 2029 2030 2031 2032 2194 2195 2203 2204 2435 2436 2506 2507 2600 2609 2620 2621 2629 2676
[26] 2694 2706 2716 2727 2744 2786 2824 3311 3401 3817 3943 3944 4204 4205 4206 4207 4433 4434 4616

# Let's save these so that we can subset c2.ind to test all gene sets with MYC in the name
myc <- grep("MYC_",names(c2.ind))
# What are these pathways called?
names(c2.ind)[myc]

 [1] "PID_MYC_ACTIV_PATHWAY"                           "PID_MYC_PATHWAY"                                
 [3] "PID_MYC_REPRESS_PATHWAY"                         "ODONNELL_TARGETS_OF_MYC_AND_TFRC_UP"            
 [5] "ODONNELL_TARGETS_OF_MYC_AND_TFRC_DN"             "CAIRO_PML_TARGETS_BOUND_BY_MYC_UP"              
 [7] "CAIRO_PML_TARGETS_BOUND_BY_MYC_DN"               "SCHLOSSER_MYC_AND_SERUM_RESPONSE_SYNERGY"       
 [9] "SCHLOSSER_MYC_TARGETS_REPRESSED_BY_SERUM"        "SCHLOSSER_MYC_TARGETS_AND_SERUM_RESPONSE_DN"    
[11] "SCHLOSSER_MYC_TARGETS_AND_SERUM_RESPONSE_UP"     "CEBALLOS_TARGETS_OF_TP53_AND_MYC_UP"            
[13] "CEBALLOS_TARGETS_OF_TP53_AND_MYC_DN"             "KIM_MYC_AMPLIFICATION_TARGETS_UP"               
[15] "KIM_MYC_AMPLIFICATION_TARGETS_DN"                "BENPORATH_MYC_TARGETS_WITH_EBOX"                
[17] "BENPORATH_MYC_MAX_TARGETS"                       "MORI_EMU_MYC_LYMPHOMA_BY_ONSET_TIME_UP"         
[19] "MORI_EMU_MYC_LYMPHOMA_BY_ONSET_TIME_DN"          "SCHUHMACHER_MYC_TARGETS_DN"                     
[21] "LEE_LIVER_CANCER_MYC_TGFA_DN"                    "LEE_LIVER_CANCER_MYC_UP"                        
[23] "MENSSEN_MYC_TARGETS"                             "LEE_LIVER_CANCER_MYC_TGFA_UP"                   
[25] "LEE_LIVER_CANCER_MYC_E2F1_DN"                    "LEE_LIVER_CANCER_MYC_E2F1_UP"                   
[27] "SCHUHMACHER_MYC_TARGETS_UP"                      "LEE_LIVER_CANCER_MYC_DN"                        
[29] "COLLER_MYC_TARGETS_UP"                           "COLLER_MYC_TARGETS_DN"                          
[31] "YU_MYC_TARGETS_UP"                               "YU_MYC_TARGETS_DN"                              
[33] "BILD_MYC_ONCOGENIC_SIGNATURE"                    "SANSOM_APC_MYC_TARGETS"                         
[35] "SOUCEK_MYC_TARGETS"                              "ELLWOOD_MYC_TARGETS_UP"                         
[37] "ELLWOOD_MYC_TARGETS_DN"                          "DANG_REGULATED_BY_MYC_UP"                       
[39] "DANG_REGULATED_BY_MYC_DN"                        "DANG_MYC_TARGETS_UP"                            
[41] "DANG_MYC_TARGETS_DN"                             "ACOSTA_PROLIFERATION_INDEPENDENT_MYC_TARGETS_UP"
[43] "ACOSTA_PROLIFERATION_INDEPENDENT_MYC_TARGETS_DN" "ALFANO_MYC_TARGETS"                             

Let’s use ROAST to see if these MYC related gene sets tend to be differentially expressed.

# Specify a design matrix without an intercept term
design <- model.matrix(~ CellType + Status, data = sampleinfo)
# Run ROAST
myc.rst <- roast(dgeObj,
                 index=c2.ind[myc],
                 design=design,
                 contrast=3,
                 nrot=999)
myc.rst[1:15,]

Each row corresponds to a single gene set.

The NGenes column gives the number of genes in each set.

The PropDown and PropUp columns contain the proportions of genes in the set that are down- and up-regulated, respectively, with absolute fold changes greater than 2.

The net direction of change is determined from the significance of changes in each direction, and is shown in the Direction column.

The PValue provides evidence for whether the majority of genes in the set are DE in the specified direction

The PValue.Mixed provides evidence for whether the majority of genes in the set are DE in any direction.

FDRs are computed from the corresponding p-values across all sets.

Challenge 1

1. Test whether the MYC signalling pathways tend to be differentially expressed between basal virgin vs lactating.
2. Look for gene sets containing “WNT” in the name and see whether they tend to be differentially expressed in basal pregnant vs lactating.

Notes

• A common application of ROAST is to use a set of DE genes that was defined from an analysis of an independent data set. ROAST can then determine whether similar changes are observed in the contrast of interest for the current data set.
• ROAST estimates p-values by simulation, so the results may change slightly between runs. More precise p-values can be obtained by increasing the number of rotations, albeit at the cost of increased computational time.
• The smallest p-value that can be reported is 1/(2nrot + 1) where nrot is the number of rotations. This lower bound can be decreased by increasing nrot.

References

Supplementary Material - CAMERA

CAMERA gene set testing using the Broad’s curated gene sets

Other databases of gene sets that are available come from the Broad Institute’s Molecular Signatures Database (MSigDB). CAMERA is good option for testing a very large number of gene sets such as the MSigDB sets, as it is very fast. It has the advantage of accounting for inter-gene correlation within each gene set (Wu and Smyth 2012).

CAMERA takes as input the DGEList object dgeObj, the indexed list of gene sets c2.ind, the design matrix, the contrast being tested, as well as some other arguments. By default, CAMERA can estimate the correlation for each gene set separately. However, in practise, it works well to set a small inter-gene correlation of about 0.05 using the inter.gene.cor argument.

# Specify a design matrix without an intercept term
design <- model.matrix(~ CellType + Status, data = sampleinfo)
#design
# Check contrasts:
print(colnames(design))

[1] "(Intercept)"     "CellTypeluminal" "Statuspregnant"  "Statusvirgin"   

# Run analysis:
gst.camera <- camera.DGEList(dgeObj,index=c2.ind,design=design,contrast=2,inter.gene.cor=0.05)

CAMERA outputs a dataframe of the resulting statistics, with each row denoting a different gene set. The output is ordered by p-value so that the most significant should be at the top. Let’s look at the top 5 gene sets:

gst.camera[1:5,]

The total number of significant gene sets at 5% FDR is:

table(gst.camera$FDR < 0.05)

FALSE  TRUE 
 4699    23 

You can write out the camera results to a csv file to open in excel.

write.csv(gst.camera,file="gst_LumVsBas.csv")

Challenge 2

1. Run camera on the pregnant vs lactating contrast.
2. Run camera on a different set of MSigDB gene sets, the hallmark datasets, mouse_H_v5.rdata. You will need to load in the hallmark gene sets, and the object will be called Mm.H in R.

Lun, Aaron T L, Yunshun Chen, and Gordon K Smyth. 2016. “It’s DE-licious: A Recipe for Differential Expression Analyses of RNA-seq Experiments Using Quasi-Likelihood Methods in edgeR.” Methods in Molecular Biology (Clifton, N.J.) 1418 (January): 391–416. doi:10.1007/978-1-4939-3578-9\_19.

Wu, D, and G K Smyth. 2012. “Camera: a competitive gene set test accounting for inter-gene correlation.” Nucleic Acids Research 40 (17). Oxford University Press: e133—–e133.

Wu, D, E Lim, F Vaillant, M L Asselin-Labat, J E Visvader, and G K Smyth. 2010. “ROAST: rotation gene set tests for complex microarray experiments.” Bioinformatics 26 (17). Oxford Univ Press: 2176–82.