Multiomics workshop part 1: Working with transcriptome data
Quentin D. Read, Matt Burnham, & Matthias Anagho-Mattanovich
July 07, 2026
Background
We are going to analyze RNAseq transcriptome data. This dataset was collected because the researchers were interested in whether gene expression changes over time after the induction of lipolysis in mouse fat cells. The RNAseq technique counts the numbers of RNA transcripts from many different genes. The number of RNA transcripts of a gene that are found within a cell are an indication of how much that gene is being expressed. We can compare the level of expression of different genes, or of the same gene under different conditions. In this case, the conditions are three different time points (0h, 4h, 24h) after induction of lipolysis of a differentiated murine adipocyte cell line (BMC).
At each time point, six replicate samples were taken, for a total of 3 time points × 6 replicates = 18 samples. (By the way, it is unclear to me from reading the paper whether there were truly 18 independent samples or whether there were three repeated measures taken on the same six experimental units. For demo purposes today, we will ignore that issue and assume we have 18 independent samples.)
In this notebook, we are going to take the more or less raw counts and go through these steps. Along the way I’ll try to explain why exactly we’re doing what we’re doing, so don’t worry if you don’t already know the rationale behind each step of the pipeline.
- Explore the data
- Add metadata
- Remove duplicate genes
- Remove genes with low expression
- Transform counts to account for library size
- Further transform to meet assumptions of the linear models
- Fit differential expression model
- Apply different cutoffs to identify differentially expressed genes
- Create graphs and tables to visualize the results
You’ll notice this is not “multi” omics yet. We are only working with the transcriptomics data. But I think it’s good to start by working with individual datasets before we get into the fancy and complicated integration stuff.
Here are the original sources this lesson is based on:
Setup
Load packages
Load packages. If you are running this code on the remote server, the
packages should already be installed. If you are running this code on
your own machine, you may need to install some of the packages. Packages
that are hosted on CRAN, the main R package repository, can be installed
using install.packages('packagename'). Packages that are
hosted on Bioconductor can be installed by first installing the
BiocManager package from CRAN and then installing Bioconductor packages
using BiocManager::install('packagename').
# CRAN packages
library(ggplot2)
library(RColorBrewer)
library(gplots)
# Bioconductor packages
library(limma)
library(edgeR)
library(Mus.musculus)
library(ComplexHeatmap)Set plot theme for any ggplot2 plots we might make during this session. Define a vector of colors so that the plots we make will have consistent colors.
theme_set(
theme_bw(base_size = 14) +
theme(panel.grid = element_blank(),
strip.background = element_blank(),
axis.text = element_text(color = 'black'))
)
custom_colors <- c("salmon", "darkolivegreen3", "turquoise3", "mediumpurple1")Import data and add metadata
The data to be imported is in a file called
transcriptome.rda. RNAseq data consists of read counts for
a bunch of different genes. The data in this file has already been
converted from a raw CSV file containing the counts, into a format that
works with the R package edgeR. It is in the format of a list with three
elements, called counts, samples, and
offset. The counts object holds the count
data, the samples object holds metadata on each of the
samples, and the offset object includes normalization
factors that have already been calculated for each of the samples.
Load the data into your workspace.
load(file = 'data/transcriptome.rda')If you are running this code on your own laptop and not on the remote server, you will need to load the data from a URL on GitHub Pages where I am hosting it, using the code below.
load(file = url('https://usda-ree-ars.github.io/SEAStatsData/omics2026/transcriptome.rda'))Look at the dimensions of the counts object. There are
> 45,000 rows (genes) and 18 columns (samples).
dim(transcriptome$counts)## [1] 45705 18The samples metadata object does not currently have
information on which time point each sample comes from. This code
manually adds the group column.
groups <- as.factor(c("tp_0h", "tp_0h", "tp_4h", "tp_24h", "tp_0h", "tp_4h", "tp_24h", "tp_0h", "tp_4h", "tp_24h", "tp_4h", "tp_24h", "tp_0h", "tp_4h", "tp_24h", "tp_0h", "tp_4h", "tp_24h"))
transcriptome$samples$group <- groups
transcriptome$samples## group lib.size norm.factors
## 0172_1 tp_0h 26501223 1
## 0172_10 tp_0h 21138097 1
## 0172_11 tp_4h 15657210 1
## 0172_12 tp_24h 8924701 1
## 0172_13 tp_0h 21819125 1
## 0172_14 tp_4h 21806566 1
## 0172_15 tp_24h 23982143 1
## 0172_16 tp_0h 11258673 1
## 0172_17 tp_4h 7010742 1
## 0172_18 tp_24h 22232438 1
## 0172_2 tp_4h 19411202 1
## 0172_3 tp_24h 20539593 1
## 0172_4 tp_0h 21339104 1
## 0172_5 tp_4h 22875462 1
## 0172_6 tp_24h 22475172 1
## 0172_7 tp_0h 18850152 1
## 0172_8 tp_4h 11455607 1
## 0172_9 tp_24h 21314510 1Take a look at the first few rows of the counts object.
You can see that the rows are named with the Ensembl gene IDs. To get
more informative IDs, look up the Ensembl gene IDs using a lookup table
that is included in the Mus.musculus package.
head(transcriptome$counts)## 0172_1 0172_10 0172_11 0172_12 0172_13 0172_14 0172_15
## ENSMUSG00000000001 1618.000 1414 1217.000 820.000 1396.000 1737 1996
## ENSMUSG00000000003 0.000 0 0.000 0.000 0.000 0 0
## ENSMUSG00000000028 1351.999 1090 779.001 379.999 1086.999 1010 1034
## ENSMUSG00000000031 4.000 6 4.000 1.000 4.000 2 8
## ENSMUSG00000000037 51.001 32 24.000 14.000 38.000 19 40
## ENSMUSG00000000049 1.000 2 3.000 0.000 0.000 0 0
## 0172_16 0172_17 0172_18 0172_2 0172_3 0172_4 0172_5
## ENSMUSG00000000001 664 217 1912 1450.000 1717.000 1352 1716
## ENSMUSG00000000003 0 0 0 0.000 0.000 0 0
## ENSMUSG00000000028 548 112 939 1045.000 844.000 1088 1163
## ENSMUSG00000000031 1 0 1 4.000 2.000 5 12
## ENSMUSG00000000037 15 1 33 17.999 24.001 41 28
## ENSMUSG00000000049 0 0 2 0.000 1.000 2 2
## 0172_6 0172_7 0172_8 0172_9
## ENSMUSG00000000001 1991.000 1203 805 1762
## ENSMUSG00000000003 0.000 0 0 0
## ENSMUSG00000000028 929.000 957 556 942
## ENSMUSG00000000031 7.000 6 0 8
## ENSMUSG00000000037 24.001 33 9 36
## ENSMUSG00000000049 0.000 1 0 0geneid <- rownames(transcriptome)
genes <- AnnotationDbi::select(
Mus.musculus,
keys = geneid,
columns = c("SYMBOL", "TXCHROM", "ENTREZID"),
keytype = "ENSEMBL"
)
dim(genes)## [1] 46051 4There are a few more genes in the new lookup table than there are in
the original dataset, so there must be a few duplicates. Remove any
duplicates by keeping only the first instance of each of the Ensembl
IDs. Then add the genes lookup table to our data list.
genes <- genes[!duplicated(genes$ENSEMBL), ]
transcriptome$genes <- genesTransform and filter data
Log counts per million transformation
The 18 samples all have a different library size, as we saw above
when we looked at the samples object. If some of the
samples were sequenced to greater depth than others, they will end up
having higher counts. We cannot directly compare the raw counts across
different samples, without first adjusting for this difference in
size.
One kind of transformation we can do is the CPM (counts per million)
transformation. In these samples, library size varies from ~7 million at
the lowest sequencing depth to ~26 million at the highest. A gene with
CPM = 1 for this sample would have ~7 counts in the low-depth and ~26
counts in the high-depth sample. The function cpm() from
edgeR implements this transformation.
cpm <- cpm(transcriptome)The log-CPM transformation is the CPM on the log base 2 scale. Because many of the counts are zero, a small constant is added before taking the log, which is based on the average library size. This avoids taking the log of zero and also smooths out the variability for genes expressed at a low level.
The summary() method prints summary statistics for each
of the 18 samples. You can see that the distribution is very similar
across the samples. (Compare to summary(cpm) for the
untransformed counts per million.)
lcpm <- cpm(transcriptome, log = TRUE)
summary(lcpm)## 0172_1 0172_10 0172_11 0172_12
## Min. :-4.292 Min. :-4.426 Min. :-4.426 Min. :-4.2230
## 1st Qu.:-3.501 1st Qu.:-3.501 1st Qu.:-3.501 1st Qu.:-3.5012
## Median :-3.501 Median :-3.501 Median :-3.501 Median :-3.5012
## Mean :-1.233 Mean :-1.130 Mean :-0.984 Mean :-0.6904
## 3rd Qu.: 1.042 3rd Qu.: 1.301 3rd Qu.: 1.705 3rd Qu.: 2.6039
## Max. :14.753 Max. :15.169 Max. :15.425 Max. :16.1232
## 0172_13 0172_14 0172_15 0172_16
## Min. :-4.426 Min. :-4.426 Min. :-4.144 Min. :-4.4255
## 1st Qu.:-3.501 1st Qu.:-3.501 1st Qu.:-3.501 1st Qu.:-3.5012
## Median :-3.501 Median :-3.501 Median :-3.501 Median :-3.5012
## Mean :-1.138 Mean :-1.141 Mean :-1.191 Mean :-0.8101
## 3rd Qu.: 1.307 3rd Qu.: 1.303 3rd Qu.: 1.207 3rd Qu.: 2.1847
## Max. :15.204 Max. :15.141 Max. :14.640 Max. :16.6784
## 0172_17 0172_18 0172_2 0172_3
## Min. :-4.4255 Min. :-4.144 Min. :-4.155 Min. :-4.198
## 1st Qu.:-3.5012 1st Qu.:-3.501 1st Qu.:-3.501 1st Qu.:-3.501
## Median :-3.5012 Median :-3.501 Median :-3.501 Median :-3.501
## Mean :-0.6537 Mean :-1.140 Mean :-1.080 Mean :-1.109
## 3rd Qu.: 2.8859 3rd Qu.: 1.379 3rd Qu.: 1.440 3rd Qu.: 1.426
## Max. :21.1341 Max. :14.771 Max. :15.170 Max. :14.901
## 0172_4 0172_5 0172_6 0172_7
## Min. :-4.172 Min. :-4.086 Min. :-4.042 Min. :-4.426
## 1st Qu.:-3.501 1st Qu.:-3.501 1st Qu.:-3.501 1st Qu.:-3.501
## Median :-3.501 Median :-3.501 Median :-3.501 Median :-3.501
## Mean :-1.128 Mean :-1.175 Mean :-1.156 Mean :-1.070
## 3rd Qu.: 1.325 3rd Qu.: 1.162 3rd Qu.: 1.330 3rd Qu.: 1.456
## Max. :14.932 Max. :14.807 Max. :14.721 Max. :15.047
## 0172_8 0172_9
## Min. :-4.4255 Min. :-4.146
## 1st Qu.:-3.5012 1st Qu.:-3.501
## Median :-3.5012 Median :-3.501
## Mean :-0.8233 Mean :-1.135
## 3rd Qu.: 2.1527 3rd Qu.: 1.341
## Max. :16.3918 Max. :14.719Filtering by expression
An obvious step is to remove genes that have all zero counts across all the samples. They contribute no information at all. The below code calculates how many genes have zero counts (but doesn’t remove them yet). Almost 20,000 of the 45,000 genes have zero counts across all 18 samples.
table(rowSums(transcriptome$counts) == 0)##
## FALSE TRUE
## 26427 19278Removing the genes that are all zeros is a no-brainer. But genes that
are expressed at very low levels across all samples probably don’t have
any biological relevance either. We want to remove them as well. But
before we do, we have to define what we mean by “a very low level” of
expression. The function filterByExpr() from the edgeR
package keeps only genes that have a CPM value above a pre-specified
cutoff in a reasonable proportion of the samples. Using all default
arguments to this function, we filter out about 30,000 genes, leaving a
little less than 15,000.
Details on how the filtering is done: the default minimum read count is 10. Convert that to a CPM using the median library size of 21 million from the dataset; 10/21 \(\approx\) 0.48. Because each time point has six replicates, the filtering keeps all genes that have a CPM of at least 0.48 in six or more of the samples.
# Makes a copy of the unfiltered data, so that we can make the figure comparing filtered vs. unfiltered later.
transcriptome_unfiltered <- transcriptome
# Decide which genes to keep and then filter the transcriptome data by that criterion
keep.exprs <- filterByExpr(transcriptome, group = transcriptome$samples$group)
transcriptome <- transcriptome[keep.exprs, , keep.lib.sizes = FALSE]
dim(transcriptome)## [1] 14790 18Let’s take a look at the distribution of log CPM values in the original unfiltered data versus the filtered data. (Don’t worry about the code used to produce this figure.)
You can see that the filtering process has removed the large proportion of genes within each sampling with little or no expression.
L <- mean(transcriptome_unfiltered$samples$lib.size) * 1e-6
M <- median(transcriptome_unfiltered$samples$lib.size) * 1e-6
lcpm.cutoff <- log2(10/M + 2/L)
lcpm_unfiltered <- cpm(transcriptome_unfiltered, log = TRUE)
lcpm_filtered <- cpm(transcriptome, log = TRUE)
nsamples <- ncol(transcriptome_unfiltered)
col <- brewer.pal(nsamples, "Set3")
par(mfrow = c(1, 2))
plot(density(lcpm_unfiltered[, 1]), col = col[1], lwd = 2, ylim = c(0,0.7), las = 2, main="", xlab="")
title(main = "A. Raw data", xlab = "Log-cpm")
abline(v = lcpm.cutoff, lty = 3)
for (i in 2:nsamples) {
den <- density(lcpm_unfiltered[, i])
lines(den$x, den$y, col = col[i], lwd = 2)
}
plot(density(lcpm_filtered[, 1]), col = col[1], lwd=2, ylim = c(0,0.26), las = 2, main="", xlab="")
title(main = "B. Filtered data", xlab = "Log-cpm")
abline(v = lcpm.cutoff, lty = 3)
for (i in 2:nsamples) {
den <- density(lcpm_filtered[, i])
lines(den$x, den$y, col = col[i], lwd = 2)
}
Normalize the remaining data
Typically the next step would be to normalize the remaining data. Normalizing means getting rid of any unwanted external factors that may influence expression of individual samples. We want to be able to make a fair comparison of the expression within each gene between the different groups. To do that, all samples need to have a similar distribution of expression values.
The figure we just made already shows us that the distributions of
log-CPM values are extremely similar between the different samples, so
in this case normalization will not have a very big effect. But it’s a
good idea to do it anyway. calcNormFactors() from edgeR
adds the normalization factors as a new column in the
samples object of our transcriptome data list. You can see
below that most of the calculated normalization factors are all pretty
close to 1 (between ~0.90 and ~1.15), indicating that normalization is
only having a modest effect. One sample has a normalization factor of
0.43, so normalization is going to have a bigger impact on that one.
transcriptome <- calcNormFactors(transcriptome, method = "TMM")## Warning in calcNormFactors.DGEList(transcriptome, method = "TMM"): object contains offsets, which take precedence over library
## sizes and norm factors (and which will not be recomputed).transcriptome$samples$norm.factors## [1] 1.0075305 0.9961379 1.0483793 1.1524165 1.0190427 1.0328989 1.1189970
## [8] 0.9128948 0.4306014 1.1252224 1.0512063 1.1015828 1.0328745 1.0489316
## [15] 1.1211031 1.0366249 0.9901253 1.0965377As it turns out, the transcriptome data we started with already has offset values calculated. These are more highly resolved normalization factors that are applied to every single observation separately, that were calculated by a previous edgeR step. Offsets take precedence over the normalization factors we just calculated.
Investigate the structure of the data
Multidimensional scaling plot
Ultimately, we are interested in looking at expression patterns across all genes to see whether there are time trends that are common across genes. Before we do a formal analysis, we can use “quick and dirty” dimension reduction methods like multidimensional scaling (MDS) or principal components analysis (PCA) to visualize our high-dimensional data in fewer dimensions. We have about 15,000 genes but we can make a multidimensional scaling plot to show similarities and dissimilarities between our samples.
The basic way MDS, and more or less all dimensionality reduction methods, work, is that it tries to reduce N-dimensional data into a smaller number of dimensions, while preserving as much information from the original data as possible. Here, our N is ~15,000. MDS takes the distances between each pair from our 18 samples in N-dimensional space and puts them into a lower-dimensional space in such a way that the distances between any two samples are preserved as well as possible. If the lower number of dimensions is 2, we can display it on a scatterplot. I’ll skip any more details on that.
par(mfrow = c(1, 1))
plotMDS(
lcpm_filtered,
labels = groups,
col = as.character(factor(groups, labels = custom_colors[1:3]))
)
The MDS plot shows us that the samples from each time point cluster very well in 2-dimensional space, with one major exception. One of the samples from the 4-hour time point is very dissimilar to all the other ones (sample 17). If you look back at the library sizes for the different samples, it is the one with by far the lowest sequencing depth. So, it seems like something went wrong with that sample. You might consider excluding it from the analysis.
Relationship between mean and variance
We are going to use a linear model approach to look at which genes are expressed at different levels in the different time points. An important assumption of a linear model that is based on normally distributed error is that there is no relationship between the mean and the variance. If that assumption is violated, we might end up drawing the wrong conclusions. Count data, like we have here, often have skewed distributions with a lot of small values at or near zero and a few high values, because it is not possible to have a count below zero. We usually deal with this issue by either transforming the data or fitting special statistical models that do not have to assume the variation around each mean is normally distributed.
In this case, we are going to use the function voom from
the limma package to deal with the mean-variance relationship issue.
(Apparently “voom” stands for “mean-variance modeling at the
observational level,” though I am not sure how they got those letters
from that.) The voom algorithm takes the log-counts per million values
(after we’ve already filtered the ones with very low counts out), fits
linear models to each gene separately to estimate the relationship
between the mean and the variance, then weights each observation to
adjust it based on its predicted variance. Observations with higher
predicted variance get a lower weight. This flattens out the
mean-variance trend.
The first step is to set up a design matrix, which is a matrix of zeros and ones that indicates which group (time point in this case) each sample belongs to. The design matrix object also includes a contrasts attribute, which indicates which of the group means have to be subtracted from each other to get the differences between each pair of means.
design <- model.matrix(~ 0 + groups)
colnames(design) <- gsub("groups", "", colnames(design))
design## tp_0h tp_24h tp_4h
## 1 1 0 0
## 2 1 0 0
## 3 0 0 1
## 4 0 1 0
## 5 1 0 0
## 6 0 0 1
## 7 0 1 0
## 8 1 0 0
## 9 0 0 1
## 10 0 1 0
## 11 0 0 1
## 12 0 1 0
## 13 1 0 0
## 14 0 0 1
## 15 0 1 0
## 16 1 0 0
## 17 0 0 1
## 18 0 1 0
## attr(,"assign")
## [1] 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$groups
## [1] "contr.treatment"contr.matrix <- makeContrasts(
Comp_0h_v_4h = tp_0h - tp_4h,
Comp_0h_v_24h = tp_0h - tp_24h,
Comp_4h_v_24h = tp_4h - tp_24h,
levels = colnames(design))
contr.matrix## Contrasts
## Levels Comp_0h_v_4h Comp_0h_v_24h Comp_4h_v_24h
## tp_0h 1 1 0
## tp_24h 0 -1 -1
## tp_4h -1 0 1Now we can just pass the filtered transcriptome data object and the
design matrix to the voom() function to estimate the
mean-variance relationship.
v <- voom(transcriptome, design, plot=TRUE)
This is a plot with mean on the x-axis and variance on the y-axis. We
see above that as the mean increases, the variance tends to decrease.
This is what we usually expect from log counts per million data.
voom() has just calculated weights for each observation
that will adjust for that.
Here are the first few rows of the weights calculated by
voom(). The 9th column stands out — again, this is the
“weird” sample we noticed earlier in the MDS plot.
head(v$weights)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1.753934 1.743821 1.744413 1.6936536 1.746672 1.756615 1.755830 1.6859297
## [2,] 1.749544 1.737217 1.729567 1.6327318 1.740630 1.748031 1.740169 1.6701506
## [3,] 1.465884 1.360404 1.141549 0.9658885 1.333331 1.330086 1.410401 0.9401023
## [4,] 1.753231 1.743977 1.737733 1.5763178 1.746746 1.754349 1.732323 1.6844632
## [5,] 1.767950 1.763623 1.758714 1.7143381 1.764127 1.766133 1.762052 1.7338848
## [6,] 1.768709 1.763859 1.749672 1.6919756 1.765009 1.759840 1.755355 1.7374208
## [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
## [1,] 1.4047076 1.753691 1.753177 1.749610 1.746292 1.758345 1.753853 1.739137
## [2,] 1.2968968 1.736233 1.743439 1.729892 1.739511 1.750333 1.736598 1.731889
## [3,] 0.6789118 1.226548 1.262895 1.344248 1.383777 1.357856 1.387398 1.343388
## [4,] 1.3230237 1.727952 1.749354 1.720155 1.747296 1.755812 1.726560 1.740019
## [5,] 1.5536456 1.760967 1.764209 1.757240 1.763856 1.767398 1.760305 1.761555
## [6,] 1.4613976 1.753104 1.757002 1.749001 1.764824 1.761360 1.753279 1.762276
## [,17] [,18]
## [1,] 1.7204320 1.750906
## [2,] 1.7008335 1.731796
## [3,] 0.9438952 1.168115
## [4,] 1.7145544 1.715150
## [5,] 1.7446048 1.759047
## [6,] 1.7283401 1.750294Differential expression analysis
Fit linear models and calculate test statistics
Now we can do differential expression analysis using the weights that
voom() calculated to adjust for the mean-variance
relationship in the data.
The lmFit() function fits the linear models for each
gene, using the inverse variance weights we just calculated using
voom(). Then, we use contrasts.fit() to take
the pairwise contrasts between each pair of time points: 0h versus 4h,
0h versus 24h, and 4h versus 24h. Next, eBayes() tests
whether each contrast is equal to zero, using a method called empirical
Bayes to bring all the variances for each gene closer to a common value.
plotSA() shows us the mean-variance trend after the
correction. Looks much better than before!
vfit <- lmFit(v, design)
vfit <- contrasts.fit(vfit, contrasts = contr.matrix)
efit <- eBayes(vfit)
plotSA(efit)
Do tests for differential expression
The function decideTests() identifies which genes are
significantly differentially expressed between each pair of time points.
By default, a false discovery rate correction is used to adjust for the
multiple tests we’re doing across many genes — if we didn’t make that
adjustment, we would get a lot of false positives because if you test
thousands of genes, some will be significant by random chance alone. (In
fact, if we use uncorrected \(p <
0.05\) as the significance threshold, 5% of the genes, or 1 in
20, would show a significant difference even if the data were just
noise. Out of almost 15,000 genes, that’s about 750. So, it’s important
to correct for multiple testing.)
summary(decideTests(efit))## Comp_0h_v_4h Comp_0h_v_24h Comp_4h_v_24h
## Down 434 146 37
## NotSig 14273 14320 14210
## Up 83 324 543The summary() method shows how many genes were
significantly “down” (less expressed) or “up” (more expressed) between
each pair of time points. For example 0h versus 4h has 434 significantly
“down” genes, meaning significantly less expressed at 0h and more
expressed at 4h, based on the way we set up the contrast matrix. 83
genes have significantly higher expression at 0h than 4h, and the
remainder are not significantly different.
Adjust tests using a cutoff
We may also be interested in a specific cutoff that we consider to be
a biologically relevant change in expression. Statistically significant
doesn’t always mean practically significant. We can supply a log fold
change cutoff to decideTests(). Now, genes are only
identified as being differentially expressed if they are both
statistically significant and the difference in expression is
above a certain cutoff.
Here we’ll use the treat() function to update the
specify lfc = 1 cutoff instead of using the empirical Bayes
adjustment. A log-2 fold change of 1 means one expression has to be at
least 2 times as high as the other to be considered a significant
difference.
tfit <- treat(vfit, lfc = 1)
dt <- decideTests(tfit)
summary(dt)## Comp_0h_v_4h Comp_0h_v_24h Comp_4h_v_24h
## Down 27 13 2
## NotSig 14750 14699 14760
## Up 13 78 28We end up with a much smaller set of genes to look at.
Visualizing differential expression
There are a few different ways to visualize differential expression.
Venn diagrams
Venn diagrams are a fairly simple, some might say oversimplified, way of visualizing differential expression.
Because we have three time points, we want to know if some genes were significantly more (or less) expressed at both 4 hours and 24 hours compared to the baseline (0h) time point. These genes are likely to be important because their change in expression relative to the baseline was sustained over a long time period.
The dt object we generated from
decideTests() with log-fold-change cutoff of 1 is a matrix
with three columns and a row for each gene. For each of the comparisons
between time points, it has a 0 if the difference was not significant, a
-1 if the difference was significantly negative, and a 1 if the
difference was significantly positive.
Here we are identifying which rows have a significant difference,
whether positive or negative, between 0 and 4 and also between 0 and 24.
There are 14 genes that fit that description. Here are their names. One
gene had an Ensembl ID but no annotation, hence why it shows up as
NA.
de.common <- which(dt[, 1] != 0 & dt[, 2] != 0)
length(de.common)## [1] 14tfit$genes$SYMBOL[de.common]## [1] "Fosb" "Btg2" "Fos" "Nr4a1" "Il1rl1" "Adcyap1r1"
## [7] "Egr3" "Egr2" "Ier2" "Alkal2" "Unc93a2" "Unc93a"
## [13] "Ttll2" NAThis Venn diagram indicates how many genes are differentially expressed only between 0 and 4 hours, only between 0 and 24 hours, and between both. The 14673 floating in the corner is the number of genes not significantly differentially expressed between either of those two pairs.
vennDiagram(dt[, 1:2], circle.col = custom_colors[c(3, 1)])
Here’s a similar Venn diagram looking at how many genes are different from the middle (4h) time point; we pull out column 1 for 0h versus 4h and column 3 for 4h versus 24h.
There are only two genes, including leptin, that are differentially expressed both between 0 and 4 and between 4 and 24 hours.
de.common <- which(dt[, 1] != 0 & dt[, 3] != 0)
tfit$genes$SYMBOL[de.common]## [1] "Nhsl1" "Lep"vennDiagram(dt[, c(1, 3)], circle.col = custom_colors[c(3, 1)])
No genes are significantly differentially expressed between all three pairs of time points.
de.common <- which(dt[, 1] != 0 & dt[, 2] != 0 & dt[, 3] != 0)
length(de.common)## [1] 0vennDiagram(dt, circle.col = custom_colors[c(3, 1, 2)])
UpSet plot
An UpSet plot is another way of showing how many genes are differentially expressed between each pair of time points, and how much overlap there is across all the pairs. It is a more information-rich alternative to Venn diagrams. (Venn diagrams with more than three circles are really hard to read. I’d recommend an UpSet plot if you have more than three groups.)
nonzero_abs_dt <- abs(dt[apply(dt, 1, function(row) !all(row ==0)), ])
grays <- c("gray30", "gray15", "gray0")
nonzero_comb_matrix <- make_comb_mat(nonzero_abs_dt)
UpSet(nonzero_comb_matrix, pt_size = unit(5, "mm"),
bg_col = "#F0F0F0", lwd = 3,
column_title = "Overlapping differential genes",
comb_col = grays[comb_degree(nonzero_comb_matrix)],
right_annotation = upset_right_annotation(nonzero_comb_matrix,
gp = gpar(col = custom_colors, fill = custom_colors),
annotation_name_side = "top",
axis_param = list(side = "top")))
Looking at individual genes
The topTreat() function shows us the genes with the most
different expression between any given pair of time points. By default
it shows us the ones with the lowest adjusted p-values, but we can
supply the argument sort.by = 'logFC' to show the ones with
the largest log-fold changes. Only a few rows of one of the comparisons
is shown here but you can explore the others. The list of genes is
similar whether p-value or log-fold change is used.
base.vs.mid <- topTreat(tfit, coef = 1, n = 20)
base.vs.mid_bylogFC <- topTreat(tfit, coef = 1, n = 20, sort.by = 'logFC')
base.vs.post <- topTreat(tfit, coef = 2, n = 20)
mid.vs.post <- topTreat(tfit, coef = 3, n = 20)
head(base.vs.mid)## ENSEMBL ENTREZID SYMBOL TXCHROM logFC
## ENSMUSG00000003545 ENSMUSG00000003545 14282 Fosb chr7 6.168529
## ENSMUSG00000021250 ENSMUSG00000021250 14281 Fos chr12 4.655665
## ENSMUSG00000026069 ENSMUSG00000026069 17082 Il1rl1 chr1 -4.795021
## ENSMUSG00000037868 ENSMUSG00000037868 13654 Egr2 chr10 3.813263
## ENSMUSG00000059201 ENSMUSG00000059201 16846 Lep chr6 3.892108
## ENSMUSG00000028001 ENSMUSG00000028001 14161 Fga chr3 -4.427314
## AveExpr t P.Value adj.P.Val
## ENSMUSG00000003545 2.469887 11.707407 5.966604e-32 8.824608e-28
## ENSMUSG00000021250 3.892827 8.657653 2.423589e-18 1.792244e-14
## ENSMUSG00000026069 1.944959 -8.048607 4.207542e-16 2.074318e-12
## ENSMUSG00000037868 3.913820 6.715447 9.396963e-12 3.474527e-08
## ENSMUSG00000059201 2.223239 6.271900 1.786607e-10 4.566935e-07
## ENSMUSG00000028001 -1.072130 -6.266241 1.852712e-10 4.566935e-07We can also use plotMD() to make a mean-difference plot,
which is basically a sideways volcano plot, for all three pairs of time
points. The x-axis on each plot shows the average of the log of
expression between the two time points, and the y-axis is the log-fold
change. Pink points above the zero line are those that are significantly
more expressed at the earlier of the two time points, and blue points
below the midline are those that are significantly more expressed at the
later time point.
plotMD(tfit, column = 1, status = dt[ ,1], main = colnames(tfit)[1], xlim = c(-5,17),
hl.col = custom_colors)
plotMD(tfit, column = 2, status = dt[ ,2], main = colnames(tfit)[2], xlim = c(-5,17),
hl.col = custom_colors)
plotMD(tfit, column = 3, status = dt[ ,3], main = colnames(tfit)[3], xlim = c(-5,17),
hl.col = custom_colors)
Heat map
Here is a heat map of the 20 genes with the most significant
differences across the three pairs of time points. The argument
scale = 'row' means that log-expression values within each
gene are standardized to a z-score. Then, a simple clustering algorithm
is applied to group similar samples together.
In this case, samples from the same time point tend to cluster together. We also see that most of the genes shown here (which are the genes with the top 20 most significant p-values) are highly expressed at time point 0 and then not as highly expressed later on. However there are also a few genes that are not expressed very much at time point 0 but increase in expression at the later time points.
base.vs.mid <- topTreat(tfit, coef = 1, n = Inf)
base.vs.post <- topTreat(tfit, coef = 2, n = Inf)
mid.vs.post <- topTreat(tfit, coef = 3, n = Inf)
topgenes_df <- rbind(base.vs.mid, base.vs.post, mid.vs.post)
topgenes_ordered_df <- topgenes_df[order(topgenes_df$adj.P.Val), ]
topgenes <- unique(c(topgenes_ordered_df$ENSEMBL))[1:20]
topgene_index <- which(v$genes$ENSEMBL %in% topgenes)
# Create a continuous color palette using colorRampPalette
heatmap_colors <- colorRampPalette(c("turquoise4", "turquoise3", "white", "salmon", "darkred"))(100)
heatmap.2(lcpm_filtered[topgene_index,],
scale = "row",
labRow = v$genes$SYMBOL[topgene_index],
labCol = factor(groups, labels = c("0h", "4h", "24h")),
col = heatmap_colors, trace = "none", density.info = "none",
dendrogram = "both")
Individual gene expression
It is also straightforward to extract log-expression values for an individual gene and plot them. For example let’s look at leptin expression. Medians are plotted as diamonds over the data.
leptin_index <- which(genes$SYMBOL == 'Lep')
leptin_data <- data.frame(time_point = factor(groups, labels = c("0h", "4h", "24h")), LCPM = lcpm[leptin_index, ])
ggplot(leptin_data, aes(x = time_point, color = time_point, y = LCPM)) +
geom_point(size = 3, alpha = 0.5) +
stat_summary(fun = 'median', geom = 'point', shape = 5, size = 5) +
theme(legend.position = 'none') +
scale_color_manual(values = custom_colors) +
labs(x = 'Time point') +
ggtitle('Leptin expression')