Chromatin Immunoprecipitation Sequencing

ChIP-Seq data for three factors, NRSF, CTCF, and FoxA1, were used in this study. ChIP-chip and ChIP-Seq (2.2 million ChIP and 2.8 million control uniquely mapped reads, simplified as 'tags') data for NRSF in Jurkat T cells were obtained from Gene Expression Omnibus (GSM210637) and Johnson et al. [8 (link)], respectively. ChIP-Seq (2.9 million ChIP tags) data for CTCF in CD4⁺ T cells were derived from Barski et al. [5 (link)].
ChIP-chip data for FoxA1 and controls in MCF7 cells were previously published [1 (link)], and their corresponding ChIP-Seq data were generated specifically for this study. Around 3 ng FoxA1 ChIP DNA and 3 ng control DNA were used for library preparation, each consisting of an equimolar mixture of DNA from three independent experiments. Libraries were prepared as described in [8 (link)] using a PCR preamplification step and size selection for DNA fragments between 150 and 400 bp. FoxA1 ChIP and control DNA were each sequenced with two lanes by the Illumina/Solexa 1G Genome Analyzer, and yielded 3.9 million and 5.2 million uniquely mapped tags, respectively.

MACS is implemented in Python and freely available with an open source Artistic License at [16 ]. It runs from the command line and takes the following parameters: -t for treatment file (ChIP tags, this is the ONLY required parameter for MACS) and -c for control file containing mapped tags; --format for input file format in BED or ELAND (output) format (default BED); --name for name of the run (for example, FoxA1, default NA); --gsize for mappable genome size to calculate λ_BG from tag count (default 2.7G bp, approximately the mappable human genome size); --tsize for tag size (default 25); --bw for bandwidth, which is half of the estimated sonication size (default 300); --pvalue for p-value cutoff to call peaks (default 1e-5); --mfold for high-confidence fold-enrichment to find model peaks for MACS modeling (default 32); --diag for generating the table to evaluate sequence saturation (default off).
In addition, the user has the option to shift tags by an arbitrary number (--shiftsize) without the MACS model (--nomodel), to use a global lambda (--nolambda) to call peaks, and to show debugging and warning messages (--verbose). If a user has replicate files for ChIP or control, it is recommended to concatenate all replicates into one input file. The output includes one BED file containing the peak chromosome coordinates, and one xls file containing the genome coordinates, summit, p-value, fold_enrichment and FDR (if control is available) of each peak. For FoxA1 ChIP-Seq in MCF7 cells with 3.9 million and 5.2 million ChIP and control tags, respectively, it takes MACS 15 seconds to model the ChIP-DNA size distribution and less than 3 minutes to detect peaks on a 2 GHz CPU Linux computer with 2 GB of RAM. Figure S6 in Additional data file 1 illustrates the whole process with a flow chart.

All data were obtained from ArrayExpress, unless stated otherwise. Human tissue RNA-seq data were obtained from: OCCAMS consortium (European Genome-Phenome Archive, EGAD00001007496). Human tissue ATAC-seq data were obtained from: E-MTAB-5169 (Britton et al., 2017 (link)), E-MTAB-6751 (Rogerson et al., 2019 (link)), and E-MTAB-8447 (Rogerson et al., 2020 (link)). The Cancer Genome Atlas OAC ATAC-seq data were obtained from the GDC data portal (https://portal.gdc.cancer.gov/; Corces et al., 2018 (link)). OE19 H3K27ac ChIP-seq was obtained from: E-MTAB-10319 (Ogden et al., 2022 (link)). GAC H3K4me1 and H3K4me3 ChIP-seq were obtained from: Gene Expression Omnibus, GSE75898 (Ooi et al., 2016 (link)). OE19 siKLF5 RNA-seq and KLF5 ChIP-seq were obtained from: E-MTAB-8446 and E-MTAB-8568, respectively (Rogerson et al., 2020 (link)). OE19 dnFOS RNA-seq was obtained from E-MTAB-10334 (Ogden et al., 2023 (link)).

CUT&Tag library generation was performed as described previously (Kaya-Okur et al., 2020 (link)) with an altered nuclear extraction step. For the nuclear extraction, OE19 cells were initially lysed in Nuclei EZ lysis buffer (Sigma-Aldrich, NUC-101) at 4°C for 10 min followed by centrifugation at 500 × g for 5 min. The subsequent clean-up was performed in a buffer composed of 10 mM Tris–HCl pH 8.0, 10 mM NaCl and 0.2% NP40 followed by centrifugation at 1300 × g for 5 min. Nuclei were then lightly cross-linked in 0.1% formaldehyde for 2 min followed by quenching with 75 mM glycine followed by centrifugation at 500 × g for 5 min. Cross-linked nuclei were resuspended in 20 mM N-2-hydroxyethylpiperazine-N'-2-ethanesulfonic acid (HEPES) pH 7.5, 150 mM NaCl, and 0.5 M spermidine at a concentration of 4–8 × 10³ / μl (2–4 × 10⁴ total). Subsequent stages were as previously described (Kaya-Okur et al., 2020 (link)). For 2–4 × 10⁴ nuclei, 0.5 μg of primary and secondary antibodies were used with 1 μl of pA-Tn5 (Epicypher, 15-1017). Antibodies used: anti-BRD4 (abcam, ab128874), anti-CTCF (Merck-Millipore, 07-729), anti-H3K27ac (abcam, ab4729), anti-H3K27me3 (Merck-Millipore, 07-449), anti-H3K4me1 (abcam, ab8895), anti-H3K4me2 (Diagenode, pAb-035-010), anti-H3K4me3 (abcam, ab8580), anti-H3K36me3 (Diagenode, pAb-058-010), anti-H4K20me1 (Diagenode, mAb-147-010), anti-PolII (abcam, ab817), anti-PolII-S2 (abcam, ab5095), anti-PolII-S5 (abcam, ab5131), and anti-Med1 (AntibodyOnline, A98044/10 UG). CUT&Tag libraries were pooled and sequenced on an Illumina HiSeq 4000 System (University of Manchester Genomic Technologies Core Facility). CUT&Tag data processing was performed as for ChIP-seq but with the MACS2 v2.1.1 (Zhang et al., 2008 (link)) but the --broad peak calling option was used for the H4K20me1, H3K27me3 and H3K36me3 marks. Fraction reads in peak (FRiP) scores for each mark were calculated using featureCounts and a stringent threshold of ≥2% was set to ensure quality of data for downstream analyses (Landt et al., 2012 (link); FRiP scores are listed in Supplementary file 12).
ChromHMM (Ernst and Kellis, 2012 (link)) was used to train an eight-state HMM using the CUT&Tag data for all marks assayed. The number of states was determined by running the model with increasing numbers of states until state separation was observed. Emission states were annotated in accordance with Roadmap Epigenomics Consortium Data (Kundaje et al., 2015 (link)).

ChIP-seq analysis was carried out as described previously (Wiseman et al., 2015 (link)). OE19 H3K27ac and GAC H3K4me1/3 ChIP-seq reads were mapped to the human genome GRCh38 (hg38) using Bowtie2 v2.3.0 (Langmead and Salzberg, 2012 (link)). Biological replicates were checked for concordance (r > 0.80). Peaks were called using MACS2 v2.1.1, using input DNA as control (Zhang et al., 2008 (link)). Mapped reads (≥q30) were retained using SAMtools (Li et al., 2009 (link)). Reads mapping to blacklisted regions were removed using BEDtools (Quinlan and Hall, 2010 (link)).