What is a gene set?

Any a priori classification of `genes’ into biologically relevant groups

• Members of same biochemical pathway
• Proteins expressed in identical cellular compartments
• Co-expressed under certain conditions
• Targets of the same regulatory elements
• On the same cytogenic band
• ...

Sets do not need to be...

• Exhaustive
• Disjoint

Collections of gene sets

Gene Ontology (GO) Annotation (GOA)

• CC Cellular Components
• BP Biological Processes
• MF Molecular Function

Pathways

• MSigDb
• KEGG (no longer freely available)
• reactome
• PantherDB
• ...

E.g., MSigDb

• c1 Positional gene sets – chromosome & cytogenic band
• c2 Curated Gene Sets from online pathway databases, publications in PubMed, and knowledge of domain experts.
• c3 motif gene sets based on conserved cis-regulatory motifs from a comparative analysis of the human, mouse, rat, and dog genomes.
• c4 computational gene sets defined by mining large collections of cancer-oriented microarray data.
• c5 GO gene sets consist of genes annotated by the same GO terms.
• c6 oncogenic signatures defined directly from microarray gene expression data from cancer gene perturbations.
• c7 immunologic signatures defined directly from microarray gene expression data from immunologic studies.

Approach 1: hypergeometric tests

Steps

1. Classify each gene as ‘differentially expressed’ DE or not, e.g., based on P < 0.05
2. Are DE genes in the set more common than DE genes not in the set?

3. Fisher hypergeometric test, via fiser.test() or GOstats

Notes

• Conditional hypergeometric to accommodate GO DAG, GOstats
• But: artificial division into two groups

Approach 2: enrichment score

• Mootha et al., 2003; modified Subramanian et al., 2005.

Steps

• Sort genes by log fold change
• Calculate running sum: incremented when gene in set, decremented when not.
• Maximum of the running sum is enrichment score ES; large ES means that genes in set are toward top of list.
• Permuting subject labels for signficance

Approach 3: category \(t\)-test

E.g., Jiang & Gentleman, 2007;

• Summarize \(t\) (or other) statistic in each set
• Test for significance by permuting the subject labels
• Much more straight-forward to implement

Expression in NEG vs BCR/ABL samples for genes in the ‘ribosome’ KEGG pathway; Category vignette.

Competitive versus self-contained null hypothesis

Goemann & Bühlmann, 2007

• Competitive null: The genes in the gene set do not have stronger association with the subject condition than other genes. (Approach 1, 2)
• Self-contained null: The genes in the gene set do not have any association with the subject condition. (Approach 3)
• Probably, self-contained null is closer to actual question of interest
• Permuting subjects (rather than genes) is appropriate

Approach 4: linear models

E.g., Hummel et al., 2008,

• Colorectal tumors have good (‘stage II’) or bad (‘stage III’) prognosis. Do genes in the p53 pathway (just one gene set!) show different activity at the two stages?
• Linear model incorporates covariates – sex of patient, location of tumor

limma

• Majewski et al., 2010 romer() and Wu & Smythe 2012 camera() for enrichment (competitive null) linear models
• Wu et al., 2010: roast(), mroast() for self-contained null linear models

Approach 5: pathway topology

E.g., Tarca et al., 2009,

• Incorporate pathway topology (e.g., interactions between gene products) into signficance testing

• Signaling Pathway Impact Analysis

• Combined evidence: pathway over-representation \(P_{NDE}\); unusual signaling \(P_{PERT}\) (equation 1 of Tarca et al.)

Evidence plot, colorectal cancer. Points: pathway gene sets. Significant after Bonferroni (red) or FDR (blue) correction.

Issues with sequence data?

• All else being equal, long genes receive more reads than short genes
• Per-gene \(P\) values proportional to gene size

E.g., Young et al., 2010, goseq

• Hypergeometric, weighted by gene size
• Substantial differences
• Better: read depth??

DE genes vs. transcript length. Points: bins of 300 genes. Line: fitted probability weighting function.

Approach 6: de novo discovery

• So far: analogous to supervised machine learning, where pathways are known in advance
• What about unsupervised discovery?

Example: Langfelder & Hovarth, WGCNA

• Weighted correlation network analysis
• Described in Langfelder & Horvath, 2008

Representing gene sets in R

• Named list(), where names of the list are sets, and each element of the list is a vector of genes in the set.
• data.frame() of set name / gene name pairs
• GSEABase

Conclusions

Gene set enrichment classifications

• Kharti et al: Over-representation analysis; functional class scoring; pathway topology
• Goemann & Bühlmann: Competitive vs. self-contained null

Selected Packages

1-6 Experimental design, ..., Analysis of gene differential expression

This (relatively old) experiment examined the effects of androgen stimulation on a human prostate cancer cell line, LNCaP (Li et al., 2008). The experiment used short (35bp) single-end reads from 4 control and 3 untreated lines. Reads were aligned to hg19 using Bowtie, and counted using ENSEMBL 54 gene models.

Input the data to edgeR ’s DGEList data structure.

library(edgeR)
path <- system.file(package="goseq", "extdata", "Li_sum.txt")
table.summary <- read.table(path, sep='\t', header=TRUE, stringsAsFactors=FALSE)
counts <- table.summary[,-1]
rownames(counts) <- table.summary[,1]
grp <- factor(rep(c("Control","Treated"), times=c(4,3)))
summarized <- DGEList(counts, lib.size=colSums(counts), group=grp)

Use a ‘common’ dispersion estimate, and compare the two groups using an exact test

disp <- estimateCommonDisp(summarized)
tested <- exactTest(disp)
topTags(tested)

## Comparison of groups: Treated-Control 
## logFC logCPM PValue FDR
## ENSG00000127954 11.557868 6.680748 2.574972e-80 1.274766e-75
## ENSG00000151503 5.398963 8.499530 1.781732e-65 4.410322e-61
## ENSG00000096060 4.897600 9.446705 7.983756e-60 1.317479e-55
## ENSG00000091879 5.737627 6.282646 1.207655e-54 1.494654e-50
## ENSG00000132437 -5.880436 7.951910 2.950042e-52 2.920896e-48
## ENSG00000166451 4.564246 8.458467 7.126763e-52 5.880292e-48
## ENSG00000131016 5.254737 6.607957 1.066807e-51 7.544766e-48
## ENSG00000163492 7.085400 5.128514 2.716461e-45 1.681014e-41
## ENSG00000113594 4.051053 8.603264 9.272066e-44 5.100255e-40
## ENSG00000116285 4.108522 7.864773 6.422468e-43 3.179507e-39

7. Comprehension

Start by extracting all P values, then correcting for multiple comparison using p.adjust(). Classify the genes as differentially expressed or not.

padj <- with(tested$table, {
 keep <- logFC != 0
 value <- p.adjust(PValue[keep], method="BH")
 setNames(value, rownames(tested)[keep])
})
genes <- padj < 0.05
table(genes)

## genes
## FALSE TRUE 
## 19535 3208

pwf <- nullp(genes,"hg19","ensGene")

## Loading hg19 length data...

## Warning in pcls(G): initial point very close to some inequality
## constraints

head(pwf)

## DEgenes bias.data pwf
## ENSG00000230758 FALSE 247 0.03757470
## ENSG00000182463 FALSE 3133 0.20436865
## ENSG00000124208 FALSE 1978 0.16881769
## ENSG00000230753 FALSE 466 0.06927243
## ENSG00000224628 FALSE 1510 0.15903532
## ENSG00000125835 FALSE 954 0.12711992

GO.wall <- goseq(pwf, "hg19", "ensGene")

## Fetching GO annotations...
## For 9751 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...

head(GO.wall)

## category over_represented_pvalue under_represented_pvalue
## 10729 GO:0044763 8.237627e-15 1
## 10708 GO:0044699 2.079753e-14 1
## 2453 GO:0005737 2.956026e-10 1
## 3004 GO:0006614 6.131543e-09 1
## 7499 GO:0031982 1.101818e-08 1
## 2372 GO:0005576 1.339207e-08 1
## numDEInCat numInCat
## 10729 1893 8355
## 10708 2054 9201
## 2453 1790 8080
## 3004 39 107
## 7499 638 2675
## 2372 669 2836
## term ontology
## 10729 single-organism cellular process BP
## 10708 single-organism process BP
## 2453 cytoplasm CC
## 3004 SRP-dependent cotranslational protein targeting to membrane BP
## 7499 vesicle CC
## 2372 extracellular region CC

GO.nobias <- goseq(pwf,"hg19","ensGene",method="Hypergeometric")

## Fetching GO annotations...
## For 9751 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...

idx <- match(GO.nobias$category, GO.wall$category)
plot(log10(GO.nobias[, "over_represented_pvalue"]) ~
 log10(GO.wall[idx, "over_represented_pvalue"]),
 xlab="Wallenius", ylab="Hypergeometric",
 xlim=c(-5, 0), ylim=c(-5, 0))
abline(0, 1, col="red", lwd=2)