Pre-filtering database for expressed transcription factors

A key consideration when using a database of known transcription factors is whether to include motifs for TFs which are not expressed. On one hand, including all motifs allows maximum discovery of potential regulators, but on the other hand may detect factors which are not relevant under your experimental conditions. Another consideration for many tools is that increasing the number of tested motifs increases the multiple testing penaly, so limiting the motif database to a smaller number of relevant candidates can be a way to improve statistical power. There is no one right answer to this question. For this vignette, I will give an example for how one could go about filtering a database to consider only expressed factors. For other analyses, this question will need to be carefully considered.

First, we load the motif database into R using read_meme. To more easily manipulate the motif metadata, we can convert this to a data.frame using to_df.

library(universalmotif)
meme_db <- read_meme(meme_db_path) %>% 
  to_df()

To simplify the vignette, we’ve preproccessed RNAseq data from our tissue of interest and shipped them with memes. These are FPKM normalized counts for each transcription factor gene at early and late developmental times.

data("example_rnaseq", package = "memes")
head(example_rnaseq, 4)

Choosing a method to decide whether a gene is expressed or not is a non-trivial topic. DO NOT BLINDLY DO WHAT I DO HERE. To keep this vignette on topic, I will use an overly simplistic cutoff of FPKM > 5 to define expressed genes. In this example, I will consider any gene that exceeds the cutoff at either timepoint as expressed.

expressed_genes <- example_rnaseq %>% 
  # For each gene, keep only those with max(FPKM) greater or equal to 5.
  dplyr::group_by(symbol) %>% 
  dplyr::filter(max(fpkm) >= 5)

Finally, I will filter the full motif database to select only those motifs corresponding to expressed gene.

meme_db_expressed <- meme_db %>% 
  # the altname slot of meme_db contains the gene symbol
  # (this is database-specific)
  dplyr::filter(altname %in% expressed_genes$symbol)

Below is a comparison of the number of motifs in the database before and after filtering.

#Number of motifs pre-filtering: 
nrow(meme_db)
#Number of motifs post-filter: 
nrow(meme_db_expressed)

Using this filtering strategy, we’ve decreased our multiple-testing penalty by ~54%.

As a convenience function, we can set a default motif database to use for each memes function by setting the R option meme_db to a universalmotif object. This can be done as follows:

# to_list() converts the database back from data.frame format to a standard `universalmotif` object.
options(meme_db = to_list(meme_db_expressed, extrainfo = FALSE))

Alternatively, you can pass a universalmotif object to the database parameter for any memes function and this will override the global option.

Examination of binding categories with AME

Next, to test for enrichment of known motifs between E93 binding categories, we first collect the sequences of each category like so:

by_binding <- summit_flank %>%
  # Get a list of chip peaks belonging to each set
  split(mcols(.)$peak_binding_description) %>%
  # look up the DNA sequence of each peak within each group
  get_sequence(dm.genome)

The data are returned as a Biostrings List, where each list entry represents the sequences of each E93 binding category.

head(by_binding)

Finally, we test each set of sequencing using AME with the runAme function. This will run AME on each set of input sequences in the Biostrings List.

ame_by_binding <- by_binding %>% 
  runAme

This returns a list object where each entry is the AME results for each category. For example, here is what the results for the ectopic peaks look like:

head(ame_by_binding$ectopic, 5)

Visualizing AME results

The plot_ame_heatmap() function provides a quick way to visualize AME results. It is built on top of ggplot2, so all ggplot2 functions can be used to further modify the plot.

By default, it uses the -log10(adjusted p-value) as the heat values. See the documentation (?plot_ame_heatmap) for additional notes on customization.

We can visualize the top 10 hits from the ectopic binding results as follows:

ame_by_binding$ectopic %>% 
  dplyr::filter(rank %in% 1:10) %>% 
  plot_ame_heatmap(group_name = "Ectopic Sites")  +
    ggtitle("Top 10 AME Hits in Ectopic Sites")

To plot results from multiple runs together, they must first be joined into a single data frame. The ame_by_binding object is a list whose names correspond to the E93 binding category. The list can be combined into a data.frame using dplyr::bind_rows. Setting .id = "binding_type" creates a new column binding_type that contains the names from the ame_by_binding list. In this way, the ame_res data.frame contains all AME results for each run, which can be distinguished by the binding_type column.

ame_res <- ame_by_binding %>% 
  dplyr::bind_rows(.id = "binding_type")
  

It is possible to aggregate results from multiple runs into a heatmap by setting the group parameter in plot_ame_heatmap(). This will stratify the y-axis by the entries in the column passed to group. This results in clustering motifs with shared hits across categories along the x-axis for quickly visualizing motifs that may be shared or distinct between multiple groups.

ame_res %>% 
  plot_ame_heatmap(group = binding_type)

Reducing redundant motif hits

Another key consideration for the above visualization is that in the FlyFactorSurvey database we used, different TFs can have multiple motif entries in the database which are all detected separately by AME. Here, when returning the top 5 hits from each group, you can see, for example, that a motif matching “Aef1” is reported 2 times within the top 5 hits of the ectopic sites. In this situation, it makes sense to summarize the data at the TF level, instead of the motif level. Note: There may be exceptions to this if, for example, a TF has multiple DNA binding sequences it can recognize, in which case having multiple hits may reflect a biological property of your sequences. You will have to handle this on a case-by-case basis for interesting hits and different motif databases. Here, we can see that at least for Aef1, the consensus sequences are very similar.

ame_res %>% 
  dplyr::filter(binding_type == "ectopic", rank %in% 1:5) %>% 
  head(5) %>% 
  dplyr::select(binding_type, rank, motif_id, motif_alt_id, consensus)

How to solve this problem will vary with different motif databases (For details on how to pre-process a motif database, see vignette("tidy_motifs")). In the tidy version of the FlyFactorSurvey database, the altnames of the motifs are set to the transcription factor gene symbol. This information is included in the AME results as the motif_alt_id column. In order to reduce the occurrance of redundant hits in the heatmap, we select the hit for each TF with the lowest p-value (i.e. the most significant hit) as follows:

ame_res %>%
  # perform the next dplyr operation on each TF within each binding type
  dplyr::group_by(binding_type, motif_alt_id) %>%
  # within each binding type, select the TF hit with the lowest adjusted p-value
  dplyr::filter(adj.pvalue == min(adj.pvalue)) %>% 
  plot_ame_heatmap(group = binding_type, id = motif_alt_id) +
    labs(y = "Binding Category",
         x = "Transcription Factor Motif")

AME Heatmap Visualization

If you have read vignette("core_ame"), the “normalized rank” heatmap was introduced as an alternative visualization method useful when AME produces very large numbers of hits, or when p-values are on very different scales between groups. The ame_compare_heatmap_methods() function can be used to visualize why.

Below is a comparison of the distribution of values when using -log10(adj.pvalue) (A) vs normalized ranks (B). Because there are relatively few hits in the results (~30), and the number of hits between groups varies more than the -log10(p-value) distributions, the “normalize” method will produce a misleading heatmap relative to the -log10(p-value) map.

ame_res %>% 
  dplyr::group_by(binding_type, motif_alt_id) %>% 
  dplyr::filter(adj.pvalue == min(adj.pvalue)) %>% 
  ame_compare_heatmap_methods(group = binding_type)

A third option for adjusting this visualization is to cap the heatmap scale values at a certain value. A good rule of thumb for selecting this cutoff is to view the cumulative distribution plot above in (A), and select a value on the x-axis which captures a majority of the data. Here, we see that -log10(adj.pvalue) < 10 captures >90% of the entopic & ectopic hits, but only ~50% of the orphan hits. By selecting 10 as the heatmap scale cap, we can improve the dynamic range of the signal for values <10, while preserving the information that many orphan sites generally have higher scores than the other categories. The heatmap scale can be capped at 10 by setting scale_max = 10. Below is a comparison of the 3 different visualizations. Notice how the intensity of the orphan vs other categories hits changes.

best_ame_hits <- ame_res %>% 
  dplyr::group_by(binding_type, motif_alt_id) %>% 
  dplyr::filter(adj.pvalue == min(adj.pvalue))
  
pval_heatmap <- best_ame_hits %>% 
  plot_ame_heatmap(group = binding_type, id = motif_alt_id) +
    labs(x = NULL,
         title = "-log10(adj.pvalue)")
  

norm_heatmap <- best_ame_hits %>% 
  plot_ame_heatmap(group = binding_type, id = motif_alt_id, value = "normalize") +
    labs(x = NULL,
         title = "normalize")

pval_scaled_heatmap <- best_ame_hits %>% 
  plot_ame_heatmap(group = binding_type, id = motif_alt_id, scale_max = 10) +
    labs(x = NULL,
         title = "-log10(adj.pvalue) (scale capped at 10)")

cowplot::plot_grid(pval_heatmap,
                   norm_heatmap,
                   pval_scaled_heatmap,
                   ncol = 1,
                   labels = "AUTO")

Test de-novo motif enrichment using AME

The above analysis suggests the motif content of the different binding categories are highly similar in sequence composition. To extend these analyses, we can use AME to test for motif enrichment of the de-novo discovered motifs within each binding category, and determine whether the motifs detected in one category are indeed enriched in another. To do this, we can provide the de-novo motifs as the AME database to test for their enrichment in each sequence category.

runAme() allows using a runDreme() results object as the database input by passing it within a list(). Naming the list() entry produces an informative motif_db name in the results data.frame.

ame_denovo_by_binding <- by_binding %>% 
  runAme(database = list("denovo_binding_motifs" = dreme_by_binding_renamed)) %>% 
  dplyr::bind_rows(.id = "binding_type") 

Plotting the heatmap of results reveals that indeed a majority of the de-novo motifs discovered within a single category are detected in all 3 categories, supporting the conclusion that orphan, ectopic, and entopic sites are highly similar in sequence content.

ame_denovo_by_binding %>% 
  plot_ame_heatmap(group = binding_type, scale_max = 10)

However, there are 2 interesting motifs which distinguish orphan and ectopic sites from entopic sites. To help identify which TFs these motifs might belong to, we can use TomTom to match them to known TF motifs.

First, we select the motif id’s that are not found in entopic sites.

entopic_motifs <- ame_denovo_by_binding %>% 
  dplyr::filter(binding_type == "entopic") %>% 
  dplyr::pull(motif_id)

ame_denovo_binding_unique <- ame_denovo_by_binding %>% 
  dplyr::filter(!(motif_id %in% entopic_motifs))

Next, we use the motif id’s from the unique AME results to select those entries in the Dreme results object, and run TomTom on that subset.

dreme_by_binding_unique <- dreme_by_binding_renamed %>% 
  dplyr::filter(name %in% ame_denovo_binding_unique$motif_id) %>% 
  runTomTom(dist = "ed")

Finally, we visualize the TomTom results to identify candidate TFs driving the presence of ectopic and orphan sites.

dreme_by_binding_unique %>% 
  view_tomtom_hits(3) %>%
  cowplot::plot_grid(plotlist = ., nrow = 1, labels = "AUTO", byrow = TRUE)

Counting the number of motifs per peak

plyranges extends dplyr-like syntax to range objects. Here we add a count for each motif per peak in the n_motifs column. We also add a column has_motif which will be a binary indicator of whether a peak contains any motifs.

summit_flank %<>% 
  # Note, if running FIMO with multiple motifs, this solution will not work
  # as it will count all motifs within the fimo-results without splitting by motif_id
  plyranges::mutate(n_motifs = plyranges::count_overlaps(., fimo_res), 
                    has_motif = n_motifs > 0)

First, we want to determine whether E93 sensitive sites are more likely to have E93 motifs in certain response types. Here we can see that sensitive decreasing sites are more likely to have E93 motifs than sensitive increasing or insensitive static sites.

summit_flank %>% 
  data.frame %>% 
  dplyr::mutate(has_match = dplyr::if_else(has_motif, "Match", "No Match")) %>% 
  ggplot(aes(e93_sensitive_behavior)) +
    geom_bar(aes(fill = forcats::fct_rev(has_match)), position = "fill") +
    scale_fill_manual(values = c("Match" = "firebrick",
                                 "No Match" = "Black")) +
    labs(fill = "E93 Motif Match",
         y = "Fraction of Sites",
         x = "Response to E93 binding")

To investigate whether E93-sensitive sites that increase, decrease, or remain static have different numbers of E93 motifs, we plot the fraction of sites with each number of motifs per group. Here we can see that in addition to being more likely to contain an E93 motif, sensitive decreasing sites are more likely to to contain 2 or more matches, where 10% contain at least 2 motifs.

summit_flank %>% 
  # currently, group operations are faster as data.frames, so we convert to data.frame
  data.frame %>%
  dplyr::group_by(e93_sensitive_behavior, n_motifs) %>% 
  dplyr::count() %>% 
  dplyr::group_by(e93_sensitive_behavior) %>% 
  dplyr::mutate(frac = n/sum(n)) %>% 
  ggplot(aes(n_motifs, frac)) +
    geom_line(aes(color = e93_sensitive_behavior), size = 1) +
    labs(y = "Fraction of Sites",
         x = "Number of E93 Motifs",
         color = "Response to E93 Binding") +
    theme_bw()

Finally, we want to assess whether the quality of E93 motifs is different between sensitivity categories. To examine this, we need to determine which motifs are found in which peaks. We use plyranges::join_overlap_intersect to return motif entries appended with peak-level metadata, like the peak id each motif is found within.

# return position of each motif match w/ peak metadata
intersect <- fimo_res %>% 
  plyranges::join_overlap_intersect(summit_flank)

We use the FIMO score as a proxy for quality, where higher scores are better matches to the motif. Here we examine only the best match (highest score) motif within each peak.

best_motifs <- intersect %>%
  # group by ChIP peak id
  plyranges::group_by(id) %>% 
  # Keep only motifs w/ the highest score within a peak
  plyranges::filter(score == max(score)) %>% 
  plyranges::ungroup()

Finally, visualize the FIMO scores for the best E93 motif within each ChIP peak across categories. This reveals that Decreasing sites appear to have higher FIMO scores, suggesting the E93 motifs within Decreasing sites are higher quality E93 motifs than those in Increasing or Static sites.

best_motifs %>% 
  data.frame %>% 
  ggplot(aes(e93_sensitive_behavior, score)) +
    geom_boxplot(aes(fill = e93_sensitive_behavior), notch = TRUE, width = 0.5) +
    guides(fill = "none") +
    labs(x = "Response to E93",
         y = "FIMO Score") +
    theme_bw()

We can test whether the distribution of FIMO scores between multiple groups is statistically significantly different using an Anderson-Darling test. For a pairwise comparison (in case you only have 2 groups to test), a Kolmogorov-Smirnov test is commonly used instead. The PMCMRplus::adAllPairsTest() performs the Anderson-Darling test and corrects for multiple testing. This reveals that Decreasing FIMO scores are significantly different from Increasing and Static scores (p < 0.05), but Increasing and Static scores are not significantly different from eachother (p > 0.05), supporting the conclusion that E93 motifs in Decreasing sites are higher quality than those in Increasing or Static sites.

best_motifs %>% 
  data.frame %>% 
  dplyr::mutate(behavior = factor(e93_sensitive_behavior)) %>% 
  PMCMRplus::adAllPairsTest(score ~ behavior, ., p.adjust.method = "fdr")

We can visualize the differences in E93 motif quality between groups by reconstructing PCM’s from the matched motif sequences. Here, it is clear that although each category contains an E93 motif match, the sequence content of those motifs differs between categories, providing a nice visual representation of the statistical analysis above.

best_motifs %>% 
  # Get sequences of each matched motif
  add_sequence(dm.genome) %>% 
  # Split by E93 response category
  split(mcols(.)$e93_sensitive_behavior) %>% 
  # Convert from GRangesList to list() to use `purrr::imap`
  as.list() %>% 
  # imap passes the list entry as .x and the name of that object to .y
  # here, .y holds the e93_sensitive_behavior name, and .x holds the motif positions
  purrr::imap(~{
    # This builds a new motif using the matched sequences from each group
    create_motif(.x$sequence, 
                 # Append the response category to the motif name
                 name = paste0("E93_", .y))
    }) %>% 
  # Prepend E93 Fly Factor motif to beginning of motif list so it is plotted at the top
  c(list("E93_FlyFactor" = e93_flyfactor), .) %>% 
  # Plot each motif
  view_motifs()

We can also build a more complex representation of these data by making a heatmap of the sequence variability for each matched sequence using plot_sequence_heatmap(). This function will plot each sequence as a row of a heatmap where each column represents the position relative to the motif start coordinate. The cells are colored by the DNA sequence found at that position for each sequence. Using this type of plot it becomes even more clear how variability at position 6 differs between Decreasing and Increasing sites, where Decreasing sites have a very consistent “C” at position 6, while Increasing sites are more likely to have non-C bases at this position.

best_motifs %>% 
  # Get sequences of each matched motif
  add_sequence(dm.genome) %>% 
  data.frame %>% 
  # Split by E93 response category
  split(.$e93_sensitive_behavior) %>% 
  # Extract the "sequence" column from each category
  purrr::map("sequence") %>% 
  # Plot a heatmap for each category
  plot_sequence_heatmap(title_hjust = 0.5) %>% 
  # Arrange into a grid
  cowplot::plot_grid(plotlist = ., nrow = 1)

Centrality of E93 motif

Next we want to visualize whether the E93 motif has altered positioning across the response categories. To do this we add the metadata for each nearest motif to our peak summits. Setting distance = TRUE in plyranges::join_nearest adds a distance column indicating the distance to the nearest joined region.

summit_nearest_e93_motif <- summit_flank %>% 
  plyranges::anchor_center() %>% 
  plyranges::mutate(width = 1) %>% 
  plyranges::join_nearest(fimo_res, distance = TRUE)

We can visualize these data as a cumuluative distribution plot. In this example, it does appear that decreasing sites tend to be closer to the peak summit than increasing or static sites.

summit_nearest_e93_motif %>% 
  data.frame %>%
  # Limit our analysis to peaks that contain motifs
  dplyr::filter(has_motif == TRUE) %>% 
  ggplot(aes(distance)) +
    stat_ecdf(aes(color = e93_sensitive_behavior), size = 1, pad = FALSE) +
    labs(x = "Distance to nearest E93 Motif",
         y = "Fraction of Sites",
         color = "E93 Response") +
    theme_linedraw()