Factor-centric chromatin analysis

Author

Jay Hesselberth

Published

October 21, 2024

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.

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.

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)

Warning: replacing previous import 'S4Arrays::makeNindexFromArrayViewport' by
'DelayedArray::makeNindexFromArrayViewport' when loading 'SummarizedExperiment'

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

# motif discovery and viz
library(BSgenome.Scerevisiae.UCSC.sacCer3)
library(rGADEM)
library(seqLogo)

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, ~ import.bw(.x, as = "GRanges")
    ),
    data_track = purrr::map2(
      big_wig, sample_type,
      ~ 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"
)

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(score, genome_lambda),
    # convert p-values to FDR
    fdr = p.adjust(pval, method = "fdr")
  )

peak_calls_sig <-
  filter(
    peak_calls,
    fdr == 0
  ) |>
  # collapse neighboring, significant sites
  bed_merge(max_dist = 20)

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

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

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

Motif discovery

Theory

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

Here we’ll apply a probabilistic approach (GADEM) 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.

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.

Practice

We’ll use the rGADEM package from Bioconductor to derive sequence motifs from the peaks we called above. This is a straightforward process:

Collect the DNA sequences within the peak windows using the BSgenome for S. cerevisiae
Provide those sequences and the genomic background to rGADEM::GADEM(), which runs uses an Expectation-Maximization (EM) approach to identify and refine motifs.
Examine the discovered motifs, and plot as a logo using seqLogo::seqLogo().

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

# takes ~2 minutes to run
gadem <- rGADEM::GADEM(
  peak_seqs,
  genome = Scerevisiae,
  verbose = 1
)

# look at the consensus motifs
consensus(gadem)

# how many consensus motifs are there?
nOccurrences(gadem)

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

pwm <- gadem@motifList[[1]]@pwm

seqLogo::seqLogo(pwm)

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?

References

GADEM: a genetic algorithm guided formation of spaced dyads coupled with an EM algorithm for motif discovery. J Comput Biol 2009 [PMC free article] [PubMed] [Google Scholar]