Overview

valr includes several functions for exploring statistical relationships between sets of intervals.

In this vignette we explore the relationship between transcription start sites and repetitive elements in the human genome.

Distance metrics

First we define a function that takes x and y intervals and computes distance statistics (using bed_reldist() and bed_absdist()) for specified groups. The value of each statistic is assigned to a .value column.

distance_stats <- function(x, y, genome, group_var, type = NA) {
  group_by_(x, .dots = group_var) %>%
    do(
       reldist = bed_reldist(., y, detail = TRUE) %>%
         select(.value = .reldist),
       absdist = bed_absdist(., y, genome) %>%
         select(.value = .absdist)
       ) %>%
    tidyr::gather_('stat', 'value', setdiff(names(.), list(group_var))) %>%
    mutate(type = type)
} 

We use the distance_stats() function to apply the bed_absdist() function to each group of data.

And the same is done for a set of shuffled group of data. bed_shuffle() is used to shuffle coordinates of the repeats within each chromosome (i.e., the coordinates change, but the chromosome stays the same.)

Now we can bind the observed and shuffled data together, and do some tidying to put the data into a format appropriate for a statistical test. This involves:

  1. unnest()ing the data frames
  2. creating groups for each repeat (name), stat (reldist or absdist) and type (obs or shf)
  3. adding unique surrogate row numbers for each group
  4. using tidyr::spread() to create two new obs and shuf columns
  5. removing rows with NA values.

Now that the data are formatted, we can use the non-parametric ks.test() to determine whether there are significant differences between the observed and shuffled data for each group. broom::tidy() is used to reformat the results of each test into a data_frame, and the results of each test are gathered to into a type column for each test type.

library(broom)

pvals <- res %>% do(twosided = tidy(ks.test(.$obs, .$shuf)),
                    less = tidy(ks.test(.$obs, .$shuf, alternative = 'less')),
                    greater = tidy(ks.test(.$obs, .$shuf, alternative = 'greater'))) %>%
  tidyr::gather(alt, type, -name, -stat) %>%
  unnest(type) %>%
  select(name:p.value) %>%
  arrange(p.value)

Histgrams of the different stats help visulaize the distribution of p.values.

ggplot(pvals, aes(p.value)) +
  geom_histogram(binwidth = 0.05) +
  facet_grid(stat ~ alt) + theme_cowplot()

We can also assess false discovery rates (q.values) using p.adjust().

pvals <-
  group_by(pvals, stat, alt) %>%
  mutate(q.value = p.adjust(p.value)) %>%
  ungroup() %>%
  arrange(q.value)

Finally we can visualize these results using stat_ecdf().

res_gather <- tidyr::gather(res, type, value, -name, -stat, -.id)

signif <- head(pvals, 5)

res_signif <-
  signif %>%
  left_join(res_gather, by = c('name', 'stat'))

ggplot(res_signif, aes(x = value, color = type)) +
  stat_ecdf() + 
  facet_grid(stat ~ name) + theme_cowplot() + scale_x_log10() +
  scale_color_brewer(palette = 'Set1')

Projection test

bed_projection() is a statistical approach to assess the relationship between two intervals based on the binomial distribution. Here, we examine the distribution of repetitive elements within the promoters of coding or non-coding genes.

First we’ll extract 5 kb regions upstream of the transcription start sites to represent the promoter regions for coding and non-coding genes.

# create intervals 5kb upstream of tss representing promoters
promoters <-
  bed_flank(genes, genome, left = 5000, strand = TRUE) %>%
  mutate(name = ifelse(grepl('NR_', name), 'non-coding', 'coding')) %>%
  select(chrom:strand)
  
# select coding and non-coding promoters
promoters_coding <- filter(promoters, name == 'coding')
promoters_ncoding <- filter(promoters, name == 'non-coding')

promoters_coding
#> # A tibble: 973 x 6
#>    chrom    start      end name   score strand
#>    <chr>    <int>    <int> <chr>  <chr> <chr> 
#>  1 chr22 16287937 16292937 coding 0     -     
#>  2 chr22 16449804 16454804 coding 0     -     
#>  3 chr22 17073700 17078700 coding 0     -     
#>  4 chr22 17302589 17307589 coding 0     -     
#>  5 chr22 17302589 17307589 coding 0     -     
#>  6 chr22 17489112 17494112 coding 0     -     
#>  7 chr22 17560848 17565848 coding 0     +     
#>  8 chr22 17560848 17565848 coding 0     +     
#>  9 chr22 17602213 17607213 coding 0     -     
#> 10 chr22 17602257 17607257 coding 0     -     
#> # ... with 963 more rows

promoters_ncoding
#> # A tibble: 294 x 6
#>    chrom    start      end name       score strand
#>    <chr>    <int>    <int> <chr>      <chr> <chr> 
#>  1 chr22 16152078 16157078 non-coding 0     +     
#>  2 chr22 16157065 16162065 non-coding 0     +     
#>  3 chr22 16157065 16162065 non-coding 0     +     
#>  4 chr22 16193009 16198009 non-coding 0     -     
#>  5 chr22 16231289 16236289 non-coding 0     -     
#>  6 chr22 16269608 16274608 non-coding 0     +     
#>  7 chr22 17077800 17082800 non-coding 0     +     
#>  8 chr22 17156430 17161430 non-coding 0     -     
#>  9 chr22 17229328 17234328 non-coding 0     -     
#> 10 chr22 17303363 17308363 non-coding 0     +     
#> # ... with 284 more rows

Next we’ll apply the bed_projection() test for each repeat class for both coding and non-coding regions.

The projection test is a two-tailed statistical test. A significant p-value indicates either enrichment or depletion of query intervals compared to the reference interval sets. A value of lower_tail = TRUE column indicates that the query intervals are depleted, whereas lower_tail = FALSE indicates that the query intervals are enriched.