From reads to genes to pathways: differential expression analysis of RNA-Seq experiments using Rsubread and the edgeR quasi-likelihood pipeline
R version: R version 4.6.1 (2026-06-24)
Bioconductor version: 3.23
Package: 1.37.0
Abstract
RNA sequencing (RNA-seq) has become a very widely used technology for profiling gene expression. One of the most common aims of RNA-seq profiling is to identify genes or molecular pathways that are differentially expressed (DE) between two or more biological conditions. This article demonstrates a computational workflow for the detection of DE genes and pathways from RNA-seq data by providing a complete analysis of an RNA-seq experiment profiling epithelial cell subsets in the mouse mammary gland. The workflow uses R software packages from the open-source Bioconductor project and covers all steps of the analysis pipeline, including alignment of read sequences, data exploration, differential expression analysis, visualization and pathway analysis. Read alignment and count quantification is conducted using the Rsubread package and the statistical analyses are performed using the edgeR package. The differential expression analysis uses the quasi-likelihood functionality of edgeR.
Introduction
RNA sequencing (RNA-seq) has become a very widely used technology for profiling transcriptional activity in biological systems. One of the most common aims of RNA-seq profiling is to identify genes or molecular pathways that are differentially expressed (DE) between two or more biological conditions. Changes in expression can then be associated with differences in biology, providing avenues for further investigation into potential mechanisms of action.
This article provides a detailed workflow for analyzing an RNA-seq study all the way from the raw reads through to differential expression and pathway analysis using Bioconductor packages (Huber et al. 2015). The article gives a complete analysis of RNA-seq data that were collected to study the effects of pregnancy and lactation on the luminal cell lineage in the mouse mammary gland (Fu et al. 2015). The pipeline uses the Rsubread package (Liao et al. 2013) for mapping reads and assigning them to genes and the edgeR package (Chen et al. 2025) for statistical analyses.
A journal article version of this workflow, including referee reports and responses, is available at F1000Research (Chen et al. 2016). Note that in edgeR version 4 (first released in Bioconductor, October 2024), low-count modelling has been further improved by replacing the saddle-point approximation with bias-corrected deviances (Chen et al. 2025). The workflow results shown here were generated with edgeR v4 and may therefore vary slightly from those reported in the journal article.
RNA-seq analysis involves a number of steps, including read alignment, read summarization, differential expression and pathway analysis. Here we use the Subread aligner (Liao et al. 2013) for mapping and featureCounts (Liao et al. 2014) for assigning reads to genes. As well as being fast and efficient, these algorithms have the advantage of having native implementations as R functions in the Rsubread package. This means that the entire analysis can be conducted efficiently within the R environment.
The workflow uses edgeR’s quasi-likelihood pipeline (edgeR-quasi) for differential expression. This statistical methodology uses negative binomial generalized linear models (McCarthy et al. 2012) but with F-tests instead of likelihood ratio tests (Lund et al. 2012). This method provides stricter error rate control than other negative binomial based pipelines, including the traditional edgeR pipelines (Robinson and Smyth 2008, 2007; McCarthy et al. 2012) or DESeq2 (Love et al. 2014). The edgeR-quasi pipeline is based on similar statistical methodology to that of the QuasiSeq package (Lund et al. 2012), which has performed well in third-party comparisons (Burden et al. 2014). Compared to QuasiSeq, the edgeR functions offer speed improvements and some additional statistical refinements (Lun, Chen, et al. 2016). The RNA-seq pipelines of the limma package also offer excellent error rate control (Law et al. 2014; Ritchie et al. 2015). While the limma pipelines are recommended for large-scale datasets, because of their speed and flexibility, the edgeR-quasi pipeline gives better performance in low-count situations (Lun and Smyth 2014, 2015; Chen et al. 2025). For the data analyzed here, the edgeR-quasi, limma-voom and limma-trend pipelines are all equally suitable and give similar results.
The analysis approach illustrated in this article can be applied to any RNA-seq study that includes some replication, but it is especially appropriate for designed experiments with multiple treatment factors and with small numbers of biological replicates. The approach assumes that RNA samples have been extracted from cells of interest under two or more treatment conditions, that RNA-seq profiling has been applied to each RNA sample and that there are independent biological replicates for at least one of the treatment conditions. The Rsubread part of the workflow takes FASTQ files of raw sequence reads as input, while the edgeR part of the pipeline takes a matrix of genewise read counts as input.
Description of the biological experiment
This workflow demonstrates a complete bioinformatics analysis of an RNA-seq study that is available from the GEO repository as series GSE60450. The RNA-seq data were collected to study the lineage of luminal cells in the mouse mammary gland and in particular how the expression profiles of the members of the lineage change upon pregnancy and lactation (Fu et al. 2015). Specifically, the study examined the expression profiles of basal stem-cell enriched cells (B) and committed luminal cells (L) in the mammary glands of virgin, pregnant and lactating mice. There are therefore six groups of RNA samples, one for each combination of cell type and mouse status. Two biological replicates were collected for each group.
This study used an Illumina Hiseq sequencer to generate about 30 million 100bp single-end reads for each sample. Subread version 1.4.4 (http://subread.sourceforge.net) was used to align the reads to the mouse mm10 genome and featureCounts was used to assign reads to Entrez Genes using RefSeq gene annotation. The FASTQ files containing the raw sequence reads were deposited to the Sequence Read Archive (SRA) repository and the read counts were deposited to GEO.
This experimental design is summarized in the table below, where the
basal and luminal cell types are abbreviated with B and
L respectively. The GEO and SRA identifiers for each RNA
sample are also shown:
> library(RnaSeqGeneEdgeRQL)
> targetsFile <- system.file("extdata", "targets.txt",
+ package="RnaSeqGeneEdgeRQL")
> targets <- read.delim(targetsFile, stringsAsFactors=FALSE)
> targets GEO SRA CellType Status
MCL1.DG GSM1480297 SRR1552450 B virgin
MCL1.DH GSM1480298 SRR1552451 B virgin
MCL1.DI GSM1480299 SRR1552452 B pregnant
MCL1.DJ GSM1480300 SRR1552453 B pregnant
MCL1.DK GSM1480301 SRR1552454 B lactating
MCL1.DL GSM1480302 SRR1552455 B lactating
MCL1.LA GSM1480291 SRR1552444 L virgin
MCL1.LB GSM1480292 SRR1552445 L virgin
MCL1.LC GSM1480293 SRR1552446 L pregnant
MCL1.LD GSM1480294 SRR1552447 L pregnant
MCL1.LE GSM1480295 SRR1552448 L lactating
MCL1.LF GSM1480296 SRR1552449 L lactatingThe experiment can be viewed as a one-way layout with six groups. For later use, we combine the treatment factors into a single grouping factor:
group
B.lactating B.pregnant B.virgin L.lactating L.pregnant L.virgin
2 2 2 2 2 2 Note. This study isolated carefully sorted cell populations obtained from genetically indentical mice under controlled laboratory conditions. As will be shown during the analysis below, the relatively small sample size (\(n=2\) in each group) here is justified by the low background variability and by the fact that the expression profiles of the different groups are distinctly different. Studies with more variability (for example on human subjects) or with smaller effect sizes may well require more biological replicates for reliable results.
Preliminary analysis
Downloading the read counts
Readers wishing to reproduce the analysis presented in this article can either download the matrix of read counts from GEO or recreate the read count matrix from the raw sequence counts. We will present first the analysis using the downloaded matrix of counts. At the end of this article we will present the R commands needed to recreate this matrix.
The following commands download the genewise read counts for the GEO
series GSE60450. The zipped tab-delimited text file
GSE60450_Lactation-GenewiseCounts.txt.gz will be downloaded
to the working R directory:
> FileURL <- paste(
+ "http://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE60450",
+ "format=file",
+ "file=GSE60450_Lactation-GenewiseCounts.txt.gz",
+ sep="&")
> download.file(FileURL, "GSE60450_Lactation-GenewiseCounts.txt.gz")The counts can then be read into a data.frame in R:
> GenewiseCounts <- read.delim("GSE60450_Lactation-GenewiseCounts.txt.gz",
+ row.names="EntrezGeneID")
> colnames(GenewiseCounts) <- substring(colnames(GenewiseCounts),1,7)
> dim(GenewiseCounts)[1] 27179 13 Length MCL1.DG MCL1.DH MCL1.DI MCL1.DJ MCL1.DK MCL1.DL MCL1.LA
497097 3634 438 300 65 237 354 287 0
100503874 3259 1 0 1 1 0 4 0
100038431 1634 0 0 0 0 0 0 0
19888 9747 1 1 0 0 0 0 10
20671 3130 106 182 82 105 43 82 16
27395 4203 309 234 337 300 290 270 560
MCL1.LB MCL1.LC MCL1.LD MCL1.LE MCL1.LF
497097 0 0 0 0 0
100503874 0 0 0 0 0
100038431 0 0 0 0 0
19888 3 10 2 0 0
20671 25 18 8 3 10
27395 464 489 328 307 342The row names of GenewiseCounts are the Entrez Gene
Identifiers. The first column contains the length of each gene, being
the total number of bases in exons and UTRs for that gene. The remaining
12 columns contain read counts and correspond to rows of
targets.
The edgeR package stores data in a simple list-based data
object called a DGEList. This object is easy to use as it
can be manipulated like an ordinary list in R, and it can also be
subsetted like a matrix. The main components of a DGEList
object are a matrix of read counts, sample information in the
data.frame format and optional gene annotation. We enter
the counts into a DGEList object using the function
DGEList in edgeR:
> library(edgeR)
> y <- DGEList(GenewiseCounts[,-1], group=group,
+ genes=GenewiseCounts[,1,drop=FALSE])
> options(digits=3)
> y$samples group lib.size norm.factors
MCL1.DG B.virgin 23227641 1
MCL1.DH B.virgin 21777891 1
MCL1.DI B.pregnant 24100765 1
MCL1.DJ B.pregnant 22665371 1
MCL1.DK B.lactating 21529331 1
MCL1.DL B.lactating 20015386 1
MCL1.LA L.virgin 20392113 1
MCL1.LB L.virgin 21708152 1
MCL1.LC L.pregnant 22241607 1
MCL1.LD L.pregnant 21988240 1
MCL1.LE L.lactating 24723827 1
MCL1.LF L.lactating 24657293 1Adding gene annotation
The Entrez Gene Ids link to gene information in the NCBI database.
The org.Mm.eg.db package can be used to complement the gene
annotation information. Here, a column of gene symbols is added to
y$genes:
> library(org.Mm.eg.db)
> y$genes$Symbol <- mapIds(org.Mm.eg.db, rownames(y),
+ keytype="ENTREZID", column="SYMBOL")
> head(y$genes) Length Symbol
497097 3634 Xkr4
100503874 3259 Gm19938
100038431 1634 Gm10568
19888 9747 Rp1
20671 3130 Sox17
27395 4203 Mrpl15Entrez Ids that no longer have official gene symbols are dropped from
the analysis. The whole DGEList object, including
annotation as well as counts, can be subsetted by rows as if it was a
matrix:
[1] 25849 12Design matrix
Linear modeling and differential expression analysis in edgeR requires a design matrix to be specified. The design matrix records which treatment conditions were applied to each samples, and it also defines how the experimental effects are parametrized in the linear models. The experimental design for this study can be viewed as a one-way layout and the design matrix can be constructed in a simple and intuitive way by:
B.lactating B.pregnant B.virgin L.lactating L.pregnant L.virgin
1 0 0 1 0 0 0
2 0 0 1 0 0 0
3 0 1 0 0 0 0
4 0 1 0 0 0 0
5 1 0 0 0 0 0
6 1 0 0 0 0 0
7 0 0 0 0 0 1
8 0 0 0 0 0 1
9 0 0 0 0 1 0
10 0 0 0 0 1 0
11 0 0 0 1 0 0
12 0 0 0 1 0 0
attr(,"assign")
[1] 1 1 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"This design matrix simply links each group to the samples that belong to it. Each row of the design matrix corresponds to a sample whereas each column represents a coefficient corresponding to one of the six groups.
Filtering to remove low counts
Genes that have very low counts across all the libraries should be removed prior to downstream analysis. This is justified on both biological and statistical grounds. From biological point of view, a gene must be expressed at some minimal level before it is likely to be translated into a protein or to be considered biologically important. From a statistical point of view, genes with consistently low counts are very unlikely be assessed as significantly DE because low counts do not provide enough statistical evidence for a reliable judgement to be made. Such genes can therefore be removed from the analysis without any loss of information.
As a rule of thumb, we require that a gene have a count of at least 10 or so in at least some libraries before it is considered to be expressed in the study. In practice, the filtering is actually based on count-per-million (CPM) values so as to avoid favoring genes that are expressed in larger libraries over those expressed in smaller libraries. In edgeR, the filtering can be accomplished using a function that takes into account the library sizes and the experimental design:
keep
FALSE TRUE
10184 15665 The DGEList object is then subsetted to retain only the
non-filtered genes:
The option keep.lib.sizes=FALSE causes the library sizes
to be recomputed after the filtering. This is generally recommended,
although the effect on the downstream analysis is usually small.
The filterByExpr function attempts to keep the maximum
number of interesting genes in the analysis, but other sensible
filtering criteria are also possible. For example
keep <- rowSums(y$counts) > 50 is a very simple
criterion that would keep genes with a total read count of more than 50.
This would give similar downstream results for this dataset to the
filtering actually used.
Another sensible filter rule would be to compute the average-log-CPM for each gene, and to choose a cutoff value heuristically:
Histogram of average log2 CPM.
Again this might give similar downstream results for this data to the filtering actually used.
Whatever the filtering rule, it should be independent of the information in the targets file. It should not make any reference to which RNA libraries belong to which group, because doing so would bias the subsequent differential expression analysis.
Library size normalization
Normalization by trimmed mean of M values (TMM) (Robinson and Oshlack 2010) is performed by
using the calcNormFactors function, which returns the
DGEList argument with only the norm.factors
changed. It calculates a set of normalization factors, one for each
sample, to eliminate composition biases between libraries. The product
of these factors and the library sizes defines the effective library
size, which replaces the original library size in all downstream
analyses.
calcNormFactors has been renamed to normLibSizes group lib.size norm.factors
MCL1.DG B.virgin 23064696 1.238
MCL1.DH B.virgin 21606688 1.214
MCL1.DI B.pregnant 23899772 1.126
MCL1.DJ B.pregnant 22473268 1.070
MCL1.DK B.lactating 21368942 1.038
MCL1.DL B.lactating 19872559 1.088
MCL1.LA L.virgin 20186601 1.371
MCL1.LB L.virgin 21482079 1.368
MCL1.LC L.pregnant 22051964 1.005
MCL1.LD L.pregnant 21818098 0.923
MCL1.LE L.lactating 24628916 0.528
MCL1.LF L.lactating 24571420 0.534The normalization factors of all the libraries multiply to unity. A normalization factor below one indicates that a small number of high count genes are monopolizing the sequencing, causing the counts for other genes to be lower than would be usual given the library size. As a result, the effective library size will be scaled down for that sample. Here we see that the luminal-lactating samples have low normalization factors. This is a sign that these samples contain a number of very highly upregulated genes.
Note. In general, we find TMM normalization to be satisfactory for almost all well-designed bulk mRNA gene expression experiments, although TMMwsp normalization may be preferable if the proportion of zero counts is unusually large. Specialized normalization methods are generally need for single-cell RNA-seq (Lun, Bach, et al. 2016) and we do not recommend TMM in that context. Another, less common, type of study requiring special treatment is that with global differential expression, with more than half of the genome differentially expressed between experimental conditions in the same direction (Wu et al. 2013). Global differential expression should generally be avoided in well designed experiments. When it can’t be avoided, then some normalization reference such as spike-ins needs to be built into the experiment for reliable normalization to be done (Risso et al. 2014).
Exploring differences between libraries
The RNA samples can be clustered in two dimensions using multi-dimensional scaling (MDS) plots. This is both an analysis step and a quality control step to explore the overall differences between the expression profiles of the different samples. Here we decorate the MDS plot to indicate the cell groups:
> pch <- c(0,1,2,15,16,17)
> colors <- rep(c("darkgreen", "red", "blue"), 2)
> plotMDS(y, col=colors[group], pch=pch[group])
> legend("topleft", legend=levels(group), pch=pch, col=colors, ncol=2)
MDS plot showing distances between expression profiles. This provides a type of unsupervised clustering of the samples. In this case, the first dimension separates samples by cell type while the second dimension corresponds to mouse status.
In the MDS plot, the distance between each pair of samples can be interpreted as the leading log-fold change between the samples for the genes that best distinguish that pair of samples. By default, leading fold-change is defined as the root-mean-square of the largest 500 log2-fold changes between that pair of samples. The figure above shows that replicate samples from the same group cluster together while samples from different groups are well separated. In other words, differences between groups are much larger than those within groups, meaning that there are likely to be statistically significant differences between the groups. The distance between basal cells on the left and luminal cells on the right is about six units on the x-axis, corresponding to a leading fold change of about 64-fold between the two cell types. The differences between the virgin, pregnant and lactating expression profiles appear to be magnified in luminal cells compared to basal.
The expression profiles of individual samples can be explored more closely with mean-difference (MD) plots. An MD plot visualizes the library size-adjusted log-fold change between two libraries (the difference) against the average log-expression across those libraries (the mean). The following command produces an MD plot that compares sample 1 to an artificial reference library constructed from the average of all the other samples:
Mean-difference plot of log2-expression in sample 1 versus the average log2-expression across all other samples. Each point represents a gene, and the red line indicates a log-ratio of zero. The majority of points cluster around the red line.
The bulk of the genes are centered around the line of zero log-fold change. The diagonal lines in the lower left of the plot correspond to genes with counts of 0, 1, 2 and so on in the first sample.
It is good practice to make MD plots for all the samples as a quality check. We now look at one of the luminal-lactating samples that was observed have a low normalization factor:
Mean-difference plot of log2-expression in sample 11 versus the average log2-expression across all other samples. The plot shows a number of genes that are both highly expressed and highly up-regulated.
For this sample, the log-ratios show noticeable positive skew, with a number of very highly upregulated genes. In particular, there are a number of points in the upper right of the plot, corresponding to genes that are both highly expressed and highly up-regulated in this sample compared to others. These genes explain why the normalization factor for this sample is well below one. By contrast, the log-ratios for sample 1 were somewhat negatively skewed, corresponding to a normalization factor above one.
Dispersion estimation
edgeR uses the negative binomial (NB) distribution to model
the read counts for each gene in each sample. The dispersion parameter
of the NB distribution accounts for variability between biological
replicates (McCarthy et al. 2012).
edgeR estimates an empirical Bayes moderated dispersion for
each individual gene. It also estimates a common dispersion, which is a
global dispersion estimate averaged over all genes, and a trended
dispersion where the dispersion of a gene is predicted from its
abundance. Dispersion estimates are most easily obtained from the
estimateDisp function:
This returns a DGEList object with additional components
(common.dispersion, trended.dispersion and
tagwise.dispersion) added to hold the estimated
dispersions. Here robust=TRUE has been used to protect the
empirical Bayes estimates against the possibility of outlier genes with
exceptionally large or small individual dispersions (Phipson et al. 2016).
The dispersion estimates can be visualized with
plotBCV:
Scatterplot of the biological coefficient of variation (BCV) against the average abundance of each gene. The plot shows the square-root estimates of the common, trended and tagwise NB dispersions.
The vertical axis of the plotBCV plot shows square-root
dispersion, also known as biological coefficient of variation
(BCV) (McCarthy et al. 2012). For RNA-seq
studies, the NB dispersions tend to be higher for genes with very low
counts. The dispersion trend tends to decrease smoothly with abundance
and to asymptotic to a constant value for genes with larger counts. From
our past experience, the asymptotic value for the BCV tends to be in
range from 0.05 to 0.2 for genetically identical mice or cell lines,
whereas somewhat larger values (\(>0.3\)) are observed for human
subjects.
The NB model can be extended with quasi-likelihood (QL) methods to account for gene-specific variability from both biological and technical sources (Lund et al. 2012; Lun, Chen, et al. 2016). Under the QL framework, the NB dispersion trend is used to describe the overall biological variability across all genes, and gene-specific variability above and below the overall level is picked up by the QL dispersion. In the QL approach, the individual (tagwise) NB dispersions are not used.
The estimation of QL dispersions is performed using the
glmQLFit function:
B.lactating B.pregnant B.virgin L.lactating L.pregnant L.virgin
497097 -11.13 -12.01 -11.22 -19.0 -19.03 -19.0
20671 -12.76 -12.51 -12.15 -14.5 -14.29 -14.1
27395 -11.27 -11.29 -11.52 -10.6 -10.86 -10.9
18777 -10.15 -10.20 -10.76 -10.1 -10.38 -10.4
21399 -9.89 -9.73 -9.78 -10.2 -9.96 -10.0
58175 -16.15 -14.84 -15.98 -13.3 -12.28 -12.0This returns a DGEGLM object with the estimated values
of the GLM coefficients for each gene. It also contains a number of
empirical Bayes (EB) statistics including the QL dispersion trend, the
squeezed QL dispersion estimates and the prior degrees of freedom (df).
The QL dispersions can be visualized by plotQLDisp:
A plot of the quarter-root QL dispersion against the average abundance
of each gene. Estimates are shown for the raw (before EB moderation),
trended and squeezed (after EB moderation) dispersions. Note that the QL
dispersions and trend shown here are relative to the NB dispersion trend
show in the previous figure produced using plotBCV().
The QL functions moderate the genewise QL dispersion estimates in the same way that the limma package moderates variances (Smyth 2004). The raw QL dispersion estimates are squeezed towards a global trend, and this moderation reduces the uncertainty of the estimates and improves testing power. The extent of the squeezing is governed by the value of the prior df estimated from the data. Large prior df estimates indicate that the QL dispersions are less variable between genes, meaning that strong EB moderation should be performed. Smaller prior df estimates indicate that the true unknown dispersions are highly variable, so weaker moderation towards the trend is appropriate.
Min. 1st Qu. Median Mean 3rd Qu. Max.
6.32 6.32 6.32 6.32 6.32 6.32 Setting robust=TRUE in glmQLFit is usually
recommended (Phipson et al. 2016). This
allows gene-specific prior df estimates, with lower values for outlier
genes and higher values for the main body of genes. This reduces the
chance of getting false positives from genes with extremely high or low
raw dispersions, while at the same time increasing statistical power to
detect differential expression for the main body of genes.
Differential expression analysis
Testing for differential expression
The next step is to test for differential expression between the
experimental groups. One of the most interesting comparisons is that
between the basal pregnant and lactating groups. The contrast
corresponding to any specified comparison can be constructed
conveniently using the makeContrasts function:
In subsequent results, a positive \(\log_2\)-fold-change (logFC) will indicate a gene up-regulated in lactating mice relative to pregnant, whereas a negative logFC will indicate a gene more highly expressed in pregnant mice.
We will use QL F-tests instead of the more usual likelihood ratio tests (LRT) as they give stricter error rate control by accounting for the uncertainty in dispersion estimation:
The top DE genes can be viewed with topTags:
Coefficient: 1*B.lactating -1*B.pregnant
Length Symbol logFC logCPM F PValue FDR
211577 2006 Mrgprf 5.15 2.75 472 7.94e-13 1.24e-08
140474 11281 Muc4 -7.18 6.06 338 4.72e-11 3.00e-07
12992 765 Csn1s2b 6.08 10.19 432 5.75e-11 3.00e-07
24117 2242 Wif1 -1.82 6.77 287 1.39e-10 5.46e-07
21953 667 Tnni2 5.75 3.87 358 1.77e-10 5.53e-07
226101 7094 Myof 2.32 6.45 311 4.11e-10 8.71e-07
381290 8292 Atp2b4 2.14 6.15 309 4.23e-10 8.71e-07
12740 1812 Cldn4 -5.32 9.88 301 4.95e-10 8.71e-07
14621 1445 Gjb4 -4.54 5.55 236 5.00e-10 8.71e-07
78896 888 Ecrg4 -2.81 6.69 229 6.08e-10 9.53e-07In order to control the false discovery rate (FDR), multiple testing correction is performed using the Benjamini-Hochberg method. The top DE gene Csn1s2b has a large positive logFC, showing that it is far more highly expressed in the basal cells of lactating than pregnant mice. This gene is indeed known to be a major source of protein in milk.
The total number of DE genes identified at an FDR of 5% can be shown
with decideTests. There are in fact more than 5000 DE genes
in this comparison:
1*B.lactating -1*B.pregnant
Down 2797
NotSig 10303
Up 2565The magnitude of the differential expression changes can be visualized with a fitted model MD plot:
MD plot showing the log-fold change and average abundance of each gene. Significantly up and down DE genes are highlighted in red and blue, respectively.
The logFC for each gene is plotted against the average abundance in
log2-CPM, i.e., logCPM in the table above. Genes that are
significantly DE are highlighted.
Differential expression above a fold-change threshold
glmQLFTest identifies differential expression based on
statistical significance regardless of how small the difference might
be. For some purposes we might be interested only in genes with
reasonably large expression changes. The above analysis found more than
5000 DE genes between the basal pregnant and lactating groups. With such
a large number of DE genes, it makes sense to narrow down the list to
genes that are more biologically meaningful.
A commonly used approach is to apply FDR and logFC cutoffs
simultaneously. However this tends to favor lowly expressed genes, and
also fails to control the FDR correctly. A better and more rigorous
approach is to modify the statistical test so as to detect expression
changes greater than a specified threshold. In edgeR, this can
be done using the glmTreat function. This function is
analogous to the TREAT method for microarrays (McCarthy and Smyth 2009) but is adapted to the
NB framework. Here we test whether the differential expression fold
changes are significantly greater than 1.5, that is, whether the logFCs
are significantly greater than \(\log_2(1.5)\):
Coefficient: 1*B.lactating -1*B.pregnant
Length Symbol logFC unshrunk.logFC logCPM PValue FDR
211577 2006 Mrgprf 5.15 5.15 2.75 1.38e-12 2.17e-08
12992 765 Csn1s2b 6.08 6.09 10.19 7.61e-11 4.23e-07
140474 11281 Muc4 -7.18 -7.34 6.06 8.10e-11 4.23e-07
21953 667 Tnni2 5.75 5.76 3.87 2.33e-10 9.13e-07
12740 1812 Cldn4 -5.32 -5.33 9.88 6.78e-10 2.12e-06
14621 1445 Gjb4 -4.54 -4.54 5.55 8.10e-10 2.12e-06
226101 7094 Myof 2.32 2.32 6.45 1.90e-09 3.40e-06
14620 1922 Gjb3 -3.58 -3.58 5.51 1.90e-09 3.40e-06
78896 888 Ecrg4 -2.81 -2.81 6.69 1.95e-09 3.40e-06
24117 2242 Wif1 -1.82 -1.82 6.77 2.18e-09 3.41e-06Note that the argument is an abbreviation for ``log-fold-change”. About 1100 genes are detected as DE with a FC significantly above 1.5 at an FDR cut-off of 5%.
1*B.lactating -1*B.pregnant
Down 702
NotSig 14566
Up 397The p-values from glmTreat are larger than those from
glmQLFTest, and the number of significantly DE genes is
fewer, because it is testing an interval null hypothesis and requires
stronger evidence for differential expression than does a conventional
test. It provides greater specificity for identifying the most important
genes with large fold changes.
The test results can be visualized in an MD plot:
MD plot showing the log-fold change and average abundance of each gene. Genes with fold-changes significantly greater than 1.5 are highlighted.
The method evaluates variability as well as the magnitude of change of expression values and therefore is not equivalent to a simple fold change cutoff. Nevertheless, all the statistically significant expression changes have logFC greater than 0.8 and almost all (97%) are greater than 0.9. These values compare to the threshold value of \(\log_2(1.5) = 0.58\). In general, an estimated logFC must exceed the TREAT threshold by a number of standard errors for it to be called significant. In other words, the whole confidence interval for the logFC must clear the threshold rather than just the estimated value itself. It is better to interpret the threshold as the FC below which we are definitely not interested in the gene rather than the FC above which we are interested in the gene.
The value of the FC threshold can be varied depending on the dataset.
In the presence of a huge number of DE genes, a relatively large FC
threshold may be appropriate to narrow down the search to genes of
interest. In the absence of DE genes, on the other hand, a small or even
no FC threshold shall be used. If the threshold level is set to zero,
then glmTreat becomes equivalent to glmQLFTest
in the workflow shown here.
In general, using glmTreat to reduce the number of DE
genes is better than simply reducing the FDR cutoff, because
glmTreat prioritizes genes with larger changes that are
likely to be more biologically significant. glmTreat can
also be used with edgeR pipelines other than quasi-likelihood,
although we don’t demonstrate that here.
Heat map clustering
Heatmaps are a popular way to display differential expression results
for publication purposes. To create a heatmap, we first convert the read
counts into log2-counts-per-million (logCPM) values. This can be done
with the cpm function:
> logCPM <- cpm(y, prior.count=2, log=TRUE)
> rownames(logCPM) <- y$genes$Symbol
> colnames(logCPM) <- paste(y$samples$group, 1:2, sep="-")The introduction of prior.count is to avoid undefined
values and to reduce the variability of the logCPM values for genes with
low counts. Larger values for prior.count shrink the logFCs
for low count genes towards zero.
We will create a heatmap to visualize the top 30 DE genes according
to the TREAT test between B.lactating and
B.pregnant. The advantage of a heatmap is that it can
display the expression pattern of the genes across all the samples.
Visualization of the results is aided by clustering together genes that
have correlated expression patterns. First we select the logCPM values
for the 30 top genes:
A heat map can then be produced by the coolmap
function:
Heat map across all the samples using the top 30 most DE genes between the basal lactating group and the basal pregnancy group.
By default, coolmap clusters genes by correlation
(highly correlated genes are closest) and clusters samples based on
Euclidean distance between the expression values. The coolmap has placed
the basal samples on the right and the luminal samples on the left. In
each case, the virgin and pregnant samples are more like each other than
the lactacting samples. As expected, the two replicate samples from each
group are clustered together.
There a number of other R functions that can create heatmaps, but
coolmap is the simplest because the default clustering
metrics are chosen to be appropriate for expression data.
Analysis of deviance
The differential expression analysis comparing two groups can be easily extended to comparisons between three or more groups. This is done by creating a matrix of independent contrasts. In this manner, users can perform a one-way analysis of deviance (ANODEV) for each gene (Dunn and Smyth 2018).
Suppose we want to compare the three groups in the luminal population, i.e., virgin, pregnant and lactating. An appropriate contrast matrix can be created as shown below, to make pairwise comparisons between all three groups:
> con <- makeContrasts(
+ L.PvsL = L.pregnant - L.lactating,
+ L.VvsL = L.virgin - L.lactating,
+ L.VvsP = L.virgin - L.pregnant, levels=design)The QL F-test is then applied to identify genes that are DE between
the three groups. This combines the three pairwise comparisons into a
single F-statistic and p-value. The top set of significant genes can be
displayed with topTags:
Warning in colnames(design)[coef] <- cn: number of items to replace is not a
multiple of replacement lengthCoefficient: LR test on 2 degrees of freedom
Length Symbol logFC.L.PvsL logFC.L.VvsL logFC.L.VvsP logCPM F
19242 2021 Ptn -1.54 7.26 8.80 7.98 2437
13645 4757 Egf -5.36 -7.22 -1.86 3.68 1198
12992 765 Csn1s2b -8.55 -11.36 -2.81 10.19 1147
11941 7050 Atp2b2 -7.37 -10.56 -3.19 6.61 1127
83492 2382 Gsdmc -4.84 -11.35 -6.51 3.17 454
52150 4089 Kcnk6 -2.42 -7.00 -4.58 5.92 984
15439 1345 Hp 1.08 5.42 4.34 4.94 994
14183 5346 Fgfr2 -1.15 3.95 5.10 7.39 949
17068 691 Ly6d 3.42 9.24 5.82 4.69 872
20856 1793 Stc2 -1.81 3.20 5.01 6.11 888
PValue FDR
19242 9.76e-17 1.53e-12
13645 8.12e-15 4.87e-11
12992 9.97e-15 4.87e-11
11941 1.24e-14 4.87e-11
83492 2.14e-14 5.78e-11
52150 2.56e-14 5.78e-11
15439 2.58e-14 5.78e-11
14183 3.18e-14 6.22e-11
17068 4.39e-14 7.01e-11
20856 4.78e-14 7.01e-11Note that the three contrasts of pairwise comparisons are linearly
dependent. Constructing the contrast matrix with any two of the
contrasts would be sufficient for an ANODEV test. If the contrast matrix
contains all three possible pairwise comparisons, then only the log-fold
changes of the first two contrasts are shown in the output of
topTags.
Complicated contrasts
The flexibility of the GLM framework makes it possible to specify arbitrary contrasts for differential expression tests. Suppose we are interested in testing whether the change in expression between lactating and pregnant mice is the same for basal cells as it is for luminal cells. In statistical terminology, this is the interaction effect between mouse status and cell type. The contrast corresponding to this testing hypothesis can be made as follows.
Then the QL F-test is conducted to identify genes that are DE under
this contrast. The top set of DE genes are viewed with
topTags.
Coefficient: -1*B.lactating 1*B.pregnant 1*L.lactating -1*L.pregnant
Length Symbol logFC logCPM F PValue FDR
19041 6277 Ppl -4.62 6.97 513 2.05e-11 3.21e-07
231991 2873 Creb5 -5.61 4.88 427 6.15e-11 3.45e-07
20512 4206 Slc1a3 5.03 3.66 425 6.60e-11 3.45e-07
67547 3707 Slc39a8 6.19 5.08 378 1.28e-10 4.61e-07
217294 1952 BC006965 -3.88 4.68 369 1.47e-10 4.61e-07
13358 1678 Slc25a1 3.47 7.50 351 1.98e-10 4.82e-07
14598 2022 Ggt1 3.18 6.38 346 2.15e-10 4.82e-07
192166 4558 Sardh 2.92 5.11 334 2.69e-10 5.05e-07
19659 2628 Rbp1 -4.39 6.83 329 2.90e-10 5.05e-07
20440 4836 St6gal1 3.35 7.61 299 5.13e-10 8.04e-07Pathway analysis
Gene ontology analysis
We now consider the problem of interpreting the differential expression results in terms of higher order biological processes or molecular pathways. One of the most common used resources is gene ontology (GO) databases, which annotate genes according to a dictionary of annotation terms. A simple and often effective way to interpret the list of DE genes is to count the number of DE genes that are annotated with each possible GO term. GO terms that occur frequently in the list of DE genes are said to be over-represented or enriched.
In edgeR, GO analyses can be conveniently conducted using
the goana function. Here were apply goana to
the output of the TREAT analysis comparing B.lactating to
B.pregant. The top most significantly enriched GO terms can
be viewed with topGO.
Term Ont N Up Down P.Up P.Down
GO:0000070 mitotic sister chromatid segregation BP 184 1 52 0.991 1.16e-27
GO:1903047 mitotic cell cycle process BP 725 8 105 0.998 1.65e-27
GO:0000280 nuclear division BP 393 6 74 0.936 1.42e-26
GO:0140014 mitotic nuclear division BP 264 2 60 0.991 3.50e-26
GO:0022402 cell cycle process BP 1193 16 137 0.999 1.38e-25
GO:0007059 chromosome segregation BP 398 3 73 0.998 1.73e-25
GO:0000819 sister chromatid segregation BP 221 1 54 0.997 2.75e-25
GO:0098813 nuclear chromosome segregation BP 291 2 61 0.995 1.33e-24
GO:0000278 mitotic cell cycle BP 838 9 107 0.999 1.75e-23
GO:0048285 organelle fission BP 453 10 74 0.716 1.23e-22
GO:0007049 cell cycle BP 1487 18 146 1.000 1.57e-20
GO:0051276 chromosome organization BP 501 5 74 0.996 5.48e-20
GO:0051301 cell division BP 589 4 76 1.000 4.69e-17
GO:0051304 chromosome separation BP 81 0 27 1.000 5.54e-17
GO:0051310 metaphase chromosome alignment BP 101 0 29 1.000 3.62e-16The goana function automatically extracts DE genes from
the tr object, and conducts overlap tests for the up- and
down-regulated DE genes separately. By default, an FDR cutoff of 5% is
used when extracting DE genes, but this can be varied. The row names of
the output are the universal identifiers of the GO terms and the
Term column gives the human-readable names of the terms.
The Ont column shows the ontology domain that each GO term
belongs to. The three domains are: biological process (BP), cellular
component (CC) and molecular function (MF). The N column
represents the total number of genes annotated with each GO term. The
Up and Down columns indicate the number of
genes within the GO term that are significantly up- and down-regulated
in this differential expression comparison, respectively. The
P.Up and P.Down columns contain the p-values
for over-representation of the GO term in the up- and down-regulated
genes, respectively. Note that the p-values are not adjusted for
multiple testing—we would usually ignore GO terms with p-values greater
than about \(10^{-5}\).
By default the output table from topGO is sorted by the
minimum of P.Up and P.Down. Other options are
available. For example, topGO(go, sort="up") lists the top
GO terms that are over-represented in the up-regulated genes. The domain
of the enriched GO terms can also be specified by users. For example,
topGO(go, ontology="BP") restricts to the top GO terms
belonging to the biological process domain while
topGO(go, ontology="MF") restricts to molecular function
terms.
The goana function uses the NCBI RefSeq annotation and
requires the use of Entrez Gene Identifiers.
KEGG pathway analysis
Another popular annotation database is the Kyoto Encyclopedia of
Genes and Genomes (KEGG). Much smaller than GO, this is a curated
database of molecular pathways and disease signatures. A KEGG analysis
can be done exactly as for GO, but using the kegga
function:
Pathway N Up Down P.Up P.Down
mmu04110 Cell cycle 153 0 30 1.00e+00 6.38e-12
mmu03008 Ribosome biogenesis in eukaryot... 75 0 19 1.00e+00 5.06e-10
mmu05150 Staphylococcus aureus infection 50 4 11 3.75e-02 1.02e-05
mmu04911 Insulin secretion 62 9 5 2.44e-05 1.44e-01
mmu01232 Nucleotide metabolism 76 1 13 8.59e-01 2.96e-05
mmu04382 Cornified envelope formation 139 13 11 5.31e-05 4.80e-02
mmu04972 Pancreatic secretion 69 9 4 5.82e-05 3.74e-01
mmu00240 Pyrimidine metabolism 51 0 10 1.00e+00 7.38e-05
mmu00100 Steroid biosynthesis 19 5 0 8.84e-05 1.00e+00
mmu04925 Aldosterone synthesis and secre... 77 9 5 1.39e-04 2.63e-01
mmu04820 Cytoskeleton in muscle cells 198 15 11 1.62e-04 2.75e-01
mmu04970 Salivary secretion 64 8 7 2.03e-04 2.39e-02
mmu05323 Rheumatoid arthritis 69 1 11 8.31e-01 2.30e-04
mmu00900 Terpenoid backbone biosynthesis 23 5 0 2.35e-04 1.00e+00
mmu04914 Progesterone-mediated oocyte ma... 82 2 12 6.19e-01 2.77e-04The output from topKEGG is the same as from
topGO except that row names become KEGG pathway IDs,
Term becomes Pathway and there is no
Ont column. Both the GO and KEGG analyses show that the
cell cycle pathway is strongly down-regulated upon lactation in mammary
stem cells.
By default, the kegga function automatically reads the
latest KEGG annotation from the Internet each time it is run. The KEGG
database uses Entrez Gene Ids, and the kegga function
assumes these are available as the row names of tr.
FRY gene set tests
The GO and KEGG analyses shown above are relatively simple analyses that rely on a list of DE genes. The list of DE genes is overlapped with the various GO and KEGG annotation terms. The results will depend on the significance threshold that is used to assess differential expression.
If the aim is to test for particular gene expression signatures or
particular pathways, a more nuanced approach is to conduct a
roast or fry gene set test (Wu et al. 2010). These functions test whether a
set of genes is DE, assessing the whole set of genes as a whole. Gene
set tests consider all the genes in the specified set and do not depend
on any pre-emptive significance cutoff. The set of genes can be chosen
to be representative of any pathway or phenotype of interest.
roast gives p-values using random rotations of the
residual space. In the edgeR context, fry is
generally recommended over roast. fry gives an
accurate analytic approximation to the results that roast
would give, with default settings, if an extremely large number of
rotations was used.
Here, suppose we are interested in two GO terms related to cytokinesis. Each GO term is used to define a set of genes annotated with that term. The names of these terms are shown below:
> library(GO.db)
> cyt.go <- c("GO:0032465", "GO:0000281")
> term <- select(GO.db, keys=cyt.go, columns="TERM")
> term GOID TERM
1 GO:0032465 regulation of cytokinesis
2 GO:0000281 mitotic cytokinesisThe first step is to extract the genes associated with each GO term from the GO database. This produces a list of two components, one for each GO term. Each component is a vector of Entrez Gene IDs for that GO term:
Suppose the comparison of interest is between the virgin and
lactating groups in the basal population. We can use fry to
test whether the cytokinesis GO terms are DE for this comparison:
> B.VvsL <- makeContrasts(B.virgin-B.lactating, levels=design)
> fry(y, index=cyt.go.genes, design=design, contrast=B.VvsL) NGenes Direction PValue FDR PValue.Mixed FDR.Mixed
GO:0032465 79 Up 0.00366 0.00731 8.00e-06 9.41e-06
GO:0000281 92 Up 0.01909 0.01909 9.41e-06 9.41e-06Each row of the output corresponds to a gene set. The
NGenes column provides the number of genes in each set. The
Direction column indicates the net direction of change. The
PValue column gives the two-sided p-value for testing
whether the set is DE as a whole, either up or down. The
PValue.Mixed column gives a p-value for testing whether
genes in the set tend to be DE, without regard to direction. The
PValue column is appropriate when genes in the set are
expected to be co-regulated, all or most changing expression in the same
direction. The PValue.Mixed column is appropriate when
genes in the set are not necessarily co-regulated or may be regulated in
different directions for the contrast in question. FDRs are calculated
from the corresponding p-values across all sets.
The results of a gene set test can be viewed in a barcode plot
produced by the barcodeplot function. Suppose visualization
is performed for the gene set defined by the GO term GO:0032465:
> res <- glmQLFTest(fit, contrast=B.VvsL)
> index <- rownames(fit) %in% cyt.go.genes[[1]]
> barcodeplot(res$table$logFC, index=index, labels=c("B.lactating","B.virgin"),
+ main=cyt.go[1])
Barcode plot showing enrichment of the GO term GO:0032465 in the basal virgin group compared to the basal lactating group. X-axis shows logFC for B.virgin vs B.lactating. Black bars represent genes annotated with the GO term. The worm shows relative enrichment.
In the plot, all genes are ranked from left to right by increasing
log-fold change for the contrast and the genes within the gene set are
represented by vertical bars, forming the barcode-like pattern. The
curve (or worm) above the barcode shows the relative local
enrichment of the bars in each part of the plot. The dotted horizontal
line indicates neutral enrichment; the worm above the dotted line shows
enrichment while the worm below the dotted line shows depletion. In this
particular barcode plot the worm shows enrichment on the right for
positive logFCs, and depletion on the left for negative logFCs. The
conclusion is that genes associated with this GO term tend to be
up-regulated in the basal cells of virgin mice compared to lactating
mice, confirming the result of the fry test above.
Camera gene set enrichment analysis
Finally we demonstrate a gene set enrichment style analysis using the Molecular Signatures Database (MSigDB) (Subramanian et al. 2005). We will use the C2 collection of the MSigDB, which is a collection of over 5000 curated gene sets, each representing the molecular signature of a particular biological process or phenotype. The MSigDB itself is purely human, but the Walter and Eliza Hall Institute (WEHI) maintains a mouse version of the database. We load the mouse version of the C2 collection from the WEHI website:
This will read into Mm.c2 a list of gene sets, each a
vector of NCBI GeneIDs. The GeneIDs can be converted to row numbers
by:
First we compare virgin stem cells to virgin luminal cells:
> BvsL.v <- makeContrasts(B.virgin - L.virgin, levels=design)
> cam <- camera(y, idx, design, contrast=BvsL.v, inter.gene.cor=0.01)
> options(digits=2)
> head(cam,14) NGenes Direction PValue
LIM_MAMMARY_STEM_CELL_UP 723 Up 8.2e-50
LIM_MAMMARY_LUMINAL_MATURE_DN 164 Up 1.5e-25
LIM_MAMMARY_STEM_CELL_DN 555 Down 4.8e-24
FARMER_BREAST_CANCER_CLUSTER_4 47 Up 7.0e-23
REACTOME_NCAM1_INTERACTIONS 46 Up 1.8e-21
NABA_BASEMENT_MEMBRANES 50 Up 8.1e-20
ROZANOV_MMP14_TARGETS_SUBSET 68 Up 1.3e-19
ANASTASSIOU_MULTICANCER_INVASIVENESS_SIGNATURE 110 Up 1.5e-19
HAEGERSTRAND_RESPONSE_TO_IMATINIB 31 Up 3.2e-19
REACTOME_COLLAGEN_BIOSYNTHESIS_AND_MODIFYING_ENZYMES 64 Up 6.9e-19
CLASPER_LYMPHATIC_VESSELS_DURING_METASTASIS_DN 66 Up 1.0e-18
REACTOME_LAMININ_INTERACTIONS 46 Up 1.3e-18
LIM_MAMMARY_LUMINAL_PROGENITOR_UP 81 Down 3.6e-18
REACTOME_ECM_PROTEOGLYCANS 86 Up 9.2e-18
FDR
LIM_MAMMARY_STEM_CELL_UP 4.5e-46
LIM_MAMMARY_LUMINAL_MATURE_DN 4.2e-22
LIM_MAMMARY_STEM_CELL_DN 8.8e-21
FARMER_BREAST_CANCER_CLUSTER_4 9.7e-20
REACTOME_NCAM1_INTERACTIONS 2.0e-18
NABA_BASEMENT_MEMBRANES 7.5e-17
ROZANOV_MMP14_TARGETS_SUBSET 1.0e-16
ANASTASSIOU_MULTICANCER_INVASIVENESS_SIGNATURE 1.1e-16
HAEGERSTRAND_RESPONSE_TO_IMATINIB 1.9e-16
REACTOME_COLLAGEN_BIOSYNTHESIS_AND_MODIFYING_ENZYMES 3.8e-16
CLASPER_LYMPHATIC_VESSELS_DURING_METASTASIS_DN 5.2e-16
REACTOME_LAMININ_INTERACTIONS 5.9e-16
LIM_MAMMARY_LUMINAL_PROGENITOR_UP 1.5e-15
REACTOME_ECM_PROTEOGLYCANS 3.6e-15With a large gene set collection, setting
inter.gene.cor = 0.01 gives a good compromise between
biological interpretability and FDR control. As expected, the mammary
stem cell and mammary luminal cell signatures from Lim et al (Lim et al. 2010) are top-ranked, and in the
expected directions.
We can visualize the top signature, combining the up and down mammary stem cell signatures to make a bi-directional signature set:
> res <- glmQLFTest(fit, contrast=BvsL.v)
> barcodeplot(res$table$logFC,
+ index=idx[["LIM_MAMMARY_STEM_CELL_UP"]],
+ index2=idx[["LIM_MAMMARY_STEM_CELL_DN"]],
+ labels=c("L.virgin","B.virgin"),
+ main="LIM_MAMMARY_STEM_CELL",
+ alpha=1)
Barcode plot showing strong enrichment of mammary stem cell signature in the stem cell vs luminal cell comparison. Red bars show up signature genes, blue bars show down genes. The worms show relative enrichment.
Packages used
This workflow depends on various packages from version 3.3 of the Bioconductor project, running on R version 3.3.0 or higher. The complete list of the packages used for this workflow are shown below:
R version 4.6.1 (2026-06-24)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 26.04.1 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
time zone: Etc/UTC
tzcode source: system (glibc)
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] RnaSeqGeneEdgeRQL_1.37.0 GO.db_3.23.1 org.Mm.eg.db_3.23.0
[4] AnnotationDbi_1.75.2 IRanges_2.47.5 S4Vectors_0.51.10
[7] Biobase_2.73.2 BiocGenerics_0.59.12 generics_0.1.4
[10] gplots_3.3.0 edgeR_4.99.6 limma_3.99.0
[13] BiocStyle_2.41.0
loaded via a namespace (and not attached):
[1] sass_0.4.10 bitops_1.1-0 KernSmooth_2.23-27
[4] gtools_3.9.5 RSQLite_3.53.3 lattice_0.23-1
[7] digest_0.6.39 caTools_1.18.4 evaluate_1.0.5
[10] grid_4.6.1 fastmap_1.2.0 blob_1.3.0
[13] jsonlite_2.0.0 DBI_1.3.0 BiocManager_1.30.27
[16] httr_1.4.9 Biostrings_2.81.9 jquerylib_0.1.4
[19] cli_3.6.6 rlang_1.3.0 crayon_1.5.3
[22] XVector_0.53.0 bit64_4.8.6 cachem_1.1.0
[25] yaml_2.3.12 BiocBaseUtils_1.15.1 tools_4.6.1
[28] memoise_2.0.1 locfit_1.5-9.12 buildtools_1.0.0
[31] vctrs_0.7.3 R6_2.6.1 png_0.1-9
[34] lifecycle_1.0.5 Seqinfo_1.3.2 KEGGREST_1.53.6
[37] bit_4.6.0 pkgconfig_2.0.3 bslib_0.12.0
[40] statmod_1.5.2 xfun_0.61 sys_3.4.3
[43] knitr_1.52 htmltools_0.5.9 rmarkdown_2.32
[46] maketools_1.3.2 compiler_4.6.1 Read alignment and quantification
Download raw sequence files from the SRA
We now revisit the question of recreating the matrix of read counts
from the raw sequence reads. Unlike the above workflow, which works for
any version of R, read alignment requires Unix or Mac OS and, in
practice, a high performance Unix server is recommended. Read alignment
and read counting require only one Bioconductor package,
Rsubread. However the fastq-dump utility from the
SRA Toolkit is also required to convert from SRA to FASTQ format. This
can be downloaded from the NCBI website (http://www.ncbi.nlm.nih.gov/Traces/sra/?view=software)
and installed on any Unix system.
The first task is to download the raw sequence files, which are
stored in SRA format on the SRA repository. The SRA files need to be
unpacked into FASTQ format using the fastq-dump utility.
The following R code makes a system call to fastq-dump to
download each SRA file and convert it to FASTQ format:
for (sra in targets$SRA) system(paste("fastq-dump", sra))The fastq-dump utility automatically downloads the
specified SRA data set from the internet. The above code will produce 12
FASTQ files, in the current working directory, with file names given by
the following vector:
all.fastq <- paste0(targets$SRA, ".fastq")Accuracy of base-calling
Sequencers typically store base-calling quality scores for each read
in the FASTQ files. Rsubread’s qualityScores
function can be used to extract these scores from any particular
file:
QS <- qualityScores("SRR1552444.fastq")The boxplot function provides a compact way to view the quality scores by position across all reads:
boxplot(QS, ylab="Quality score", xlab="Base position",
main="SRR1552444.fastq", cex=0.25, col="orange")
Figure 12: Boxplots of quality scores by base position for the first FASTQ file.
The vertical axis shows the Phred quality score, equal to \(-10\log_{10}(p)\) where \(p\) is the probability of an erroneous call. The maximum possible value is 40, and all values above 10 correspond to extremely small error probabilities. The horizontal axis shows position within a read. The file contains 100bp single-end reads, so the scale is from 1 to 100. The plot displays a compact boxplot at each base position. As is very commonly observed, the quality scores are best in the middle of the reads and decrease slightly towards the start and end of the reads. However the quality remains generally good even near the ends of the reads: the scores would need to be very much lower than this before they would cause problems for the alignment. Similar plots can be made for each of the FASTQ files.
Build a genome index
Before the sequence reads can be aligned, we need to build an index for the GRCm38/mm10 (Dec 2011) build of the mouse genome. Most laboratories that use Rsubread regularly will already have an index file prepared, as this is a once-off operation for each genome release. If you are using Rsubread for mouse for the first time, then the latest mouse genome build can be downloaded from the NCBI webpage (https://www.ncbi.nlm.nih.gov/assembly/GCF_000001635.26). The full URL is:
ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/635/GCF_000001635.26_GRCm38.p6/GCF_000001635.26_GRCm38.p6_genomic.fna.gz(Note that this link is for patch 6 of the RefSeq mm10 build, which is valid at the time of writing in May 2018. The link will change as new patches are released periodically.) An index can then be built by:
library(Rsubread)
buildindex(basename = "mm10",
reference = "GCF_000001635.26_GRCm38.p6_genomic.fna.gz")Aligning reads
The sequence reads can now be aligned to the mouse genome using the
align function:
all.bam <- sub(".fastq", ".bam", all.fastq)
align(index="mm10", readfile1=all.fastq, input_format="FASTQ",
output_file=all.bam)This produces a set of BAM files containing the read alignments for
each RNA library. The mapping proportions can be summarized by the
propmapped function:
> propmapped(all.bam)
Samples NumTotal NumMapped PropMapped
1 SRR1552450.bam 30109290 26577308 0.883
2 SRR1552451.bam 28322351 24794251 0.875
3 SRR1552452.bam 31688348 27937620 0.882
4 SRR1552453.bam 29614284 26074034 0.880
5 SRR1552454.bam 27225012 24381742 0.896
6 SRR1552455.bam 25433157 22813815 0.897
7 SRR1552444.bam 27919481 23927833 0.857
8 SRR1552445.bam 29731031 25487822 0.857
9 SRR1552446.bam 29879070 25500318 0.853
10 SRR1552447.bam 29245388 25187577 0.861
11 SRR1552448.bam 31425424 27326500 0.870
12 SRR1552449.bam 31276061 27204156 0.870Ideally, the proportion of mapped reads should be above 80%. By default, only reads with unique mapping locations are reported by Rsubread as being successfully mapped. Restricting to uniquely mapped reads is recommended, as it avoids spurious signal from non-uniquely mapped reads derived from, e.g., repeat regions.
Quantifying read counts for each gene
The read counts for each gene can be quantified using the
featureCounts function in Rsubread. Conveniently,
the Rsubread package includes inbuilt NCBI RefSeq annotation of
the mouse and human genomes. featureCounts generates a
matrix of read counts for each gene in each sample:
fc <- featureCounts(all.bam, annot.inbuilt="mm10")The output is a simple list, containing the matrix of counts
(counts), a data frame of gene characteristics
(annotation), a vector of file names (targets)
and summary mapping statistics (stat):
> names(fc)
[1] "counts" "annotation" "targets" "stat" The row names of fc$counts are the Entrez gene
identifiers for each gene. The column names are the output file names
from align, which we simplify here for brevity:
> colnames(fc$counts) <- rownames(targets)The first six rows of the counts matrix are shown below.
> head(fc$counts)
MCL1.DG MCL1.DH MCL1.DI MCL1.DJ MCL1.DK MCL1.DL MCL1.LA
497097 438 299 65 237 354 287 0
100503874 1 0 1 1 0 4 0
100038431 0 0 0 0 0 0 0
19888 1 1 0 0 0 0 10
20671 106 181 82 104 43 83 16
27395 309 232 339 290 291 270 558
MCL1.LB MCL1.LC MCL1.LD MCL1.LE MCL1.LF
497097 0 0 0 0 0
100503874 0 0 0 0 0
100038431 0 0 0 0 0
19888 3 10 2 0 0
20671 25 18 8 3 10
27395 468 488 332 312 344Finally, a DGEList object can be assembled by:
y <- DGEList(fc$counts, group=group)
y$genes <- fc$annotation[, "Length", drop=FALSE]Data and software availability
This html document was generated as a vignette of the RnaSeqGeneEdgeRQL package, which can be installed by running
install.packages("BiocManager")
BiocManager::install("RnaSeqGeneEdgeRQL")at the R prompt. Except for the targets file
targets.txt, all data analyzed in the workflow is read
automatically from public websites as part of the code. All software
used is publicly available as part of Bioconductor (Release 3.3 or
later), except for the fastq-dump utility, which can be
downloaded from NCBI website as described in the text. The article
includes the complete code necessary to reproduce the analyses
shown.
The original Bioconductor 3.3 version of this workflow, including a pdf article version and a corresponding Rnw file, is available from http://f1000research.com/articles/5-1438 and http://bioinf.wehi.edu.au/edgeR/F1000Research2016 (Chen et al. 2016).
Grant information
This work was supported by the National Health and Medical Research Council (Fellowship 1058892 and Program 1054618 to G.K.S, Independent Research Institutes Infrastructure Support to the Walter and Eliza Hall Institute) and by a Victorian State Government Operational Infrastructure Support Grant.
Acknowledgments
The authors thank Wei Shi and Yang Liao for advice with Rsubread and Yifang Hu for creating the mouse version of the MSigDB. We also wish to acknowledge the early developers of the quasi-likelihood testing approach. Davis McCarthy authored the first versions of the quasi-likelihood functions in edgeR in February 2011 and Steve Lund and Dan Nettleton worked independently on QuasiSeq around the same time. We thank Steve Lund and Dan Nettleton for a valuable and enjoyable collaboration that lead to Lund et al (Lund et al. 2012).
How to cite
Please cite this workflow article as
Chen, Yunshun, A. T. L. Lun, and G. K. Smyth. 2016. “From Reads to Genes to Pathways: Differential Expression Analysis of RNA-Seq Experiments Using Rsubread and the edgeR Quasi-Likelihood Pipeline.” F1000Research 5:1438.