Differential Expression and Multiple Comparisons

Biological versus technical variability

The t-test takes into account within population variability. In experimental biology, this variability can be divided into two main classes: technical and biological.

In general, the variability we observe across biological units, such as individuals, within a population is referred to as biological. We refer to the variability we observe across measurements of the same biological unit, such a aliquots from the same biological sample, as technical. Because newly developed measurement technologies are common in genomics, technical replicates are used many times to assess experimental data. By generating measurements from samples that are designed to be the same, we are able to measure and assess technical variability. We also use the terminology biological replicates and technical replicates to refer to samples from which we can measure biological and technical variability respectively.

It is important not to confuse biological and technical variability when performing statistical inference as the interpretation is quite different. For example, when analyzing data from technical replicates, the population is just the one sample from which these come from as opposed to more general population such as healthy humans or control mice.

Simple T-Tests Code:

library(genefilter)
rowttests(assay_data, grouping_factor)

Gene Set Analysis: Inference on Coordinated Expression Changes

Link

Here, we will explore software for testing differential expression in a set of genes. These tests differ from the gene-by-gene tests we saw previously. Again, the gene set testing software we will use lives in the limma package.

We download an experiment from the GEO website, using the getGEO function from the GEOquery package:

library(GEOquery)
g <- getGEO("GSE34313")
e <- g[1](/abishpius/R-for-Computational-Biology/wiki/1)

This dataset is hosted by GEO at the following link: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE34313

The experiment is described in the paper by Masuno 2011.

Briefly, the investigators applied a glucocorticoid hormone to cultured human airway smooth muscle. The glucocorticoid hormone is used to treat asthma, as it reduces the inflammation response, however it has many other effects throughout the different tissues of the body.

The groups are defined in the characteristics_ch1.2 variable:

e$condition <- e$characteristics_ch1.2
levels(e$condition) <- c("dex24","dex4","control")
table(e$condition)

By examining boxplots, we can guess that the data has already been normalized somehow, and on the GEO site the investigators report that they normalized using Agilent software.

We will subset to the control samples and the samples treated with dexamethasone (the hormone) after 4 hours.

boxplot(exprs(e), range=0)
names(fData(e))
lvls <- c("control", "dex4")
es <- e[,e$condition %in% lvls]
es$condition <- factor(es$condition, levels=lvls)

The following lines run the linear model in limma. We note that the top genes are common immune-response genes (CSF2, LIF, CCL2, IL6). Also present is FKBP5, a gene which regulates and is regulated by the protein which receives the glucocorticoid hormone.

library(limma)
design <- model.matrix(~ es$condition)
fit <- lmFit(es, design=design)
fit <- eBayes(fit)
tt <- topTable(fit, coef=2, genelist=fData(es)$GENE_SYMBOL)
tt

We will use the ROAST method for gene set testing. We can test a single gene set by looking up the genes which contain a certain GO ID, and providing this to the roast function. We will show how to get such lists of genes associated with a GO ID in the next chunk.

The roast function performs an advanced statistical technique, rotation of residuals, in order to generate a sense of the null distribution for the test statistic. The test statistic in this case is the summary of the scores from each gene. The tests are self-contained because only the summary for a single set is used, whereas other gene set tests might compare a set to all the other genes in the dataset, e.g., a competitive gene set test.

The result here tells us that the immune response genes are significantly down-regulated, and additionally, mixed up and down.

# Immune response
idx <- grep("GO:0006955", fData(es)$GO_ID)
length(idx)
r1 <- roast(es, idx, design)
# ?roast
r1

Testing multiple gene sets

We can also use the mroast function to perform multiple roast tests. First we need to create a list, which contains the indices of genes in the ExpressionSet for each of a number of gene sets. We will use the org.Hs.eg.db package to gather the gene set information.

# BiocManager::install("org.Hs.eg.db")
library(org.Hs.eg.db)
org.Hs.egGO2EG
go2eg <- as.list(org.Hs.egGO2EG)
head(go2eg)

The following code unlists the list, then gets matches for each Entrez gene ID to the index in the ExpressionSet. Finally, we rebuild the list.

govector <- unlist(go2eg)
golengths <- sapply(go2eg, length)
head(fData(es)$GENE)
idxvector <- match(govector, fData(es)$GENE)
table(is.na(idxvector))
idx <- split(idxvector, rep(names(go2eg), golengths))
go2eg[1](/abishpius/R-for-Computational-Biology/wiki/1)
fData(es)$GENE[idx[1](/abishpius/R-for-Computational-Biology/wiki/1)]

We need to clean this list such that there are no NA values. We also clean it to remove gene sets which have less than 10 genes.

idxclean <- lapply(idx, function(x) x[!is.na(x)])
idxlengths <- sapply(idxclean, length)
idxsub <- idxclean[idxlengths > 10]
length(idxsub)

The following line of code runs the multiple ROAST test. This can take about 3 minutes.

r2 <- mroast(es, idxsub, design)
head(r2)
r2 <- r2[order(r2$PValue.Mixed),]

We can use the GO.db annotation package to extract the GO terms for the top results, by the mixed test.

# BiocManager::install("GO.db")
library(GO.db)
columns(GO.db)
keytypes(GO.db)
GOTERM[rownames(r2)[1](/abishpius/R-for-Computational-Biology/wiki/rownames(r2)[1)]
r2tab <- select(GO.db, keys=rownames(r2)[1:10],
                columns=c("GOID","TERM","DEFINITION"), 
                keytype="GOID")
r2tab[,1:2]

We can also look for the top results using the standard p-value and in the up direction.

r2 <- r2[order(r2$PValue),]
r2tab <- select(GO.db, keys=rownames(r2)[r2$Direction == "Up"][1:10],
                columns=c("GOID","TERM","DEFINITION"), 
                keytype="GOID")
r2tab[,1:2]

Again but for the down direction.

r2tab <- select(GO.db, keys=rownames(r2)[r2$Direction == "Down"][1:5],
                columns=c("GOID","TERM","DEFINITION"), 
                keytype="GOID")
r2tab[,1:2]