Fluent genomics with plyranges and tximeta (original) (raw)
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%.
da_peaks <- peaks %>%
filter(da_padj < 0.01)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?
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 inother_genes:
size <- length(de_genes)
slice(other_genes, sample.int(n(), size))## GRanges object with 822 ranges and 3 metadata columns:
## seqnames ranges strand | gene_id de_log2FC
## <Rle> <IRanges> <Rle> | <character> <numeric>
## [1] chr16 566995-584136 + | ENSG00000007541.16 -0.5232899
## [2] chr2 23385217-23708611 + | ENSG00000119771.14 0.0722953
## [3] chr12 111405948-111451623 + | ENSG00000111252.10 -0.5550827
## [4] chr9 97674909-97697357 - | ENSG00000136936.10 -0.4198125
## [5] chr8 143579697-143599541 - | ENSG00000104529.17 -0.1558568
## ... ... ... ... . ... ...
## [818] chr1 155308748-155320666 + | ENSG00000160752.14 -0.468680
## [819] chr13 41566837-41961120 - | ENSG00000102763.17 0.333812
## [820] chrX 2903970-2929351 - | ENSG00000006756.15 0.123042
## [821] chr19 984271-998438 + | ENSG00000065268.10 -0.148756
## [822] chr8 16107878-16567490 - | ENSG00000038945.14 0.104844
## de_padj
## <numeric>
## [1] 4.10137e-03
## [2] 8.95877e-03
## [3] 8.00407e-03
## [4] 1.82673e-04
## [5] 7.81870e-18
## ... ...
## [818] 3.31400e-05
## [819] 2.74332e-03
## [820] 5.60234e-13
## [821] 2.99917e-04
## [822] 1.50568e-04
## -------
## seqinfo: 25 sequences (1 circular) from hg38 genomeWe 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.
# set a seed for the results
set.seed(2019-08-02)
boot_genes <- replicate(10,
slice(other_genes, sample.int(n(), size)),
simplify = FALSE)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## GRanges object with 9042 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
## ... ... ... ... . ... ...
## [9038] chr2 197705369-197708387 + | ENSG00000247626.4 -0.0862514
## [9039] chr5 90515611-90529584 - | ENSG00000176018.12 0.0161995
## [9040] chr12 6666477-6689572 - | ENSG00000126746.17 0.0158683
## [9041] chr19 984271-998438 + | ENSG00000065268.10 -0.1487556
## [9042] chr12 47782722-47833132 - | ENSG00000061273.17 -0.3302549
## de_padj resample origin
## <numeric> <integer> <factor>
## [1] 6.85285e-14 0 de
## [2] 1.58739e-05 0 de
## [3] 1.00655e-11 0 de
## [4] 5.36800e-09 0 de
## [5] 2.41452e-18 0 de
## ... ... ... ...
## [9038] 1.25864e-04 10 not_de
## [9039] 6.59016e-05 10 not_de
## [9040] 1.15330e-13 10 not_de
## [9041] 2.99917e-04 10 not_de
## [9042] 9.17085e-05 10 not_de
## -------
## seqinfo: 25 sequences (1 circular) from hg38 genomeExpanding 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 ofstart 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.
all_genes <- all_genes %>%
anchor_5p() %>%
mutate(width = 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.
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.
genes_olap_peaks <- all_genes %>%
join_overlap_left(da_peaks)
genes_olap_peaks## GRanges object with 30449 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
## ... ... ... ... . ... ...
## [30445] chr12 47823132-47843131 - | ENSG00000061273.17 -0.330255
## [30446] chr12 47823132-47843131 - | ENSG00000061273.17 -0.330255
## [30447] chr12 47823132-47843131 - | ENSG00000061273.17 -0.330255
## [30448] chr12 47823132-47843131 - | ENSG00000061273.17 -0.330255
## [30449] chr12 47823132-47843131 - | ENSG00000061273.17 -0.330255
## de_padj resample origin peak_id da_log2FC da_padj
## <numeric> <integer> <factor> <character> <numeric> <numeric>
## [1] 6.85285e-14 0 de ATAC_peak_21236 -0.546582 1.15274e-04
## [2] 1.58739e-05 0 de ATAC_peak_231183 1.453297 9.73225e-17
## [3] 1.00655e-11 0 de ATAC_peak_193578 0.222371 3.00939e-11
## [4] 1.00655e-11 0 de ATAC_peak_193579 -0.281615 7.99889e-05
## [5] 1.00655e-11 0 de ATAC_peak_193580 0.673705 7.60043e-15
## ... ... ... ... ... ... ...
## [30445] 9.17085e-05 10 not_de ATAC_peak_61944 -0.864385 2.15479e-10
## [30446] 9.17085e-05 10 not_de ATAC_peak_61945 -1.263839 2.83702e-21
## [30447] 9.17085e-05 10 not_de ATAC_peak_61946 -0.854639 2.02495e-16
## [30448] 9.17085e-05 10 not_de ATAC_peak_61947 -0.421568 5.22070e-04
## [30449] 9.17085e-05 10 not_de ATAC_peak_61948 -0.538309 1.46133e-08
## -------
## seqinfo: 25 sequences (1 circular) from hg38 genomeNow 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 viareduce_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:
library(ggplot2)
gene_peak_max_lfc %>%
filter(peak_count > 0) %>%
as.data.frame() %>%
ggplot(aes(origin, peak_max_lfc)) +
geom_boxplot()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:
origin_peak_lfc <- genes_olap_peaks %>%
group_by(origin) %>%
summarize(
peak_count = sum(!is.na(da_padj)) / n_distinct(resample),
lfc1_peak_count =sum(abs(da_log2FC) > 1, na.rm=TRUE)/ n_distinct(resample),
lfc2_peak_count = sum(abs(da_log2FC) > 2, na.rm=TRUE)/ n_distinct(resample)
)
origin_peak_lfc## DataFrame with 2 rows and 4 columns
## origin peak_count lfc1_peak_count lfc2_peak_count
## <factor> <numeric> <numeric> <numeric>
## 1 not_de 2625.7 406.8 35.7
## 2 de 3726.0 1181.0 246.0Here 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:
origin_peak_lfc %>%
as.data.frame() %>%
tidyr::pivot_longer(cols = -origin) %>%
tidyr::pivot_wider(names_from = origin, values_from = value) %>%
mutate(enrichment = de / not_de)## # A tibble: 3 × 4
## name not_de de enrichment
## <chr> <dbl> <dbl> <dbl>
## 1 peak_count 2626. 3726 1.42
## 2 lfc1_peak_count 407. 1181 2.90
## 3 lfc2_peak_count 35.7 246 6.89The 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).
genes_olap_peaks %>%
group_by(gene_id, origin, resample) %>%
reduce_ranges_directed(
lfc1 = sum(abs(da_log2FC) > 1, na.rm=TRUE),
lfc2 = sum(abs(da_log2FC) > 2, na.rm=TRUE)
) %>%
group_by(origin) %>%
summarize(
lfc1_gene_count = sum(lfc1 > 0) / n_distinct(resample),
lfc1_peak_count = sum(lfc1) / n_distinct(resample),
lfc2_gene_count = sum(lfc2 > 0) / n_distinct(resample),
lfc2_peak_count = sum(lfc2) / n_distinct(resample)
)## 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 300 406.8 33.4 35.7
## 2 de 564 1181.0 162.0 246.0To do this for many thresholds is cumbersome and would create a lot of duplicate code. Instead we create a single function calledcount_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 thethresholds 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 calledvalue. 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.
genes_peak_all_thresholds <- genes_olap_peaks %>%
group_by(gene_id, origin, resample) %>%
reduce_ranges_directed(
value = count_if_above_threshold(da_log2FC, thresholds),
threshold = list(thresholds)
)
genes_peak_all_thresholds## GRanges object with 9042 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
## ... ... ... ... . ... ...
## [9038] chr2 197695369-197715368 + | ENSG00000247626.4 not_de
## [9039] chr5 90519584-90539583 - | ENSG00000176018.12 not_de
## [9040] chr12 6679572-6699571 - | ENSG00000126746.17 not_de
## [9041] chr19 974271-994270 + | ENSG00000065268.10 not_de
## [9042] chr12 47823132-47843131 - | ENSG00000061273.17 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,...
## ... ... ... ...
## [9038] 10 6,5,5,... 0.0658243,0.1184840,0.1711436,...
## [9039] 10 2,2,1,... 0.0658243,0.1184840,0.1711436,...
## [9040] 10 2,2,2,... 0.0658243,0.1184840,0.1711436,...
## [9041] 10 3,3,3,... 0.0658243,0.1184840,0.1711436,...
## [9042] 10 7,7,7,... 0.0658243,0.1184840,0.1711436,...
## -------
## seqinfo: 25 sequences (1 circular) from hg38 genomeNow we can expand these list-columns into a long GRanges object using theexpand_ranges() function. This function will unlist the value andthreshold columns and lengthen the resulting GRanges object. To compute the peak and gene counts for each threshold, we apply the same summarization as before:
origin_peak_all_thresholds <- genes_peak_all_thresholds %>%
expand_ranges() %>%
group_by(origin, threshold) %>%
summarize(
gene_count = sum(value > 0) / n_distinct(resample),
peak_count = sum(value) / n_distinct(resample)
)
origin_peak_all_thresholds## DataFrame with 200 rows and 4 columns
## origin threshold gene_count peak_count
## <factor> <numeric> <numeric> <numeric>
## 1 not_de 0.0658243 778.1 2625.3
## 2 not_de 0.1184840 768.0 2548.1
## 3 not_de 0.1711436 753.8 2391.3
## 4 not_de 0.2238033 739.3 2183.3
## 5 not_de 0.2764629 714.6 1961.2
## ... ... ... ... ...
## 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 0Again 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:
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 rows containing missing values or values outside the scale range
## (`geom_line()`).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.
origin_threshold_counts %>%
ggplot(aes(x = threshold,
y = count + 1,
color = origin,
linetype = type)) +
geom_line() +
scale_y_log10()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.