DiffBind Analysis

Imports

# if (!require("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")

# Load packages
library(DiffBind)

Sample Sheets

Create sample sheets and comparison contrasts for erythroid ATAC-seq and erythroid RAD21 ChIP-seq

Erythroid ATAC-seq

ery_atac_bam_dir <- "Example_Data/Erythroid_ATAC/BAMs/"
ery_atac_peak_dir <- "Example_Data/Erythroid_ATAC/LanceOtron_Peaks"

ery_atac_bam_dir <- "/ceph/project/Wellcome_Discovery/datashare/towilson/Public_Data/Human/Erythroid/50_Donors_Erythroid/ZEN_T_Wilson_2026/ATAC/hg38/BAMs/"
ery_atac_peak_dir <- "/ceph/project/Wellcome_Discovery/datashare/towilson/Public_Data/Human/Erythroid/50_Donors_Erythroid/ZEN_T_Wilson_2026/ATAC/hg38/LanceOtron_Peaks"

# Set sample name GSM
gsm_ids <- paste0("GSM", 9320314:9320327)

# Sample names
samples <- c(paste0("Don002_ATAC_d13_rep", 1:7),
             paste0("Don003_ATAC_d13_rep", 1:7))

# Erythroid ATAC-seq sample sheet
ery_atac_samples <- data.frame(SampleID = c(paste0("Don002_ATAC_d13_rep", 1:7),
                                            paste0("Don003_ATAC_d13_rep", 1:7)),
                               Tissue = rep("CD34_erythroid", 14),
                               Factor = c(rep("donor_2", 7), 
                                          rep("donor_3", 7)),
                               Condition = rep("healthy", 14),
                               Treatment = rep("bulk", 14),
                               Replicate = c(1:7, 1:7),
                               bamReads = file.path(ery_atac_bam_dir,
                                                    paste0(gsm_ids, "_", samples, ".bam")),
                               Peaks = file.path(ery_atac_peak_dir,
                                                 paste0(gsm_ids, "_", samples, "_L-tron.bed")),
                               PeakCaller = rep("bed", 14))

# Set samples and contrast for donor 2 replicate comparison
ery_don2_samples <- ery_atac_samples[(ery_atac_samples["Factor"] == "donor_2"),][1:6,]
ery_don2_samples["Factor"] = c(rep("group_1", 3), rep("group_2", 3))
ery_don2_contrast <- list(Factor = "group_1")

# Set samples and contrast for donor 2 vs donor 3 replicate comparison
ery_don2_vs_don3_samples <- ery_atac_samples
ery_don2_vs_don3_contrast <- list(Factor = "donor_2")

Erythroid RAD21 ChIP-seq

ery_rad21_bam_dir <- "Example_Data/Erythroid_RAD21/BAMs"
ery_rad21_peak_dir <- "Example_Data/Erythroid_RAD21/LanceOtron_Peaks"

# Erythroid RAD21 ChIP-seq sample sheet
ery_rad21_samples <- data.frame(SampleID = c(paste0("Don001_RAD21_d13_rep", 1:3),
                                             paste0("Don002_RAD21_d13_rep", 1:3),
                                             "Don030_RAD21_d13_rep_1",
                                             "Don030_RAD21_d13_rep_3"),
                                Tissue = rep("CD34_erythroid", 8),
                                Factor = c(rep("donor_1", 3),
                                           rep("donor_2", 3),
                                           rep("donor_30", 2)),
                                Condition = rep("healthy", 8),
                                Treatment = rep("bulk", 8),
                                Replicate = c(1:3, 1:3, 1, 3),
                                bamReads = c(file.path(ery_rad21_bam_dir, paste0("Don001_RAD21_d13_rep", 1:3, ".bam")),
                                             file.path(ery_rad21_bam_dir, paste0("Don002_RAD21_d13_rep", 1:3, ".bam")),
                                             file.path(ery_rad21_bam_dir, "Don030_RAD21_d13_rep1.bam"),
                                             file.path(ery_rad21_bam_dir, "Don030_RAD21_d13_rep3.bam")),
                                Peaks = c(file.path(ery_rad21_peak_dir, paste0("Don001_RAD21_d13_rep", 1:3, ".bed")),
                                          file.path(ery_rad21_peak_dir, paste0("Don002_RAD21_d13_rep", 1:3, ".bed")),
                                          file.path(ery_rad21_peak_dir, paste0("Don030_RAD21_d13_rep1", ".bed")),
                                          file.path(ery_rad21_peak_dir, paste0("Don030_RAD21_d13_rep3", ".bed"))),
                            PeakCaller = rep("bed", 8))

# Set the contrast
contrast <- list(Factor = "donor_1")

Extract LIB, RLE and TMM Scaling Factors

Extract scaling factors to rescale bigWigs by full library size (LIB), RLE and TMM normalisation

# Set sample sheets per dataset
sample_sheets <- list("Erythroid_ATAC" = ery_atac_samples, "Erythroid_RAD21" = ery_rad21_samples)
# Autosomal and sex chromosomes
chromosomes <- c(paste0("chr", 1:22), "chrX", "chrY")
# Normalisation methods to test
norm_methods <- c(DBA_NORM_LIB, DBA_NORM_RLE, DBA_NORM_TMM)

scaling_factors <- list()

for (data in c("Erythroid_ATAC", "Erythroid_RAD21")) {
  # Create DBA object
  dba_object <- dba(sampleSheet = sample_sheets[[data]])
  
  # Find chromosomes to remove
  peakset <- dba.peakset(dba_object, bRetrieve = TRUE)
  greylist <- peakset[!(seqnames(peakset) %in% chromosomes)]
  
  # Apply ENCODE blacklist and keep only standard chromosomes
  dba_object <- dba.blacklist(dba_object, blacklist = DBA_BLACKLIST_HG38, greylist = greylist)
  
  # Re-center peaks and count reads
  dba_object <- dba.count(dba_object)
  
  # Normalise by each method
  for (norm in norm_methods) {
    dba_object <- dba.normalize(dba_object, normalize = norm)
    # Extract the scaling factors
    norm_info <- dba.normalize(dba_object, bRetrieve = TRUE)
    scalers <- norm_info$norm.factors
    names(scalers) <- dba_object$samples$SampleID
    scaling_factors[[data]][[norm]] <- scalers
  }
}

Error in if (file.info(peaks)$size > 0) { : 
  missing value where TRUE/FALSE needed

# View scaling factors per dataset and normalisation method
scaling_factors

Differential Accessibility Analysis

ZEN Scaling Factors

# Scaling factors per sample for ZEN normalisation
ery_atac_zen <- c(Don002_ATAC_d13_rep1 = 0.2795583724430078,
                  Don002_ATAC_d13_rep2 = 0.24549477454191823,
                  Don002_ATAC_d13_rep3 = 0.1651864697590268,
                  Don002_ATAC_d13_rep4 = 0.09541567508135224,
                  Don002_ATAC_d13_rep5 = 0.16076437035969782,
                  Don002_ATAC_d13_rep6 = 0.1551683098985989,
                  Don002_ATAC_d13_rep7 = 0.143536842576444,
                  Don003_ATAC_d13_rep1 = 0.2582363354917962,
                  Don003_ATAC_d13_rep2 = 0.25053836607363844,
                  Don003_ATAC_d13_rep3 = 0.17587857961294,
                  Don003_ATAC_d13_rep4 = 0.11435196335506045,
                  Don003_ATAC_d13_rep5 = 0.10978229659749278,
                  Don003_ATAC_d13_rep6 = 0.1178921511723179,
                  Don003_ATAC_d13_rep7 = 0.0877359143756193)

Run DiffBind

# Create DBA object
sheet <- ery_don2_vs_don3_samples
contrast <- ery_don2_vs_don3_contrast
  
dba_object <- dba(sampleSheet = sheet)
dba_object

14 Samples, 101677 sites in matrix (396696 total):

# Clustering by peak calls
dba.plotHeatmap(dba_object)

# Find chromosomes to remove
chromosomes <- c(paste0("chr", 1:22), "chrX", "chrY")
peakset <- dba.peakset(dba_object, bRetrieve = TRUE)
greylist <- peakset[!(seqnames(peakset) %in% chromosomes)]

# Apply ENCODE blacklist and keep only standard chromosomes
dba_object <- dba.blacklist(dba_object, blacklist = DBA_BLACKLIST_HG38, greylist = greylist)

# Re-center peaks and count reads
dba_object <- dba.count(dba_object)
dba_object

14 Samples, 65625 sites in matrix:

# DBA_NORM_LIB (default), DBA_NORM_RLE (DESeq2), DBA_NORM_TMM (edgeR) or custom vector (ZEN)
norm_method <- "ZEN"
scalars <- ery_atac_zen
  
if (norm_method == "ZEN") {
  # Select scaling factors for the specific replicates
  norm_sub <- scalars[names(scalars) %in% sheet$SampleID]
} else {
  norm_sub <- norm_method
}

dba_object <- dba.normalize(dba_object, normalize = norm_sub)
norm_info <- dba.normalize(dba_object, bRetrieve = TRUE)
# View the scaling factors
norm_info

$norm.method
[1] "user"

$norm.factors
Don002_ATAC_d13_rep1 Don002_ATAC_d13_rep2 Don002_ATAC_d13_rep3 Don002_ATAC_d13_rep4 Don002_ATAC_d13_rep5 Don002_ATAC_d13_rep6 
          0.27955837           0.24549477           0.16518647           0.09541568           0.16076437           0.15516831 
Don002_ATAC_d13_rep7 Don003_ATAC_d13_rep1 Don003_ATAC_d13_rep2 Don003_ATAC_d13_rep3 Don003_ATAC_d13_rep4 Don003_ATAC_d13_rep5 
          0.14353684           0.25823634           0.25053837           0.17587858           0.11435196           0.10978230 
Don003_ATAC_d13_rep6 Don003_ATAC_d13_rep7 
          0.11789215           0.08773591 

$lib.method
[1] "full"

$lib.sizes
 [1] 22105785 32805795 49429595 38990930 16896174 12451638 27350017 22723762 33719954 61605892 19314496 14912872 12046132
[14] 25140080

$filter.value
[1] 1

# Clustering by count scores
plot(dba_object)

# Set the contrast
dba_object <- dba.contrast(dba_object, reorderMeta = contrast)
dba_object <- dba.analyze(dba_object)
dba.show(dba_object, bContrasts = TRUE)

plot(dba_object, contrast = 1)

results <- dba.report(dba_object, contrast = 1)
results

GRanges object with 145 ranges and 6 metadata columns:
        seqnames              ranges strand |      Conc Conc_donor_3 Conc_donor_2      Fold     p-value         FDR
           <Rle>           <IRanges>  <Rle> | <numeric>    <numeric>    <numeric> <numeric>   <numeric>   <numeric>
    590     chr1   25285078-25285478      * |   8.91678      0.00000      9.91678  -8.79572 2.77988e-21 1.82430e-16
  51349     chr6   79614820-79615220      * |   9.26869      5.13925     10.22687  -4.93362 3.44500e-17 7.62786e-13
    588     chr1   25279697-25280097      * |   8.60247      2.21490      9.59383  -7.13503 3.48702e-17 7.62786e-13
   2415     chr1 120850923-120851323      * |   6.98049      7.98049      0.00000   7.24046 5.86071e-16 9.61523e-12
  43780     chr4   68864486-68864886      * |   8.84793      4.97343      9.79789  -4.67314 4.17293e-15 4.98268e-11
    ...      ...                 ...    ... .       ...          ...          ...       ...         ...         ...
  54010     chr7   26556232-26556632      * |   6.65763      7.40864      5.00028  1.194607 8.61785e-05   0.0401097
   2670     chr1 153617552-153617952      * |   6.62390      4.96119      7.37592 -1.123389 8.71989e-05   0.0402988
   8093    chr11   16596208-16596608      * |   6.47146      4.62599      7.25532 -1.230685 9.04683e-05   0.0415174
   3180     chr1 171683840-171684240      * |   6.55215      5.14395      7.25101 -0.875443 1.06984e-04   0.0487557
  16477    chr14   21796913-21797313      * |   6.38132      4.73246      7.13073 -0.986939 1.08884e-04   0.0492796
  -------
  seqinfo: 23 sequences from an unspecified genome; no seqlengths

output_file <- "Erythroid_ATAC/DiffBind_Erythroid_Donor2_vs_Donor3_RLE.csv"

results_df <- as(results, "data.frame")
results_df <- results_df[c("seqnames", "start", "end", "Fold", "p-value", "FDR")]
names(results_df)[names(results_df) == "seqnames"] <- "chrom"
print(paste("Significant:", nrow(results_df[results_df$FDR < 0.05, ])))
write.csv(results_df, output_file, row.names = FALSE)
results_df