Enrichment visualization and differential analysis of chromatin accessibility (ATAC-seq) in MCF7 vs. MCF10A

Mai Tam Pham – 29 March 2025

Objectives

We only examine chr22:17,000,000 – 52,000,000 (GRCh38)

• Task 1:

Detect differential accessibility regions in MCF7 vs. MCF10A

Tools: csaw R/Bioconductor package (https://bioconductor.org/books/release/csawBook/)

• Task 2:

Visualize enrichment ATAC-seq signals of

• Differential accessibility regions in MCF7 vs. MCF10A (Visualization example)
• Transcription start site (TSS) and gene regions (transcription start/end site - TSS/TES) (Real-case example)

Tools:

deeptools (https://deeptools.readthedocs.io/en/latest/)

IGV (https://igv.org/doc/desktop/#)

​

Note: we can do the same thing with ChIP-seq data

​

Datasets

We only examine chr22:17,000,000 – 52,000,000 (hg19 - Gencode)

For task 1:

• Input: Pre-processed .bam/.bai files (2 reps for each cell line)
• Output: Differential accessibility regions (direction – upReg/downReg, logFC, FDR…)

For task 2:

• Input:
• Processed .bam/.bai files (2 reps for each cell line)
• BED3 file contains differential accessibility regions
• BED6 file contain gene regions (TSS/TES) with strand specified
• Output:
• Score matrix
• Profile/heatmap plot

ATAC-seq overview

https://hutchdatascience.org/Choosing_Genomics_Tools/atac-seq-1.html

• Nucleosome-free regions (open chromation)
• Short fragments (50-100bp)
• Active enhancers/promoters, TF-binding sites, TF footprint
• Nucleosome-bound regions
• Long fragments (200-600bp)
• Poised promoters/enhancers
• Closed chromatin (heterochromatin)
• Very low/no reads
• Telomere, centromeres, inactive genes

ATAC-seq analysis workflow

Raw FASTQ

Processed FASTQ

Raw BAM

Processed BAM

• QC - FastQC
• Adapter (Nextera Transposase sequence)/read trimming - Cutadapt

• Alignment – bowtie2

• Remove reads mapping to ChrM + multimapped reads – samtools view
• Mark duplicates – Picard MarkDuplicates
• Remove reads in blacklist regions – bedtools intersect
• Shift read coordinates – alignmentSieve --ATACshift (optional for TFBS analysis)

Downstream (fundamental) analysis

Peak calling

Differential accessibility regions

Enrichment ATAC-seq signal visualization

bamToBigWig

macs2

csaw

deeptools bamCoverage

• deeptools computeMatrix
• deeptools plotProfile/plotHeatmap

Task 1

Task 2

Other (advanced) analysis:

• Peak annotation (MACS2, CHIPseekers)
• Known/De novo motif discovery (HOMER)
• Peak calling/Nucleosome positioning detection (HMMRATAC – HMM-based FL modeling)
• TF footprint (HINT-ATAC – HMM-based modeling)

Code and results – Task 1 – Differential accessibility regions

Step 1: Load sorted and indexed BAM files (require .bai files in the same location as .bam)

Step 2: Prepare design matrix and contrast expression (for differential analysis)

• DiffBind, DBChiP: perform differential binding analysis based on read counts for peak intervals identified by MACS
• csaw: uses read counts of sliding windows across the genome (call differential directly from .bam files)

Note: These tools are designed for calling differential binding of ChIP-seq, but can also use for ATAC-seq

Step 3: Summarize window-level read counts

SE reads

PE reads

• readParam(): specify read features (e.g. minMQS, ignore dup reads…) to be counted
• windowCounts(): use sliding-window approach to count fragments for a set of .bam files
• Choice of window width: 5-10bp for TF, >=150bp for histone ChIP-seq or ATAC-seq (to account for the length of DNA wrapped up nucleosome)
• assay(data): a matrix with col as library (sample), rows as fragment counts for each window
• rowRanges(data): genomic locations of each window
• data$totals: total number of fragment counts of each library

• windowCounts(): fragment counts for background (2000bp in width)
• to estimate background abundance

​

• filterWindowsGlobal(): Compare background abundance (binned) with window-based abundance (data)
• to determine which window is driven by background enrichment
• reduce false positives

​

• log2(3): cutoff filtering option

Filtering options:

By count size, by proportion, by local enrichment, by global enrichment, with negative controls, by prior information

Normalization is necessary to remove library-specific biases (e.g. different sequencing depth) prior to differential analysis

Step 4: Detect differential accessibility regions

Feature

(Window)

Window 1

Tagwise (feature-specific) dispersion follows Chi-squared distribution

NB distribution X~NB(mean,dispersion)

Var(X) = mean(X)+mean^2(X)*1/dispersion

Mean1, var1, dispersion1

Mean2, var2, dispersion2

GLM-NB distribution model is used to account for overdispersion (the nature of count data)

F (QL F-statistic): higher means stronger differential

Code and results – Task 2.1 Visualization example of enrichment signals of differential accessibility regions

deeptools (https://deeptools.readthedocs.io/en/latest/)

Step 1: Prepare score (.bw/.bigwig) files from .bam files (require .bai files at the same location as .bam files

Step 2: Prepare region files (BED3/BED6 with strand specified)

Step 3: Construct a score matrix (rows as regions, columns as bin-based scores)

Step 4: Create heatmap plot based on the computed matrix score

Scenario 1: Enrichment heatmap plot for only 2 samples (MCF7_rep1 and MCF10A_rep1)

Scenario 2: Enrichment heatmap plot for 4 samples (MCF7_rep1, MCF7_rep2, MCF10A_rep1, MCF10A_rep2)

Code and results – Task 2.2 Visualization of ATAC-seq enrichment signals in TSS and gene regions

Scenario 1: Enrichment ATAC-seq signal around TSS

Scenario 2: Enrichment ATAC-seq signal along gene regions (TSS/TES)

Task 2.3 ATAC-seq peak/alignment visualization using IGV

TSS+/-8000bp

TSS+/-500bp

TSS+/-200bp

TSS+/-200bp

References

• https://github.com/CebolaLab/ATAC-seq
• https://github.com/harvardinformatics/ATAC-seq

​