Fluent genomics with plyranges and tximeta

Stuart Lee, Michael Lawrence, Michael Love

Abstract

We construct a simple workflow for fluent genomics data analysis using the R/Bioconductor ecosystem. This involves three core steps: import the data into an appropriate abstraction, model the data with respect to the biological questions of interest, and integrate the results with respect to their underlying genomic coordinates. Here we show how to implement these steps to integrate published RNA-seq and ATAC-seq experiments on macrophage cell lines. Using tximeta, we import RNA-seq transcript quantifications into an analysis-ready data structure, called the SummarizedExperiment, that contains the ranges of the reference transcripts and metadata on their provenance. Using SummarizedExperiments to represent the ATAC-seq and RNA-seq data, we model differentially accessible (DA) chromatin peaks and differentially expressed (DE) genes with existing Bioconductor packages. Using plyranges we then integrate the results to see if there is an enrichment of DA peaks near DE genes by finding overlaps and aggregating over log-fold change thresholds. The combination of these packages and their integration with the Bioconductor ecosystem provide a coherent framework for analysts to iteratively and reproducibly explore their biological data.

1 Introduction

In this workflow, we examine a subset of the RNA-seq and ATAC-seq data from Alasoo et al. (2018), a study that involved treatment of macrophage cell lines from a number of human donors with interferon gamma (IFNg), Salmonella infection, or both treatments combined. Alasoo et al. (2018) examined gene expression and chromatin accessibility in a subset of 86 successfully differentiated induced pluripotent stem cells (iPSC) lines, and compared baseline and response with respect to chromatin accessibility and gene expression at specific quantitative trait loci (QTL). The authors found that many of the stimulus-specific expression QTL were already detectable as chromatin QTL in naive cells, and further hypothesize about the nature and role of transcription factors implicated in the response to stimulus.

We will perform a much simpler analysis than the one found in Alasoo et al. (2018), using their publicly available RNA-seq and ATAC-seq data (ignoring the genotypes). We will examine the effect of IFNg stimulation on gene expression and chromatin accessibility, and look to see if there is an enrichment of differentially accessible (DA) ATAC-seq peaks in the vicinity of differentially expressed (DE) genes. This is plausible, as the transcriptomic response to IFNg stimulation may be mediated through binding of regulatory proteins to accessible regions, and this binding may increase the accessibility of those regions such that it can be detected by ATAC-seq.

An overview of the fluent genomics workflow. First, we import data as a SummarizedExperiment object, which enables interoperability with downstream analysis packages. Then we model our assay data, using the existing Bioconductor packages DESeq2 and limma. We take the results of our models for each assay with respect to their genomic coordinates, and integrate them. First, we compute the overlap between the results of each assay, then aggregate over the combined genomic regions, and finally summarize to compare enrichment for differentially expressed genes to non differentially expressed genes. The final output can be used for downstream visualization or further transformation.

Figure 1.1: An overview of the fluent genomics workflow. First, we import data as a SummarizedExperiment object, which enables interoperability with downstream analysis packages. Then we model our assay data, using the existing Bioconductor packages DESeq2 and limma. We take the results of our models for each assay with respect to their genomic coordinates, and integrate them. First, we compute the overlap between the results of each assay, then aggregate over the combined genomic regions, and finally summarize to compare enrichment for differentially expressed genes to non differentially expressed genes. The final output can be used for downstream visualization or further transformation.

Throughout the workflow (Figure 1.1), we will use existing Bioconductor infrastructure to understand these datasets. In particular, we will emphasize the use of the Bioconductor packages plyranges and tximeta. The plyranges package fluently transforms data tied to genomic ranges using operations like shifting, window construction, overlap detection, etc. It is described by Lee, Cook, and Lawrence (2019) and leverages underlying core Bioconductor infrastructure (Lawrence et al. 2013; Huber et al. 2015) and the tidyverse design principles Wickham et al. (2019).

The tximeta package described by Love et al. (2019) is used to read RNA-seq quantification data into R/Bioconductor, such that the transcript ranges and their provenance are automatically attached to the object containing expression values and differential expression results.

1.1 Experimental Data

The data used in this workflow is available from two packages: the macrophage Bioconductor ExperimentData package and from the workflow package fluentGenomics.

The macrophage package contains RNA-seq quantification from 24 RNA-seq samples, a subset of the RNA-seq samples generated and analyzed by Alasoo et al. (2018). The paired-end reads were quantified using Salmon (Patro et al. 2017), using the Gencode 29 human reference transcripts (Frankish, GENCODE-consoritum, and Flicek 2018). For more details on quantification, and the exact code used, consult the vignette of the macrophage package. The package also contains the Snakemake file that was used to distribute the Salmon quantification jobs on a cluster (Köster and Rahmann 2012).

The fluentGenomics package contains functionality to download and generate a cached SummarizedExperiment object from the normalized ATAC-seq data provided by Alasoo and Gaffney (2017). This object contains all 145 ATAC-seq samples across all experimental conditions as analyzed by Alasoo et al. (2018). The data can be also be downloaded directly from the Zenodo deposition.

The following code loads the path to the cached data file, or if it is not present, will create the cache and generate a SummarizedExperiment using the the BiocFileCache package (Shepherd and Morgan 2019).

We can then read the cached file and assign it to an object called atac.

A precise description of how we obtained this SummarizedExperiment object can be found in section 2.2.

2 Import Data as a SummarizedExperiment

2.1 Using tximeta to import RNA-seq quantification data

First, we specify a directory dir, where the quantification files are stored. You could simply specify this directory with:

where the path is relative to your current R session. However, in this case we have distributed the files in the macrophage package. The relevant directory and associated files can be located using system.file.

Information about the experiment is contained in the coldata.csv file. We leverage the dplyr and readr packages (as part of the tidyverse) to read this file into R (Wickham et al. 2019). We will see later that plyranges extends these packages to accommodate genomic ranges.

## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## Rows: 24 Columns: 13
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (10): names, sample_id, line_id, condition_name, macrophage_harvest, sal...
## dbl  (3): replicate, ng_ul_mean, rna_auto
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## # A tibble: 24 × 5
##    names          id     line   condition   files                               
##    <chr>          <chr>  <fct>  <fct>       <chr>                               
##  1 SAMEA103885102 diku_A diku_1 naive       /home/biocbuild/bbs-3.15-bioc/R/lib…
##  2 SAMEA103885347 diku_B diku_1 IFNg        /home/biocbuild/bbs-3.15-bioc/R/lib…
##  3 SAMEA103885043 diku_C diku_1 SL1344      /home/biocbuild/bbs-3.15-bioc/R/lib…
##  4 SAMEA103885392 diku_D diku_1 IFNg_SL1344 /home/biocbuild/bbs-3.15-bioc/R/lib…
##  5 SAMEA103885182 eiwy_A eiwy_1 naive       /home/biocbuild/bbs-3.15-bioc/R/lib…
##  6 SAMEA103885136 eiwy_B eiwy_1 IFNg        /home/biocbuild/bbs-3.15-bioc/R/lib…
##  7 SAMEA103885413 eiwy_C eiwy_1 SL1344      /home/biocbuild/bbs-3.15-bioc/R/lib…
##  8 SAMEA103884967 eiwy_D eiwy_1 IFNg_SL1344 /home/biocbuild/bbs-3.15-bioc/R/lib…
##  9 SAMEA103885368 fikt_A fikt_3 naive       /home/biocbuild/bbs-3.15-bioc/R/lib…
## 10 SAMEA103885218 fikt_B fikt_3 IFNg        /home/biocbuild/bbs-3.15-bioc/R/lib…
## # … with 14 more rows

After we have read the coldata.csv file, we select relevant columns from this table, create a new column called files, and transform the existing line and condition columns into factors. In the case of condition, we specify the “naive” cell line as the reference level. The files column points to the quantifications for each observation – these files have been gzipped, but would typically not have the ‘gz’ ending if used from Salmon directly. One other thing to note is the use of the pipe operator,%>%, which can be read as “then”, i.e. first read the data, then select columns, then mutate them.

Now we have a table summarizing the experimental design and the locations of the quantifications. The following lines of code do a lot of work for the analyst: importing the RNA-seq quantification (dropping inferential replicates in this case), locating the relevant reference transcriptome, attaching the transcript ranges to the data, and fetching genome information. Inferential replicates are especially useful for performing transcript-level analysis, but here we will use a point estimate for the per-gene counts and perform gene-level analysis.

The result is a SummarizedExperiment object.

## importing quantifications
## reading in files with read_tsv
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 
## found matching linked transcriptome:
## [ GENCODE - Homo sapiens - release 29 ]
## useHub=TRUE: checking for TxDb via 'AnnotationHub'
## snapshotDate(): 2022-04-21
## found matching TxDb via 'AnnotationHub'
## loading from cache
## Loading required package: GenomicFeatures
## Loading required package: AnnotationDbi
## 
## Attaching package: 'AnnotationDbi'
## 
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## generating transcript ranges
## class: RangedSummarizedExperiment 
## dim: 205870 24 
## metadata(6): tximetaInfo quantInfo ... txomeInfo txdbInfo
## assays(3): counts abundance length
## rownames(205870): ENST00000456328.2 ENST00000450305.2 ...
##   ENST00000387460.2 ENST00000387461.2
## rowData names(3): tx_id gene_id tx_name
## colnames(24): SAMEA103885102 SAMEA103885347 ... SAMEA103885308
##   SAMEA103884949
## colData names(4): names id line condition

On a machine with a working internet connection, the above command works without any extra steps, as the tximeta function obtains any necessary metadata via FTP, unless it is already cached locally. The tximeta package can also be used without an internet connection, in this case the linked transcriptome can be created directly from a Salmon index and gtf.

Because tximeta knows the correct reference transcriptome, we can ask tximeta to summarize the transcript-level data to the gene level using the methods of Soneson, Love, and Robinson (2015).

## loading existing TxDb created: 2022-04-29 12:03:19
## obtaining transcript-to-gene mapping from database
## generating gene ranges
## summarizing abundance
## summarizing counts
## summarizing length

One final note is that the start of positive strand genes and the end of negative strand genes is now dictated by the genomic extent of the isoforms of the gene (so the start and end of the reduced GRanges). Another alternative would be to either operate on transcript abundance, and perform differential analysis on transcript (and so avoid defining the TSS of a set of isoforms), or to use gene-level summarized expression but to pick the most representative TSS based on isoform expression.

2.2 Importing ATAC-seq data as a SummarizedExperiment object

The SummarizedExperiment object containing ATAC-seq peaks can be created from the following tab-delimited files from Alasoo and Gaffney (2017):

To begin, we read in the sample metadata, following similar steps to those we used to generate the coldata table for the RNA-seq experiment:

The ATAC-seq counts have already been normalized with cqn (Hansen, Irizarry, and Wu 2012) and log2 transformed. Loading the cqn-normalized matrix of log2 transformed read counts takes ~30 seconds and loads an object of ~370 Mb. We set the column names so that the first column contains the rownames of the matrix, and the remaining columns are the sample identities from the atac_coldata object.

We read in the peak metadata (locations in the genome), and convert it to a GRanges object. The as_granges() function automatically converts the data.frame into a GRanges object. From that result, we extract the peak_id column and set the genome information to the build “GRCh38”. We know this from the Zenodo entry.

Finally, we construct a SummarizedExperiment object. We place the matrix into the assays slot as a named list, the annotated peaks into the row-wise ranges slot, and the sample metadata into the column-wise data slot:

3 Model assays

3.1 RNA-seq differential gene expression analysis

We can easily run a differential expression analysis with DESeq2 using the following code chunks (Love, Huber, and Anders 2014). The design formula indicates that we want to control for the donor baselines (line) and test for differences in gene expression on the condition. For a more comprehensive discussion of DE workflows in Bioconductor see Love et al. (2016) and Law et al. (2018).

## using counts and average transcript lengths from tximeta

The model is fit with the following line of code:

## estimating size factors
## using 'avgTxLength' from assays(dds), correcting for library size
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

Below we set the contrast on the condition variable, indicating we are estimating the \(\log_2\) fold change (LFC) of IFNg stimulated cell lines against naive cell lines. We are interested in LFC greater than 1 at a nominal false discovery rate (FDR) of 1%.

To see the results of the expression analysis, we can generate a summary table and an MA plot:

## 
## out of 17806 with nonzero total read count
## adjusted p-value < 0.01
## LFC > 1.00 (up)    : 502, 2.8%
## LFC < -1.00 (down) : 247, 1.4%
## outliers [1]       : 0, 0%
## low counts [2]     : 0, 0%
## (mean count < 3)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
Visualization of DESeq2 results as an “MA plot”. Genes that have an adjusted p-value below 0.01 are colored red.

Figure 3.1: Visualization of DESeq2 results as an “MA plot”. Genes that have an adjusted p-value below 0.01 are colored red.

We now output the results as a GRanges object, and due to the conventions of plyranges, we construct a new column called gene_id from the row names of the results. Each row now contains the genomic region (seqnames, start, end, strand) along with corresponding metadata columns (the gene_id and the results of the test). Note that tximeta has correctly identified the reference genome as “hg38”, and this has also been added to the GRanges along the results columns. This kind of book-keeping is vital once overlap operations are performed to ensure that plyranges is not comparing across incompatible genomes.

## GRanges object with 17806 ranges and 7 metadata columns:
##           seqnames              ranges strand |            gene_id   baseMean
##              <Rle>           <IRanges>  <Rle> |        <character>  <numeric>
##       [1]     chrX 100627109-100639991      - | ENSG00000000003.14    171.571
##       [2]    chr20   50934867-50958555      - | ENSG00000000419.12    967.751
##       [3]     chr1 169849631-169894267      - | ENSG00000000457.13    682.433
##       [4]     chr1 169662007-169854080      + | ENSG00000000460.16    262.963
##       [5]     chr1   27612064-27635277      - | ENSG00000000938.12   2660.102
##       ...      ...                 ...    ... .                ...        ...
##   [17802]    chr10   84167228-84172093      - |  ENSG00000285972.1   10.04746
##   [17803]     chr6   63572012-63583587      + |  ENSG00000285976.1 4586.34617
##   [17804]    chr16   57177349-57181390      + |  ENSG00000285979.1   14.29653
##   [17805]     chr8 103398658-103501895      - |  ENSG00000285982.1   27.76296
##   [17806]    chr10   12563151-12567351      + |  ENSG00000285994.1    6.60409
##           log2FoldChange     lfcSE      stat    pvalue      padj
##                <numeric> <numeric> <numeric> <numeric> <numeric>
##       [1]     -0.2822450 0.3005710   0.00000 1.0000000  1.000000
##       [2]      0.0391223 0.0859708   0.00000 1.0000000  1.000000
##       [3]      1.2846179 0.1969067   1.44545 0.1483329  1.000000
##       [4]     -1.4718762 0.2186916  -2.15772 0.0309493  0.409728
##       [5]      0.6754781 0.2360530   0.00000 1.0000000  1.000000
##       ...            ...       ...       ...       ...       ...
##   [17802]      0.5484518  0.444319         0         1         1
##   [17803]     -0.0339296  0.188005         0         1         1
##   [17804]      0.3123477  0.522700         0         1         1
##   [17805]      0.9945187  1.582373         0         1         1
##   [17806]      0.2539975  0.595751         0         1         1
##   -------
##   seqinfo: 25 sequences (1 circular) from hg38 genome

From this, we can restrict the results to those that meet our FDR threshold and select (and rename) the metadata columns we’re interested in:

We now wish to extract genes for which there is evidence that the LFC is not large. We perform this test by specifying an LFC threshold and an alternative hypothesis (altHypothesis) that the LFC is less than the threshold in absolute value. To visualize the result of this test, you can run results without format="GRanges", and pass this object to plotMA as before. We label these genes as other_genes and later as “non-DE genes”, for comparison with our de_genes set.

3.2 ATAC-seq peak differential abundance analysis

The following section describes the process we have used for generating a GRanges object of differential peaks from the ATAC-seq data in Alasoo et al. (2018).

The code chunks for the remainder of this section are not run.

For assessing differential accessibility, we run limma (Smyth 2004), and generate the a summary of LFCs and adjusted p-values for the peaks:

We now take the rowRanges of the SummarizedExperiment and attach the LFCs and adjusted p-values from limma, so that we can consider the overlap with differential expression. Note that we set the genome build to “hg38” and restyle the chromosome information to use the “UCSC” style (e.g. “chr1”, “chr2”, etc.). Again, we know the genome build from the Zenodo entry for the ATAC-seq data.

The final GRanges object containing the DA peaks is included in the workflow package and can be loaded as follows:

## GRanges object with 296220 ranges and 3 metadata columns:
##            seqnames              ranges strand |          peak_id da_log2FC
##               <Rle>           <IRanges>  <Rle> |      <character> <numeric>
##        [1]     chr1          9979-10668      * |      ATAC_peak_1  0.266185
##        [2]     chr1         10939-11473      * |      ATAC_peak_2  0.322177
##        [3]     chr1         15505-15729      * |      ATAC_peak_3 -0.574160
##        [4]     chr1         21148-21481      * |      ATAC_peak_4 -1.147066
##        [5]     chr1         21864-22067      * |      ATAC_peak_5 -0.896143
##        ...      ...                 ...    ... .              ...       ...
##   [296216]     chrX 155896572-155896835      * | ATAC_peak_296216 -0.834629
##   [296217]     chrX 155958507-155958646      * | ATAC_peak_296217 -0.147537
##   [296218]     chrX 156016760-156016975      * | ATAC_peak_296218 -0.609732
##   [296219]     chrX 156028551-156029422      * | ATAC_peak_296219 -0.347678
##   [296220]     chrX 156030135-156030785      * | ATAC_peak_296220  0.492442
##                da_padj
##              <numeric>
##        [1] 9.10673e-05
##        [2] 2.03435e-05
##        [3] 3.41708e-08
##        [4] 8.22299e-26
##        [5] 4.79453e-11
##        ...         ...
##   [296216] 1.33546e-11
##   [296217] 3.13015e-01
##   [296218] 3.62339e-09
##   [296219] 6.94823e-06
##   [296220] 7.07664e-13
##   -------
##   seqinfo: 23 sequences from hg38 genome; no seqlengths

4 Integrate ranges

4.1 Finding overlaps with plyranges

We have already used plyranges a number of times above, to filter, mutate, and select on GRanges objects, as well as ensuring the correct genome annotation and style has been used. The plyranges package provides a grammar for performing transformations of genomic data (Lee, Cook, and Lawrence 2019). Computations resulting from compositions of plyranges “verbs” are performed using underlying, highly optimized range operations in the GenomicRanges package (Lawrence et al. 2013).

For the overlap analysis, we filter the annotated peaks to have a nominal FDR bound of 1%.

We now have GRanges objects that contain DE genes, genes without strong signal of DE, and DA peaks. We are ready to answer the question: is there an enrichment of DA ATAC-seq peaks in the vicinity of DE genes compared to genes without sufficient DE signal?

4.2 Down sampling non-differentially expressed genes

As plyranges is built on top of dplyr, it implements methods for many of its verbs for GRanges objects. Here we can use slice to randomly sample the rows of the other_genes. The sample.int function will generate random samples of size equal to the number of DE-genes from the number of rows in other_genes:

## GRanges object with 749 ranges and 3 metadata columns:
##         seqnames              ranges strand |            gene_id   de_log2FC
##            <Rle>           <IRanges>  <Rle> |        <character>   <numeric>
##     [1]     chr7     5898725-5926550      + | ENSG00000122674.11   0.2031642
##     [2]    chr12   22625075-22690665      + | ENSG00000139163.15   0.2575908
##     [3]     chr1   43974487-43978295      + | ENSG00000117410.13   0.1148944
##     [4]    chr19   55635016-55643469      + |  ENSG00000213015.8  -0.0892649
##     [5]    chr16   21402237-21448567      - | ENSG00000169246.16   0.1878874
##     ...      ...                 ...    ... .                ...         ...
##   [745]     chrX   15319451-15335580      - | ENSG00000165195.15  0.22848693
##   [746]     chr5 134194334-134226142      - |  ENSG00000113575.9  0.06052272
##   [747]    chr12      990509-1495933      + | ENSG00000082805.19  0.00646525
##   [748]    chr11   61130256-61161617      - | ENSG00000167987.10 -0.09338839
##   [749]    chr17   17810399-17837002      - | ENSG00000072310.16  0.30044350
##             de_padj
##           <numeric>
##     [1] 7.06653e-07
##     [2] 4.69348e-06
##     [3] 4.53819e-13
##     [4] 2.61090e-08
##     [5] 6.06624e-07
##     ...         ...
##   [745] 2.09461e-07
##   [746] 8.24597e-24
##   [747] 8.07219e-09
##   [748] 1.09401e-12
##   [749] 2.45385e-04
##   -------
##   seqinfo: 25 sequences (1 circular) from hg38 genome

We can repeat this many times to create many samples via replicate. By replicating the sub-sampling multiple times, we minimize the variance on the enrichment statistics induced by the sampling process.

This creates a list of GRanges objects as a list, and we can bind these together using the bind_ranges function. This function creates a new column called “resample” on the result that identifies each of the input GRanges objects:

Similarly, we can then combine the boot_genes GRanges, with the DE GRanges object. As the resample column was not present on the DE GRanges object, this is given a missing value which we recode to a 0 using mutate()

## GRanges object with 8239 ranges and 5 metadata columns:
##          seqnames              ranges strand |            gene_id de_log2FC
##             <Rle>           <IRanges>  <Rle> |        <character> <numeric>
##      [1]     chr1 196651878-196747504      + | ENSG00000000971.15   4.98711
##      [2]     chr6   46129993-46146699      + |  ENSG00000001561.6   1.92722
##      [3]     chr4   17577192-17607972      + | ENSG00000002549.12   2.93373
##      [4]     chr7 150800403-150805120      + |  ENSG00000002933.8   3.16722
##      [5]     chr4   15778275-15853230      + | ENSG00000004468.12   5.40894
##      ...      ...                 ...    ... .                ...       ...
##   [8235]    chr17   43527844-43579620      - | ENSG00000175832.12 -0.240918
##   [8236]    chr17   18260534-18266552      + | ENSG00000177427.12 -0.166059
##   [8237]    chr20   63895182-63936031      + | ENSG00000101152.10  0.250539
##   [8238]     chr1   39081316-39487177      + | ENSG00000127603.25 -0.385054
##   [8239]     chr8   41577187-41625001      + | ENSG00000158669.11  0.155922
##              de_padj  resample   origin
##            <numeric> <integer> <factor>
##      [1] 1.37057e-13         0       de
##      [2] 3.17478e-05         0       de
##      [3] 2.01310e-11         0       de
##      [4] 1.07360e-08         0       de
##      [5] 4.82905e-18         0       de
##      ...         ...       ...      ...
##   [8235] 9.91611e-03        10   not_de
##   [8236] 9.12051e-05        10   not_de
##   [8237] 1.74085e-09        10   not_de
##   [8238] 2.65539e-03        10   not_de
##   [8239] 2.96375e-17        10   not_de
##   -------
##   seqinfo: 25 sequences (1 circular) from hg38 genome

4.3 Expanding genomic coordinates around the transcription start site

Now we would like to modify our gene ranges so they contain the 10 kilobases on either side of their transcription start site (TSS). There are many ways one could do this, but we prefer an approach via the anchoring methods in plyranges. Because there is a mutual dependence between the start, end, width, and strand of a GRanges object, we define anchors to fix one of start and end, while modifying the width. As an example, to extract just the TSS, we can anchor by the 5’ end of the range and modify the width of the range to equal 1.

Anchoring by the 5’ end of a range will fix the end of negatively stranded ranges, and fix the start of positively stranded ranges.

We can then repeat the same pattern but this time using anchor_center() to tell plyranges that we are making the TSS the midpoint of a range that has total width of 20kb, or 10kb both upstream and downstream of the TSS.

4.4 Use overlap joins to find relative enrichment

We are now ready to compute overlaps between RNA-seq genes (our DE set and bootstrap sets) and the ATAC-seq peaks. In plyranges, overlaps are defined as joins between two GRanges objects: a left and a right GRanges object. In an overlap join, a match is any range on the left GRanges that is overlapped by the right GRanges. One powerful aspect of the overlap joins is that the result maintains all (metadata) columns from each of the left and right ranges which makes downstream summaries easy to compute.

To combine the DE genes with the DA peaks, we perform a left overlap join. This returns to us the all_genes ranges (potentially with duplication), but with the metadata columns from those overlapping DA peaks. For any gene that has no overlaps, the DA peak columns will have NA’s.

## GRanges object with 27766 ranges and 8 metadata columns:
##           seqnames              ranges strand |            gene_id de_log2FC
##              <Rle>           <IRanges>  <Rle> |        <character> <numeric>
##       [1]     chr1 196641878-196661877      + | ENSG00000000971.15   4.98711
##       [2]     chr6   46119993-46139992      + |  ENSG00000001561.6   1.92722
##       [3]     chr4   17567192-17587191      + | ENSG00000002549.12   2.93373
##       [4]     chr4   17567192-17587191      + | ENSG00000002549.12   2.93373
##       [5]     chr4   17567192-17587191      + | ENSG00000002549.12   2.93373
##       ...      ...                 ...    ... .                ...       ...
##   [27762]     chr1   39071316-39091315      + | ENSG00000127603.25 -0.385054
##   [27763]     chr1   39071316-39091315      + | ENSG00000127603.25 -0.385054
##   [27764]     chr8   41567187-41587186      + | ENSG00000158669.11  0.155922
##   [27765]     chr8   41567187-41587186      + | ENSG00000158669.11  0.155922
##   [27766]     chr8   41567187-41587186      + | ENSG00000158669.11  0.155922
##               de_padj  resample   origin          peak_id da_log2FC     da_padj
##             <numeric> <integer> <factor>      <character> <numeric>   <numeric>
##       [1] 1.37057e-13         0       de  ATAC_peak_21236 -0.546582 1.15274e-04
##       [2] 3.17478e-05         0       de ATAC_peak_231183  1.453297 9.73225e-17
##       [3] 2.01310e-11         0       de ATAC_peak_193578  0.222371 3.00939e-11
##       [4] 2.01310e-11         0       de ATAC_peak_193579 -0.281615 7.99889e-05
##       [5] 2.01310e-11         0       de ATAC_peak_193580  0.673705 7.60043e-15
##       ...         ...       ...      ...              ...       ...         ...
##   [27762] 2.65539e-03        10   not_de   ATAC_peak_5357 -1.058236 3.69052e-16
##   [27763] 2.65539e-03        10   not_de   ATAC_peak_5358 -1.314112 6.44280e-26
##   [27764] 2.96375e-17        10   not_de ATAC_peak_263396 -0.904080 8.19577e-13
##   [27765] 2.96375e-17        10   not_de ATAC_peak_263397  0.364738 2.08835e-08
##   [27766] 2.96375e-17        10   not_de ATAC_peak_263399  0.317387 1.20088e-08
##   -------
##   seqinfo: 25 sequences (1 circular) from hg38 genome

Now we can ask, how many DA peaks are near DE genes relative to “other” non-DE genes? A gene may appear more than once in genes_olap_peaks, because multiple peaks may overlap a single gene, or because we have re-sampled the same gene more than once, or a combination of these two cases.

For each gene (that is the combination of chromosome, the start, end, and strand), and the “origin” (DE vs not-DE) we can compute the distinct number of peaks for each gene and the maximum peak based on LFC. This is achieved via reduce_ranges_directed, which allows an aggregation to result in a GRanges object via merging neighboring genomic regions. The use of the directed suffix indicates we’re maintaining strand information. In this case, we are simply merging ranges (genes) via the groups we mentioned above. We also have to account for the number of resamples we have performed when counting if there are any peaks, to ensure we do not double count the same peak:

We can then filter genes if they have any peaks and compare the peak fold changes between non-DE and DE genes using a boxplot:

A boxplot of maximum LFCs for DA peaks for DE genes compared to non-DE genes where genes have at least one DA peak.

Figure 4.1: A boxplot of maximum LFCs for DA peaks for DE genes compared to non-DE genes where genes have at least one DA peak.

In general, the DE genes have larger maximum DA fold changes relative to the non-DE genes.

Next we examine how thresholds on the DA LFC modify the enrichment we observe of DA peaks near DE or non-DE genes. First, we want to know how the number of peaks within DE genes and non-DE genes change as we change threshold values on the peak LFC. As an example, we could compute this by arbitrarily chosen LFC thresholds of 1 or 2 as follows:

## DataFrame with 2 rows and 4 columns
##     origin peak_count lfc1_peak_count lfc2_peak_count
##   <factor>  <numeric>       <numeric>       <numeric>
## 1   not_de     2391.8           369.5            32.5
## 2   de         3416.0          1097.0           234.0

Here we see that DE genes tend to have more DA peaks near them, and that the number of DA peaks decreases as we increase the DA LFC threshold (as expected). We now show how to compute the ratio of peak counts from DE compared to non-DE genes, so we can see how this ratio changes for various DA LFC thresholds.

For all variables except for the origin column we divide the first row’s values by the second row, which will be the enrichment of peaks in DE genes compared to other genes. This requires us to reshape the summary table from long form back to wide form using the tidyr package. First we pivot the results of the peak_count columns into name-value pairs, then pivot again to place values into the origin column. Then we create a new column with the relative enrichment:

## # A tibble: 3 × 4
##   name            not_de    de enrichment
##   <chr>            <dbl> <dbl>      <dbl>
## 1 peak_count      2392.   3416       1.43
## 2 lfc1_peak_count  370.   1097       2.97
## 3 lfc2_peak_count   32.5   234       7.2

The above table shows that relative enrichment increases for a larger LFC threshold.

Due to the one-to-many mappings of genes to peaks, it is unknown if we have the same number of DE genes participating or less, as we increase the threshold on the DA LFC. We can examine the number of genes with overlapping DA peaks at various thresholds by grouping and aggregating twice. First, the number of peaks that meet the thresholds are computed within each gene, origin, and resample group. Second, within the origin column, we compute the total number of peaks that meet the DA LFC threshold and the number of genes that have more than zero peaks (again averaging over the number of resamples).

## DataFrame with 2 rows and 5 columns
##     origin lfc1_gene_count lfc1_peak_count lfc2_gene_count lfc2_peak_count
##   <factor>       <numeric>       <numeric>       <numeric>       <numeric>
## 1   not_de           271.2           369.5            30.3            32.5
## 2   de               515.0          1097.0           151.0           234.0

To do this for many thresholds is cumbersome and would create a lot of duplicate code. Instead we create a single function called count_above_threshold that accepts a variable and a vector of thresholds, and computes the sum of the absolute value of the variable for each element in the thresholds vector.

The above function will compute the counts for any arbitrary threshold, so we can apply it over possible LFC thresholds of interest. We choose a grid of one hundred thresholds based on the range of absolute LFC values in the da_peaks GRanges object:

The peak counts for each threshold are computed as a new list-column called value. First, the GRanges object has been grouped by the gene, origin, and the number of resamples columns. Then we aggregate over those columns, so each row will contain the peak counts for all of the thresholds for a gene, origin, and resample. We also maintain another list-column that contains the threshold values.

## GRanges object with 8239 ranges and 5 metadata columns:
##          seqnames              ranges strand |            gene_id   origin
##             <Rle>           <IRanges>  <Rle> |        <character> <factor>
##      [1]     chr1 196641878-196661877      + | ENSG00000000971.15       de
##      [2]     chr6   46119993-46139992      + |  ENSG00000001561.6       de
##      [3]     chr4   17567192-17587191      + | ENSG00000002549.12       de
##      [4]     chr7 150790403-150810402      + |  ENSG00000002933.8       de
##      [5]     chr4   15768275-15788274      + | ENSG00000004468.12       de
##      ...      ...                 ...    ... .                ...      ...
##   [8235]    chr17   43569620-43589619      - | ENSG00000175832.12   not_de
##   [8236]    chr17   18250534-18270533      + | ENSG00000177427.12   not_de
##   [8237]    chr20   63885182-63905181      + | ENSG00000101152.10   not_de
##   [8238]     chr1   39071316-39091315      + | ENSG00000127603.25   not_de
##   [8239]     chr8   41567187-41587186      + | ENSG00000158669.11   not_de
##           resample         value                         threshold
##          <integer> <IntegerList>                     <NumericList>
##      [1]         0     1,1,1,... 0.0658243,0.1184840,0.1711436,...
##      [2]         0     1,1,1,... 0.0658243,0.1184840,0.1711436,...
##      [3]         0     6,6,6,... 0.0658243,0.1184840,0.1711436,...
##      [4]         0     4,4,4,... 0.0658243,0.1184840,0.1711436,...
##      [5]         0  11,11,11,... 0.0658243,0.1184840,0.1711436,...
##      ...       ...           ...                               ...
##   [8235]        10     1,1,1,... 0.0658243,0.1184840,0.1711436,...
##   [8236]        10     3,3,2,... 0.0658243,0.1184840,0.1711436,...
##   [8237]        10     5,5,5,... 0.0658243,0.1184840,0.1711436,...
##   [8238]        10     3,3,3,... 0.0658243,0.1184840,0.1711436,...
##   [8239]        10     3,3,3,... 0.0658243,0.1184840,0.1711436,...
##   -------
##   seqinfo: 25 sequences (1 circular) from hg38 genome

Now we can expand these list-columns into a long GRanges object using the expand_ranges() function. This function will unlist the value and threshold columns and lengthen the resulting GRanges object. To compute the peak and gene counts for each threshold, we apply the same summarization as before:

## DataFrame with 200 rows and 4 columns
##       origin threshold gene_count peak_count
##     <factor> <numeric>  <numeric>  <numeric>
## 1     not_de 0.0658243      708.0     2391.4
## 2     not_de 0.1184840      698.8     2320.6
## 3     not_de 0.1711436      686.2     2178.6
## 4     not_de 0.2238033      672.4     1989.4
## 5     not_de 0.2764629      650.4     1785.8
## ...      ...       ...        ...        ...
## 196       de   5.06849          2          2
## 197       de   5.12115          0          0
## 198       de   5.17381          0          0
## 199       de   5.22647          0          0
## 200       de   5.27913          0          0

Again we can compute the relative enrichment in LFCs in the same manner as before, by pivoting the results to long form then back to wide form to compute the enrichment. We visualize the peak enrichment changes of DE genes relative to other genes as a line chart:

## Warning: Removed 4 row(s) containing missing values (geom_path).
A line chart displaying how relative enrichment of DA peaks change between DE genes compared to non-DE genes as the absolute DA LFC threshold increases.

Figure 4.2: A line chart displaying how relative enrichment of DA peaks change between DE genes compared to non-DE genes as the absolute DA LFC threshold increases.

We computed the sum of DA peaks near the DE genes, for increasing LFC thresholds on the accessibility change. As we increased the threshold, the number of total peaks went down (likewise the mean number of DA peaks per gene). It is also likely the number of DE genes with a DA peak nearby with such a large change went down. We can investigate this with a plot that summarizes many of the aspects underlying the enrichment plot above.

A line chart displaying how gene and peak counts change as the absolute DA LFC threshold increases. Lines are colored according to whether they represent a gene that is DE or not. Note the x-axis is on a \(\log_{10}\) scale.

Figure 4.3: A line chart displaying how gene and peak counts change as the absolute DA LFC threshold increases. Lines are colored according to whether they represent a gene that is DE or not. Note the x-axis is on a \(\log_{10}\) scale.

5 Discussion

We have shown that by using plyranges and tximeta (with support of Bioconductor and tidyverse ecosystems) we can fluently iterate through the biological data science workflow: from import, through to modeling, and data integration.

There are several further steps that would be interesting to perform in this analysis; for example, we could modify window size around the TSS to see how it affects enrichment, and vary the FDR cut-offs for both the DE gene and DA peak sets. We could also have computed variance in addition to the mean of the bootstrap set, and so drawn an interval around the enrichment line.

Finally, our workflow illustrates the benefits of using appropriate data abstractions provided by Bioconductor such as the SummarizedExperiment and GRanges. These abstractions provide users with a mental model of their experimental data and are the building blocks for constructing the modular and iterative analyses we have shown here. Consequently, we have been able to interoperate many decoupled R packages (from both Bioconductor and the tidyverse) to construct a seamless end-to-end workflow that is far too specialized for a single monolithic tool.

6 Software Availability

The workflow materials, including this article can be fully reproduced following the instructions found at the Github repository sa-lee/fluentGenomics. Moreover, the development version of the workflow and all downstream dependencies can be installed using the BiocManager package by running:

# development version from Github
BiocManager::install("sa-lee/fluentGenomics")
# version available from Bioconductor
BiocManager::install("fluentGenomics")

This article and the analyses were performed with R (R Core Team 2019) using the rmarkdown (Allaire et al. 2019), and knitr (Xie 2019, 2015) packages.

6.1 Session Info

## R version 4.2.0 RC (2022-04-19 r82224)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.15-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.15-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB              LC_COLLATE=C              
##  [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       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] ggplot2_3.3.5               plyranges_1.16.0           
##  [3] DESeq2_1.36.0               GenomicFeatures_1.48.0     
##  [5] AnnotationDbi_1.58.0        SummarizedExperiment_1.26.0
##  [7] Biobase_2.56.0              GenomicRanges_1.48.0       
##  [9] GenomeInfoDb_1.32.1         IRanges_2.30.0             
## [11] S4Vectors_0.34.0            BiocGenerics_0.42.0        
## [13] MatrixGenerics_1.8.0        matrixStats_0.62.0         
## [15] readr_2.1.2                 dplyr_1.0.9                
## [17] tximeta_1.14.0              fluentGenomics_1.8.0       
## 
## loaded via a namespace (and not attached):
##   [1] colorspace_2.0-3              rjson_0.2.21                 
##   [3] ellipsis_0.3.2                XVector_0.36.0               
##   [5] farver_2.1.0                  bit64_4.0.5                  
##   [7] interactiveDisplayBase_1.34.0 fansi_1.0.3                  
##   [9] xml2_1.3.3                    splines_4.2.0                
##  [11] tximport_1.24.0               cachem_1.0.6                 
##  [13] geneplotter_1.74.0            knitr_1.39                   
##  [15] jsonlite_1.8.0                Rsamtools_2.12.0             
##  [17] annotate_1.74.0               dbplyr_2.1.1                 
##  [19] png_0.1-7                     shiny_1.7.1                  
##  [21] BiocManager_1.30.17           compiler_4.2.0               
##  [23] httr_1.4.2                    assertthat_0.2.1             
##  [25] Matrix_1.4-1                  fastmap_1.1.0                
##  [27] lazyeval_0.2.2                cli_3.3.0                    
##  [29] later_1.3.0                   htmltools_0.5.2              
##  [31] prettyunits_1.1.1             tools_4.2.0                  
##  [33] gtable_0.3.0                  glue_1.6.2                   
##  [35] GenomeInfoDbData_1.2.8        rappdirs_0.3.3               
##  [37] Rcpp_1.0.8.3                  jquerylib_0.1.4              
##  [39] vctrs_0.4.1                   Biostrings_2.64.0            
##  [41] rtracklayer_1.56.0            xfun_0.30                    
##  [43] stringr_1.4.0                 mime_0.12                    
##  [45] lifecycle_1.0.1               restfulr_0.0.13              
##  [47] ensembldb_2.20.0              XML_3.99-0.9                 
##  [49] AnnotationHub_3.4.0           scales_1.2.0                 
##  [51] zlibbioc_1.42.0               vroom_1.5.7                  
##  [53] hms_1.1.1                     promises_1.2.0.1             
##  [55] ProtGenerics_1.28.0           parallel_4.2.0               
##  [57] RColorBrewer_1.1-3            AnnotationFilter_1.20.0      
##  [59] yaml_2.3.5                    curl_4.3.2                   
##  [61] memoise_2.0.1                 sass_0.4.1                   
##  [63] biomaRt_2.52.0                stringi_1.7.6                
##  [65] RSQLite_2.2.12                BiocVersion_3.15.2           
##  [67] highr_0.9                     genefilter_1.78.0            
##  [69] BiocIO_1.6.0                  filelock_1.0.2               
##  [71] BiocParallel_1.30.0           rlang_1.0.2                  
##  [73] pkgconfig_2.0.3               bitops_1.0-7                 
##  [75] evaluate_0.15                 lattice_0.20-45              
##  [77] purrr_0.3.4                   labeling_0.4.2               
##  [79] GenomicAlignments_1.32.0      bit_4.0.4                    
##  [81] tidyselect_1.1.2              magrittr_2.0.3               
##  [83] bookdown_0.26                 R6_2.5.1                     
##  [85] generics_0.1.2                DelayedArray_0.22.0          
##  [87] DBI_1.1.2                     pillar_1.7.0                 
##  [89] withr_2.5.0                   survival_3.3-1               
##  [91] KEGGREST_1.36.0               RCurl_1.98-1.6               
##  [93] tibble_3.1.6                  crayon_1.5.1                 
##  [95] utf8_1.2.2                    BiocFileCache_2.4.0          
##  [97] tzdb_0.3.0                    rmarkdown_2.14               
##  [99] progress_1.2.2                locfit_1.5-9.5               
## [101] grid_4.2.0                    blob_1.2.3                   
## [103] digest_0.6.29                 xtable_1.8-4                 
## [105] tidyr_1.2.0                   httpuv_1.6.5                 
## [107] munsell_0.5.0                 bslib_0.3.1

6.2 Author Contributions

All authors contributed to the writing and development of the workflow.

6.3 Competing interests

The authors declare that they have no competing interests.

6.4 Funding

SL is supported by an Australian Government Research Training Program (RTP) scholarship with a top up scholarship from CSL Limited.

MIL’s contribution is supported by NIH grant R01 HG009937.

I confirm that the funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

6.5 Acknowledgements

The authors would like to thank all participants of the Bioconductor 2019 and BiocAsia 2019 conferences who attended and provided feedback on early versions of this workflow paper.

References

Alasoo, Kaur, and Daniel Gaffney. 2017. “Processed read counts from macrophage RNA-seq and ATAC-seq experiments.” Zenodo. https://doi.org/10.5281/zenodo.1188300.

Alasoo, K, J Rodrigues, S Mukhopadhyay, AJ Knights, AL Mann, K Kundu, HIPSCI-Consortium, C Hale, Dougan G, and DJ Gaffney. 2018. “Shared genetic effects on chromatin and gene expression indicate a role for enhancer priming in immune response.” Nature Genetics 50: 424–31. https://doi.org/10.1038/s41588-018-0046-7.

Allaire, JJ, Yihui Xie, Jonathan McPherson, Javier Luraschi, Kevin Ushey, Aron Atkins, Hadley Wickham, Joe Cheng, Winston Chang, and Richard Iannone. 2019. Rmarkdown: Dynamic Documents for R. https://github.com/rstudio/rmarkdown.

Frankish, A, GENCODE-consoritum, and P Flicek. 2018. “GENCODE reference annotation for the human and mouse genomes.” Nucleic Acids Research. https://doi.org/10.1093/nar/gky955.

Hansen, Kasper D., Rafael A. Irizarry, and Zhijin Wu. 2012. “Removing technical variability in RNA-seq data using condition quantile normalization.” Biostatistics 13 (2): 204–16. https://doi.org/10.1093/biostatistics/kxr054.

Huber, Wolfgang, Vincent J Carey, Robert Gentleman, Simon Anders, Marc Carlson, Benilton S Carvalho, Hector Corrada Bravo, et al. 2015. “Orchestrating High-Throughput Genomic Analysis with Bioconductor.” Nature Methods 12 (2): 115–21. https://doi.org/10.1038/nmeth.3252.

Köster, Johannes, and Sven Rahmann. 2012. “Snakemake—a scalable bioinformatics workflow engine.” Bioinformatics 28 (19): 2520–2. https://doi.org/10.1093/bioinformatics/bts480.

Law, Charity W, Monther Alhamdoosh, Shian Su, Xueyi Dong, Luyi Tian, Gordon K Smyth, and Matthew E Ritchie. 2018. “RNA-seq Analysis Is Easy as 1-2-3 with Limma, Glimma and edgeR.” F1000 Research 5 (1408): 1408. https://doi.org/10.12688/f1000research.9005.3.

Lawrence, Michael, Wolfgang Huber, Hervé Pagès, Patrick Aboyoun, Marc Carlson, Robert Gentleman, Martin Morgan, and Vincent Carey. 2013. “Software for Computing and Annotating Genomic Ranges.” PLoS Comput. Biol. 9. https://doi.org/10.1371/journal.pcbi.1003118.

Lee, Stuart, Dianne Cook, and Michael Lawrence. 2019. “Plyranges: A Grammar of Genomic Data Transformation.” Genome Biology 20 (1): 4. https://doi.org/10.1186/s13059-018-1597-8.

Love, Michael I, Simon Anders, Vladislav Kim, and Wolfgang Huber. 2016. “RNA-Seq Workflow: Gene-Level Exploratory Analysis and Differential Expression.” F1000 Research 4 (1070): 1070. https://doi.org/10.12688/f1000research.7035.2.

Love, Michael I., Wolfgang Huber, and Simon Anders. 2014. “Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.” Genome Biology 15 (12): 550. https://doi.org/10.1186/s13059-014-0550-8.

Love, Michael I, Charlotte Soneson, Peter F Hickey, Lisa K Johnson, N Tessa Pierce, Lori Shepherd, Martin Morgan, and Rob Patro. 2019. “Tximeta: Reference Sequence Checksums for Provenance Identification in RNA-seq.” bioRxiv, September, 777888. https://doi.org/10.1101/777888.

Patro, R, G Duggal, MI Love, RA Irizarry, and C Kingsford. 2017. “Salmon Provides Fast and Bias-Aware Quantification of Transcript Expression.” Nature Methods 14: 417–19. https://doi.org/10.1038/nmeth.4197.

R Core Team. 2019. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Shepherd, Lori, and Martin Morgan. 2019. BiocFileCache: Manage Files Across Sessions.

Smyth, Gordon K. 2004. “Linear Models and Empirical Bayes Methods for Assessing Differential Expression in Microarray Experiments.” Statistical Applications in Genetics and Molecular Biology 3 (1).

Soneson, Charlotte, Michael I. Love, and Mark Robinson. 2015. “Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences.” F1000Research 4 (1521). https://doi.org/10.12688/f1000research.7563.1.

Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.

Xie, Yihui. 2015. Dynamic Documents with R and Knitr. 2nd ed. Boca Raton, Florida: Chapman; Hall/CRC. https://yihui.name/knitr/.

———. 2019. Knitr: A General-Purpose Package for Dynamic Report Generation in R. https://yihui.name/knitr/.