Fluent genomics with plyranges and tximeta

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.


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.
Throughout the workflow (Figure 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 et al.  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.

Experimental data
The data used in this workflow is available from two packages: the macrophage Bioconductor ExperimentData package and from the workflow package fluentGenomics (Lee & Love, 2020).
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 et al., 2019). 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 & Rahmann, 2012).
The fluentGenomics package (Lee & Love, 2020) contains functionality to download and generate a cached SummarizedExperiment object from the normalized ATAC-seq data provided by Alasoo & 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 & Morgan, 2019).

library(fluentGenomics) path_to_se <-cache_atac_se()
We can then read the cached file and assign it to an object called atac.

atac <-readRDS(path_to_se)
A precise description of how we obtained this SummarizedExperiment object can be found in Importing ATAC-seq data as a SummarizedExperiment object.

Import data as a SummarizedExperiment
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: dir <-"/path/to/quant/files" 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. dir <-system.file("extdata", package="macrophage") 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. 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.

makeLinkedTxome(
indexDir=file.path(dir, "gencode.v29_salmon_0.12.0"), source="Gencode", organism="Homo sapiens", release="29", genome="GRCh38", fasta="ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_29/ gencode.v29.transcripts.fa.gz", gtf=file.path(dir, "gencode.v29.annotation.gtf.gz"), # local version write=FALSE ) 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 et al. (2015). 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.
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 & Gaffney (2017): • The sample metadata: ATAC_sample_metadata.txt.gz (<1M) • The matrix of normalized read counts: ATAC_cqn_matrix.txt.gz (109M) • The annotated peaks: ATAC_peak_metadata.txt.gz (5.6M) 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: atac_coldata <-read_tsv("ATAC_sample_metadata.txt.gz") %>% select( sample_id, donor, condition = condition_name ) %>% mutate(condition = relevel(factor(condition), "naive")) The ATAC-seq counts have already been normalized with cqn (Hansen et al., 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.
atac_mat <-read_tsv("ATAC_cqn_matrix.txt.gz", skip = 1, col_names =c("rownames", atac_coldata[["sample_id"]])) rownames <-atac_mat[["rownames"]] atac_mat <-as.matrix(atac_mat[,-1]) rownames(atac_mat) <-rownames 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: . library(plyranges) peaks_df <-read_tsv("ATAC_peak_metadata.txt.gz", col_types = c("cidciicdc") ) peaks_gr <-peaks_df %>% as_granges(seqnames = chr) %>% select(peak_id=gene_id) %>% set_genome_info(genome = "GRCh38") 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: atac <-SummarizedExperiment(assays = list(cqndata=atac_mat), rowRanges=peaks_gr, colData=atac_coldata) Model assays RNA-seq differential gene expression analysis We can easily run a differential expression analysis with DESeq2 using the following code chunks (  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.
suppressPackageStartupMessages(library(plyranges)) de_genes <-results(dds, contrast=c("condition","IFNg","naive"), lfcThreshold=1, format="GRanges") %>% names_to_column("gene_id") de_genes ## GRanges object with 17806 ranges and 7 metadata columns: 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.
The code chunks for the remainder of this section are not run.
# set a seed for the results set.seed(2019-08-02) boot_genes <-replicate(10, slice(other_genes, sample.int(n(), size)), 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: boot_genes <-bind_ranges(boot_genes, .id = "resample") 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() all_genes <-bind_ranges( de=de_genes, not_de = boot_genes, .id="origin" ) %>% mutate( origin = factor(origin, c("not_de", "de")), resample = ifelse(is.na(resample), 0L,as.integer(resample)) ) all_genes 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.
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 20 kb, or 10 kb both upstream and downstream of the TSS.
all_genes <-all_genes %>% anchor_center() %>% mutate(width=2*1e4) 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. 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: gene_peak_max_lfc <-genes_olap_peaks %>% group_by(gene_id, origin) %>% reduce_ranges_directed( peak_count = sum(!is.na(da_padj)) / n_distinct(resample), peak_max_lfc = max(abs(da_log2FC)) ) 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 ( Figure 3 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 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 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. count_if_above_threshold <-function(var, thresholds) { lapply(thresholds, function(.) sum(abs(var) > ., na.rm = TRUE)) } 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: thresholds <-da_peaks %>% mutate(abs_lfc = abs(da_log2FC)) %>% with( seq(min(abs_lfc), max(abs_lfc), length.out = 100) ) 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. 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: 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 (Figure 4): origin_threshold_counts <-origin_peak_all_thresholds %>% as.data.frame() %>% tidyr::pivot_longer(cols = -c(origin, threshold), names_to = c("type", "var"), names_sep = "_", values_to = "count") %>% select(-var) origin_threshold_counts %>% filter(type == "peak") %>% tidyr::pivot_wider(names_from = origin, values_from = count) %>% mutate(enrichment = de / not_de) %>% ggplot(aes(x = threshold, y = enrichment)) + geom_line() + labs(x = "logFC threshold", y = "Relative Enrichment") ## Warning: Removed 4 row(s) containing missing values (geom_path).

Figure 4. 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 ( Figure 5) that summarizes many of the aspects underlying the enrichment plot above.

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.

Data availability
All data underlying the results are available as part of the article and no additional source data are required.

Author contributions
All authors contributed to the writing and development of the workflow. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R.  2. It is unclear who the audience for this paper is. As written it seems to presume a lot of technical knowledge about Bioconductor, R, and specific Bioconductor packages. It would not make a good starting point for someone who doesn't have large familiarity with all of the above. This is too bad because it has the bones of something that could be useful to people with less expertise if some additional explanatory text were provided and clearer choices made about naming variables. There is also a danger that this will be used as a cookbook without understanding.
I have marked some examples of things that could use more explanation but this is non-exhaustive. It would be worthwhile for the authors to examine the manuscript again and ask whether all of the large number of concepts, jargon, and names introduced are explained.
3. The goal of the paper appeared to be clear after reading the paper. The introduction, figure 1, and following sections, however, leave the impression that the authors tried to reproduce results from another publication which is not the case. The paper would benefit from clarifying this toward the end of the introduction. 4. Abstract: it is not clear which problem does the workflow solve and why is it important to solve it. For example, how does the import, model, and integrate structure increases result reproducibility? Is this structure new? Why would someone use it? We think that the abstract would benefit from more justification and moving the text specific to the tools used (tximeta, SummarizedExperiment) to the introduction.
5. Abstract: The sentence "Using tximeta, ..." and the following can be confusing. Why is it only using of RNA-seq in the first sentence, and both RNA-seq and ATAC-seq in the following one? 6. p3: Figure 1: The meaning of some of the objects is not immediately apparent. For example: coldata, se, dds, boot_genes. This figure does not make a great introduction for someone who is not already familiar with this workflow because so many of the details are mystifying. You should either simplify or explain more. 7. p3: Figure 1: Text in boxes is blurry.
8. p3: meaning of "fluent" should be described briefly here and not assumed 9. Introduction: This section would benefit from a more thorough description of the tools used in the workflow and explain why they are important. For example: What is Bioconductor? Why is plyranges workflow and explain why they are important. For example: What is Bioconductor? Why is plyranges relevant for this analysis? What are tidyverse design principles and why are they relevant here?
10. Introduction: Authors should also explain in the introduction why the "import, model, integrate" structure benefits users.
11. The authors should consistently use the same tense (we examine … we will …).
13. It would be useful to note that "genomic ranges" are often referred to as intervals elsewhere such as in BEDTools.
14. Please cite Tximeta at the first occurrence. 59. p13: "This creates a list of _GRanges_ objects as a list". It would be clearer to say that this creates a _list_ of 10 _GRanges_ objects. 60. What are "verbs" for GRanges 61. p13: The importance of doing the bootstrap procedure here is really unclear 62. p14: Why is there a circular sequence? I assume this is the mitochondrial chromosome? Does its presence in the background disturb the analysis at all? 63. p14: "Now we would like to modify our gene ranges so they contain…" Why? 64. p14: should be 5 prime not 5 apostrophe 65. p15: the significance of the clarification that each "gene" is just a tuple of chromosome, start, end, strand is unclear 66. p15: A lot of na.rm or is.na from here outward. Would be worth explaining why this is.  Figure 4 is non-monotonic 70. p20: Saying the x-axis is on a log10 scale is a bit confusing because the y-axis is also on a log10 scale. It is better to make log scales more obvious by having non-uniform tick marks and grid lines. Also I don't think the x-axis is on a log10 scale, log fold change here is measured with log2. (This would be more obvious if in the text you called it \log2 fold change consistently instead of the ambiguous "LFC".) 71. p21: It is unfortunate that the main point of this paper is to show that you can do something "fluently", while what that means is never discussed or justified in the paper. 72. p21: "There are several further steps that would be interesting to perform", Are these steps easy to implement with the tools you described? If so, can you explain how easier it is with your approach? 73. p21: "Finally, our workflow illustrates the benefits of using appropriate data abstractions", This is probably obvious for R power users but not for other readers. Can you point at a specific section that shows that a user will benefit from using your approach? 74. p21: The data availability statement is confusing. Clearly not all data are part of the article, they are part of external packages or articles (some specified without a particular URL or accession number, just reference to a paper). 75. p21: instructions for installing BiocManager would be helpful, as would listing the *minimum* version of R or any other relevant packages. Also should mention that you need to library(BiocManager) first 76. p21: It is unclear why you would want the development version, especially since it seems to be an older version number than what's on Bioconductor 77. p23: Items in the bibliography are not consistently referenced with specificity in a visible way. DOI should be provided for all Zenodo entries not just some. Some other documents that are not referenced via traditional publishing (journal volume and page number, or monograph from a publisher) should probably have DOI or URL specified also. Discretionary recommendations: 78. p4: `path_to_se` is a poor name for a variable. It is not orthogonal to the other variable names used here and has little explanatory power. I recommend `atac_filename7 9. p4: `dir` also is not descriptive. I recommend `quant_dirname`. 80. p4: Would suggest more descriptive names than colfile/coldata, orthogonal to whatever naming scheme you use (which means either changes to names here or to use variations of the names I used above) 81. p5: This transcript is a lot to digest, and it's hard to know which parts of the explanatory text refer to which part of the transcript, which has multiple commands that go back a page previous. I suggest breaking it up more with explanatory text to the extent possible. Piping is convenient for a programmer but not for explanation. This comment applies to many of the subsequent transcripts too 82. p5: Get the display of coldata so that it doesn't exceed one line per row 82. p5: Get the display of coldata so that it doesn't exceed one line per row 83. p6: It's not a good idea to call your experiment object merely `se` when it is one of multiple SummarizedExperiments here 84. p6: Please use a consistent style for spacing around = in R code 85. p7: I suggest not giving the same name to different objects (`atac_mat` is first something like a data.frame and then a matrix. 86. p7: Suggest not shadowing built-in functions with variable names like `rownames`. 87. p7: "We know this" antecedent unclear 88. p8: I suggest running all the R code in here through a linter. This will ensure you have consistent style. 89. p8: res is not a specific variable name 90. p10: Why switch between "other genes" and "non DE genes"? Just call them `not_de_genes` from the start. (Or switch to `genes_de` and `genes_not_de`.) 91. p11: Why is library(fluentGenomics) loaded again? 92. p11: Using library(fluentGenomics) adds a data object called `peaks` to your environment? This seems undesirable 93. p12: The text would be a bit more readable if you eschewed the abbreviations "DA", "DE", and "LFC" in the main text. 94. p13: here and elsewhere: integer literals (e.g. 10L) should be used instead of float literals for integers 95. p14: why `2*1e4`? 20000L would be clearer 96. p16: Figure 3: Please specify what the components of the box plot mean. There is no universal box plot.

Are sufficient details provided to allow replication of the method development and its use by others? No
If any results are presented, are all the source data underlying the results available to ensure full reproducibility? No 1.
(as indeed the original study authors must have). But showing the code to generate a PCA would be useful; along with interpretation of the results The code is shown to remove genes with low counts prior to the DESeq analysis. I thought this was no longer necessary as DESeq2 incorporates it's own filtering. Or is this a memory-saving exercise? Some justification of why 10 counts and 6 samples was chosen would be beneficial to those new to Bioconductor and RNA-seq.