In this tutorial we use differential expression analysis to identify taxa (e.g., individual OTUs) that have abundances that vary significantly between two groups. In Tutorial 8 we used PERMANOVA to determine if bacterial community composition varies between groups. Thus the two analyses can be thought of as complementing each other, PERMANOVA reveals differences at the community level, and differential expression analysis reveals differences at the OTU level.
DESeq2 is a Bioconductor package. To install DESeq2 we need to install Bioconductor first. Follow these instructions to install Bioconductor and DESeq2.
Setup and reading in the clean data
Let’s start by cleaning up our work space and loading the required packages.
rm(list =ls())gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 580966 31.1 1325185 70.8 NA 669434 35.8
Vcells 1070401 8.2 8388608 64.0 102400 1851726 14.2
library(tidyverse)
Warning: package 'tidyr' was built under R version 4.3.2
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.4.4 ✔ tibble 3.2.1
✔ lubridate 1.9.3 ✔ tidyr 1.3.1
✔ purrr 1.0.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(DESeq2)
Warning: package 'DESeq2' was built under R version 4.3.1
Loading required package: S4Vectors
Warning: package 'S4Vectors' was built under R version 4.3.2
Warning: package 'BiocGenerics' was built under R version 4.3.1
Attaching package: 'BiocGenerics'
The following objects are masked from 'package:lubridate':
intersect, setdiff, union
The following objects are masked from 'package:dplyr':
combine, intersect, setdiff, union
The following objects are masked from 'package:stats':
IQR, mad, sd, var, xtabs
The following objects are masked from 'package:base':
anyDuplicated, aperm, append, as.data.frame, basename, cbind,
colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
table, tapply, union, unique, unsplit, which.max, which.min
Attaching package: 'S4Vectors'
The following objects are masked from 'package:lubridate':
second, second<-
The following objects are masked from 'package:dplyr':
first, rename
The following object is masked from 'package:tidyr':
expand
The following object is masked from 'package:utils':
findMatches
The following objects are masked from 'package:base':
expand.grid, I, unname
Loading required package: IRanges
Warning: package 'IRanges' was built under R version 4.3.1
Attaching package: 'IRanges'
The following object is masked from 'package:lubridate':
%within%
The following objects are masked from 'package:dplyr':
collapse, desc, slice
The following object is masked from 'package:purrr':
reduce
Loading required package: GenomicRanges
Warning: package 'GenomicRanges' was built under R version 4.3.1
Loading required package: GenomeInfoDb
Warning: package 'GenomeInfoDb' was built under R version 4.3.2
Loading required package: SummarizedExperiment
Warning: package 'SummarizedExperiment' was built under R version 4.3.1
Loading required package: MatrixGenerics
Warning: package 'MatrixGenerics' was built under R version 4.3.1
Warning: package 'Biobase' was built under R version 4.3.1
Welcome to Bioconductor
Vignettes contain introductory material; view with
'browseVignettes()'. To cite Bioconductor, see
'citation("Biobase")', and for packages 'citation("pkgname")'.
Attaching package: 'Biobase'
The following object is masked from 'package:MatrixGenerics':
rowMedians
The following objects are masked from 'package:matrixStats':
anyMissing, rowMedians
Now we will read in the clean data. We will not be using ordination coordinates, so we can use the metadata that does not include them.
In the Lake Erie data PERMANOVA tutorial, we determined that the station where a sample was collected significantly explained some of the variability in the bacterial community compositions in the Lake Erie data. In the NMDS plot from that tutorial it appears that the community composition of the samples from the offshore station is different from the composition of the samples from the two nearshore stations.
A DESeq2 analysis compares count data between two groups of samples. We will group samples from the two nearshore stations together and compare samples from the offshore station to the two nearshore stations. To do this we will add a new variable to the metadata called Station_location with levels “nearshore” and “offshore”. An alternative approach to obtaining two groups to compare would be to filter the data so that samples from the offshore station and just one of the nearshore stations were retained.
erie_meta <- erie_meta %>%mutate(Station_location = Station %>%fct_collapse(nearshore =c("nearshore1", "nearshore2")))# if we wanted to filter out the second nearshore station instead we could do this# erie_meta <- erie_meta %>%# filter(Station == "nearshore1" | Station == "offshore")
Pre-filtering the data by removing OTUs with very few reads from the data
At this point it is possible that there are OTUs in our filtered data set that have zero reads. This is due to the fact that we filtered the data to remove some samples (using a minimum sequencing depth and possibly to reduce the dataset to two groups). There is no reason to retain these OTUs in our analysis.
In fact, it is helpful to remove any OTUs with very few reads (not just 0 reads) from the dataset. These OTUs do not impact the results of the analysis. Filtering them out at this point reduces the size of the objects we are working with, requiring less memory.
We will use 10 reads as a threshold, and only keep OTUs with at least 10 reads in the data set.
We’ll start by creating a vector that lists all of the OTUs that have at least 10 reads in the filtered data.
Now we remove any OTUs that have less than 10 reads from the taxonomy table and the shared table by only keeping the ones in the list. At the same time, we will also drop the two extra columns that we won’t use (label and numOtus).
The DESeq function we will use for our analysis expects the data to be in the form of a DESeqDataSet object. One way to generate such an object is using a count matrix. The count matrix must have a row for each OTU and column for each sample. The row names should be the OTU names, and the column names should be the sample IDs. Our data are currently stored in a tibble (not a matrix) and have the opposite orientation.
The DESeq function matches metadata with samples using their positions in the data. Thus we need to make sure that the samples in the metadata are ordered the same as in the shared table. We will accomplish this by joining the two data sets. However, before we do that, we will reduce the metadata to two variables–SampleID and Station_location, which is all we need for the analysis.
erie_sl <- erie_meta %>%select(SampleID,Station_location)erie_share_sl <- erie %>%inner_join(erie_sl, by ="SampleID")
Let’s create the count matrix. We’ll print the first 2 rows to the screen.
Now we can create the DESeqDataSet object. We will use the DESeqDataSetFromMatrix command. The design option indicates the variable that indicates the groups we wish to compare.
Finally, we will identify which level of the Station_location variable should be considered the reference level (e.g., the control). Otherwise, DESeq defaults to using the first factor level, which is determined alphabetically by default. We will arbitrarily choose nearshore as our reference level.
-- replacing outliers and refitting for 155 genes
-- DESeq argument 'minReplicatesForReplace' = 7
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing
Then we use the results function and print some of the results. We format the results as a tibble. The first column of the results table shows the OTUs. The baseMean column shows the average normalized counts across all of the samples.
The log2FoldChange column shows the log2 fold change (lfc). The fold change is ratio of the concentrations (taking into account differences in sequencing depth) of OTUs between the two groups. That is, if the fold change is 4 (a 4-fold change) for a particular OTU, then that otu is estimated as having 4 times the concentration in the treatment group than in the control group. The log2 fold change is simply the log base 2 of the fold change. Thus, for example, if the lfc is 3, that implies the fold change is 8.
Both a p-value and adjustedp-value are reported for each OTU. The null hypothesis being tested is that there is no difference in the concentration of that OTU between the two groups. Recall that the p-value is the probability of obtaining a result at least as extreme as the observed one if the null hypothesis is true. Typically, if you are conducting a single test of significance you define a significance level (\(\alpha\)) ahead of time that indicates what you are willing to accept as a type 1 error rate (the probability of rejecting a true null hypothesis). In biology it is common to use \(\alpha = 0.05\). When you conduct multiple tests of significance, the probability that at least one of the tests will result in rejecting a true null hypothesis can become much larger than \(\alpha\) (especially when conducting hundreds or thousands of tests). The adjusted p-value accounts for this by adjusting the p-value for each test so that the overall false positive rate (or the false discovery rate, depending on the adustment method) is approximately the same as it would be for a single test. The false discovery rate is the proportion of the rejected null hypotheses (false positives + true positives) that are actually false. See, for example Jafari and Ansari-Pour (2018) for more details.
Let’s focus on the first row of the output table. It tells us that Otu00002 is estimated to have \(2^{1.411729}\approx 2.66\) times higher concentration in offshore locations (the treatment) than in nearshore locations (the control). Based on the adjusted p-value we would reject the null hypothesis that the concentrations are the same in both types of locations at a significance level of \(\alpha=0.05\).
Statistical significance vs. biological significance
A statistically significant difference is not always biologically relevant. Although larger differences between groups result in smaller p-values, it is not always the case that small p-values imply large differences between groups. Small differences are often statistically significant when the sample size is large or when the variation within groups is small. To be of biological interest (i.e., to be considered biologically significant) the differences between groups must be large enough to be considered relevant. We may, for example, decide that we are only interested in OTUs that are at least twice as abundant in one location type than the other (i.e., OTUs that have a log2 fold change value that is at least 1 or at most -1).
If we have such a threshold in mind then we can test whether OTUs are significantly above or significantly below the threshold. For example, if we use an log2 fold change threshold of 1 for biological significance, then each null hypothesis states that the concentration of the OTU is between half and twice as large in the treatment as it is in the control.
We can implement this using DESeq by specifying appropriate options for the results function. (By default lcfThreshold = 0).
This time there are only 31 significant results. These are OTUs we should probably pay attention to!
erie_dds_res_lfc1 %>% dplyr::count(padj <0.05)
# A tibble: 3 × 2
`padj < 0.05` n
<lgl> <int>
1 FALSE 611
2 TRUE 31
3 NA 1148
In some cases we may wish to identify OTUs that have fold change values with magnitudes that are significantly below a threshold. These are the OTUs that we would be confident are not changing very much between treatments. We can also use the results function to find these OTUs.