Skip to main content

STEP 8: Analysis

This page covers advanced downstream analysis of ATAC-seq and ChIP-seq data, including peak annotation, motif discovery, differential analysis, and functional interpretation.

Motif Analysis

Using HOMER

Find enriched transcription factor binding motifs:

module load homer

# Find motifs in peaks (de novo + known)
findMotifsGenome.pl peaks.bed mm39 motif_output/ \
  -size 200 \
  -mask

# Use background regions
findMotifsGenome.pl peaks.bed mm39 motif_output/ \
  -size 200 \
  -bg background.bed

# Find specific motif enrichment
findMotifsGenome.pl peaks.bed mm39 motif_output/ \
  -size 200 \
  -mknown /path/to/motif/database/

Output includes: - HTML report with top motifs - Position weight matrices - Known motif matches - De novo discovered motifs

Using MEME Suite

# Extract sequences under peaks
bedtools getfasta -fi genome.fa -bed peaks.bed -fo peak_sequences.fa

module load meme

# Run MEME for de novo motif discovery
meme peak_sequences.fa \
  -dna \
  -oc meme_output \
  -mod zoops \
  -nmotifs 10 \
  -minw 6 -maxw 20 \
  -maxsize 1000000

# OR Run FIMO to scan for known motifs
fimo --oc fimo_output --thresh 1e-3 \
  motif_database.meme \ # Download from MEME Suite or JASPAR
  peak_sequences.fa

Motif Enrichment in R

library(motifmatchr)
library(TFBSTools)
library(JASPAR2020)

# Get JASPAR motifs
opts <- list(species=9606, all_versions=FALSE)
motifs <- getMatrixSet(JASPAR2020, opts)

# Match motifs in peaks
peak_sequences <- getSeq(BSgenome.Mmusculus.UCSC.mm39, peaks)
motif_matches <- matchMotifs(motifs, peak_sequences, genome="mm39")

# Get motif counts
motif_counts <- motifCounts(motif_matches)

Functional Enrichment Analysis

Gene Ontology Enrichment

library(clusterProfiler)
library(org.Mm.eg.db)

# Get genes near peaks
genes <- anno_df$geneId

# GO enrichment
ego <- enrichGO(gene = genes,
                OrgDb = org.Mm.eg.db,
                ont = "BP",
                pAdjustMethod = "BH",
                pvalueCutoff = 0.05,
                qvalueCutoff = 0.2)

# Visualize
dotplot(ego, showCategory=20)
barplot(ego, showCategory=15)
cnetplot(ego, showCategory=5)

Using bedtools

Command-line peak annotation:

# Find nearest gene
bedtools closest -a peaks.bed -b genes.bed -d > peaks_nearest_gene.bed

# Intersect with promoters
bedtools intersect -a peaks.bed -b promoters.bed -wa -wb > peaks_in_promoters.bed

# Count peaks overlapping each gene
bedtools intersect -a genes.bed -b peaks.bed -c > genes_with_peak_counts.bed

Differential Binding Analysis

Using DiffBind (R)

Compare binding between conditions:

library(DiffBind)

# Create sample sheet (CSV file)
# Columns: SampleID, Tissue, Factor, Condition, Replicate, bamReads, Peaks
samples <- read.csv("sample_sheet.csv")

# Create DiffBind object
dba <- dba(sampleSheet=samples)

# Count reads in consensus peaks
dba <- dba.count(dba, bUseSummarizeOverlaps=TRUE)

# Establish contrast
dba <- dba.contrast(dba, categories=DBA_CONDITION, minMembers=2)

# Perform differential analysis (DESeq2 method)
dba <- dba.analyze(dba, method=DBA_DESEQ2)

# Generate report
report <- dba.report(dba, th=0.05)

# Save results
write.csv(as.data.frame(report), "differential_peaks.csv", row.names=FALSE)

Visualize DiffBind Results

# PCA plot
dba.plotPCA(dba, attributes=DBA_CONDITION)

# Correlation heatmap
dba.plotHeatmap(dba)

# MA plot
dba.plotMA(dba)

# Volcano plot
dba.plotVolcano(dba)

# Binding heatmap
dba.plotHeatmap(dba, contrast=1, correlations=FALSE)

Enrichment for Differential Peaks

# Get genes near gained peaks (higher in treatment)
gained_peaks <- subset(report, Fold > 1 & FDR < 0.05)
gained_genes <- unique(gained_peaks$geneId)

# Get genes near lost peaks (lower in treatment)
lost_peaks <- subset(report, Fold < -1 & FDR < 0.05)
lost_genes <- unique(lost_peaks$geneId)

# Compare GO enrichment
ego_gained <- enrichGO(gene=gained_genes, OrgDb=org.Mm.eg.db, ont="BP")
ego_lost <- enrichGO(gene=lost_genes, OrgDb=org.Mm.eg.db, ont="BP")

# Compare
compareCluster(list(Gained=gained_genes, Lost=lost_genes),
               fun="enrichGO",
               OrgDb=org.Mm.eg.db)

Integration with RNA-seq

Correlate Binding with Expression

library(DESeq2)

# Load RNA-seq differential expression results
rna_results <- read.csv("rna_deseq2_results.csv")

# Load ChIP-seq peak annotations
chip_results <- read.csv("annotated_peaks.csv")

# Merge by gene
merged <- merge(chip_results, rna_results, by="gene_id")

# Correlate peak score with expression change
cor.test(merged$peak_score, merged$log2FoldChange)

# Plot
plot(merged$peak_score, merged$log2FoldChange,
     xlab="ChIP Peak Score", ylab="RNA log2FC")

Identify Direct Targets

# Genes with both ChIP peak AND expression change
direct_targets <- subset(merged, 
                        peak_score > 50 & 
                        abs(log2FoldChange) > 1 & 
                        padj < 0.05)

# Export
write.csv(direct_targets, "direct_target_genes.csv")

Chromatin State Analysis

Peak Overlap with Histone Marks

# Overlap ATAC peaks with H3K27ac (active enhancers)
bedtools intersect -a atac_peaks.bed -b H3K27ac_peaks.bed -u > active_enhancers.bed

# Overlap with H3K4me3 (active promoters)
bedtools intersect -a atac_peaks.bed -b H3K4me3_peaks.bed -u > active_promoters.bed

# Overlap with H3K27me3 (repressed regions)
bedtools intersect -a atac_peaks.bed -b H3K27me3_peaks.bed -u > polycomb_regions.bed

Best Practices

  1. Always annotate peaks to understand genomic context
  2. Validate top findings experimentally (ChIP-qPCR, reporter assays)
  3. Integrate with other data (RNA-seq, other ChIP-seq)
  4. Use proper controls for differential analysis
  5. Validate motifs with known biology
  6. Consider batch effects in multi-sample analyses

Additional Resources