Factor-centric chromatin analysis

Srinivas Ramachandran

RNA Bioscience Initiative | CU Anschutz

2026-09-25

Where do transcription factors bind in the genome?

Today we’ll look at where two yeast transcription factors bind in the genome using CUT&RUN.

Where do transcription factors bind in the genome?

Techniques like CUT&RUN require an affinity reagent (e.g., an antibody) that uniquely recognizes a transcription factor in the cell.

This antibody is added to permeabilized cells, and the antibody associates with the epitope. A separate reagent, a fusion of Protein A (which binds IgG) and micrococcal nuclease (MNase) then associates with the antibody. Addition of calcium activates MNase, and nearby DNA is digested. These DNA fragments are then isolated and sequenced to identify sites of TF association in the genome.

Where do transcription factors bind in the genome?

Fig 1a, Skene et al.

Data download and pre-processing

CUT&RUN data were downloaded from the NCBI GEO page for Skene et al.

I selected the 16 second time point for S. cerevisiae Abf1 and Reb1 (note the paper combined data from the 1-32 second time points).

BED files containing mapped DNA fragments were separated by size and converted to bigWig with:

# separate fragments by size
awk '($3 - $2 <= 120)' Abf1.bed > CutRun_Abf1_lt120.bed
awk '($3 - $2 => 150)' Abf1.bed > CutRun_Abf1_gt150.bed

# for each file with the different sizes
bedtools genomecov -i Abf1.bed -g sacCer3.chrom.sizes -bg > Abf1.bg
bedGraphToBigWig Abf1.bg sacCer3.chrom.sizes Abf1.bw

The bigWig files are available here in the data/ directory.

CUT&RUN analysis

library(tidyverse)
library(here)
library(valr)

# genome viz
library(TxDb.Scerevisiae.UCSC.sacCer3.sgdGene)
library(Gviz)
library(rtracklayer)

# motif discovery and viz
library(BSgenome.Scerevisiae.UCSC.sacCer3)
library(memes)
library(universalmotif)

track_start <- 90000
track_end <- 150000

# genes track
sgd_genes_trk <-
  GeneRegionTrack(
    TxDb.Scerevisiae.UCSC.sacCer3.sgdGene,
    chromosome = "chrII",
    start = track_start,
    end = track_end,
    background.title = "white",
    col.title = "black",
    fontsize = 16
  )

# signal tracks
track_info <-
  tibble(
    file_name = c(
      "CutRun_Reb1_lt120.bw",
      "CutRun_Abf1_lt120.bw",
      "CutRun_Reb1_gt150.bw",
      "CutRun_Abf1_gt150.bw"
    ),
    sample_type = c(
      "Reb1_Short",
      "Abf1_Short",
      "Reb1_Long",
      "Abf1_Long"
    )
  ) |>
  mutate(
    file_path = here("data/block-dna", file_name),
    big_wig = purrr::map(
      file_path,
      \(x) import.bw(x, as = "GRanges")
    ),
    data_track = purrr::map2(
      big_wig,
      sample_type,
      \(x, y) {
        DataTrack(
          x,
          name = y,
          background.title = "white",
          col.title = "black",
          col.axis = "black",
          fontsize = 16
        )
      }
    )
  ) |>
  dplyr::select(sample_type, big_wig, data_track)

# x-axis track
x_axis_trk <- GenomeAxisTrack(
  col = "black",
  col.axis = "black",
  fontsize = 16
)

plotTracks(
  c(
    sgd_genes_trk,
    track_info$data_track,
    x_axis_trk
  ),
  from = track_start,
  to = track_end,
  chromosome = "chrII",
  transcriptAnnotation = "gene",
  shape = "arrow",
  type = "histogram"
)

Genome browser tracks for Abf1 and Reb1 CUT&RUN data.

Peak calling

abf1_tbl <- read_bigwig(here("data/block-dna/CutRun_Abf1_lt120.bw"))
total_reads <- 16e6

genome <- read_genome(here("data/block-dna/sacCer3.chrom.sizes"))
genome_size <- sum(genome$size)

genome_lambda <- total_reads / genome_size
peak_calls <-
  abf1_tbl |>
  # define single-base sites
  mutate(
    midpoint = start + round((end - start) / 2),
    start = midpoint,
    end = start + 1,
    # use the poisson to calculate a p-value with the genome-wide lambda
    pval = dpois(value, genome_lambda),
    # convert p-values to FDR
    fdr = p.adjust(pval, method = "fdr")
  )

Peak calling

peak_calls_sig <-
  filter(
    peak_calls,
    fdr == 0
  ) |>
  # collapse neighboring, significant sites
  bed_merge(max_dist = 20) |>
  # drop narrow "spike" intervals too short for motif discovery
  filter(end - start >= 10)

filter(
  peak_calls_sig,
  chrom == "chrII" &
    start >= track_start &
    end <= track_end
)

# A tibble: 5 × 3
  chrom  start    end
  <chr>  <dbl>  <dbl>
1 chrII 100248 100289
2 chrII 101292 101393
3 chrII 124916 124949
4 chrII 136181 136264
5 chrII 141070 141121

How do proteins recognize specific locations in the genome to bind?

Motif discovery

Theory

There are two major approaches to defining sequence motifs enriched in a sample: enumerative and probabilistic approaches.

Theory

Here we’ll apply a probabilistic approach (MEME) to discover motifs in a collection of DNA sequences. During the RNA block, you’ll learn about k-mer analysis, which is a form of enumerative approach.

In each case, the goal is to define a set of sequence motifs that are encriched in a set of provided sequences (i.e., peaks from CUT&RUN data) relative to a genomic background.

Theory

Motifs are expressed in a Position Weight Matrix, which captures the propensities for a position to be a particular nucleotide in a sequence motif.

These PWMs can be represented as sequence logos, visually represent the amount of information provided by the motif, typically using “information content”, expressed in bits.

Theory

LexA sequence motif

Practice

We’ll use the memes and universalmotif packages from Bioconductor to derive sequence motifs from the peaks we called above.

Note: MEME (the underlying motif discovery software) isn’t installed in this environment, so this step was run ahead of time on a server where it is. The output is provided for you to import.

Collect the DNA sequences within the peak windows using the BSgenome for S. cerevisiae
Provide those sequences to memes::runMeme() (run ahead of time here), which uses an Expectation-Maximization (EM) approach to identify and refine motifs.
Import the pre-computed output with universalmotif::read_meme(), and examine the discovered motifs.
Plot the top motif as a logo using universalmotif::view_motifs().

peak_calls_gr <-
  GRanges(
    seqnames = peak_calls_sig$chrom,
    ranges = IRanges(peak_calls_sig$start, peak_calls_sig$end)
  )

peak_seqs <- BSgenome::getSeq(
  # provided by BSgenome.Scerevisiae.UCSC.sacCer3
  Scerevisiae,
  peak_calls_gr
)

names(peak_seqs) <- as.character(peak_calls_gr)
Biostrings::writeXStringSet(peak_seqs, here("data/block-dna/abf1_peaks.fasta"))

# MEME Suite is not installed in this environment, so this step was run
# ahead of time on a server where it is installed. This is the command
# that was used to generate the motifs (not run here):
#
# meme_results <- memes::runMeme(
#   peak_seqs,
#   outdir = "meme_out",
#   nmotifs = 5,
#   minw = 6,
#   maxw = 20
# )

# import the pre-computed MEME output provided for this exercise
meme_motifs <- universalmotif::read_meme(
  here("data/block-dna/abf1_meme_out/meme.txt")
)

# read_meme() returns a single motif object (not a list) when there's only
# one motif -- normalize to a list either way
if (!is.list(meme_motifs)) meme_motifs <- list(meme_motifs)

# convert to a data frame for easy inspection (sorted by significance)
meme_df <- universalmotif::to_df(meme_motifs)

# look at the consensus motifs
meme_df

Plot the logo

Now let’s look at the sequence logo for the top hit.

# meme_df is sorted by significance, so the first row is the top hit
top_motif <- meme_df$motif[[1]]

universalmotif::view_motifs(top_motif)

Plot the logo

Questions

Does this motif make sense, based on what you know about the requirements and specificity of DNA binding by transcription factors?
How might you confirm that a specific sequence (that conforms to a motif) is bound directly by a transcription factor?