Low HER2 enables dedifferentiation and transformation of normal breast epithelial cells via chromatin opening

Overexpression of the human epidermal growth factor 2 (HER2) protein in breast cancer patients is a predictor of poor prognosis and resistance to therapies. Despite significant advances in the development of targeted therapies and improvements in the 5-year survival rate of metastatic HER2-positive breast cancer patients, a better understanding of the disease at an early stage is needed to prevent its progression. Here, we used an inducible breast cancer transformation system that allows investigation of early molecular changes at high temporal resolution. HER2 overexpression to similar levels as those observed in a subtype of HER2 positive breast cancer patients induced transformation of MCF10A cells and resulted in gross morphological changes, increased anchorage-independent growth of cells, and altered transcriptional programme of genes associated with oncogenic transformation. Global phosphoproteomic analysis during the first few hours of HER2 induction predominantly detected an increase in protein phosphorylation. Intriguingly, this correlated with a wave of chromatin opening, as measured by ATAC-seq on acini isolated from 3D cell culture. We observed that HER2 overexpression leads to reprogramming of many distal regulatory regions and promotes reprogramming-associated heterogeneity. We found that a subset of cells acquired a dedifferentiated breast stem-like phenotype, making them likely candidates for malignant transformation. Our data show that this population of cells, which counterintuitively enriches for relatively low HER2 protein abundance and increased chromatin accessibility, possesses transformational drive, resulting in increased anchorage-independent growth in vitro compared to cells not displaying a stem-like phenotype. Our data provide a discovery platform for signalling to chromatin pathways in HER2-driven cancers, offering an opportunity for biomarker discovery and identification of novel drug targets.

10 *Corresponding author 11 12 Overexpression of the human epidermal growth factor 2 (HER2) protein in breast cancer patients is a 13 predictor of poor prognosis and resistance to therapies. Despite significant advances in the 14 development of targeted therapies and improvements in the 5-year survival rate of metastatic 15 HER2-positive breast cancer patients, a better understanding of the disease at an early stage is 16 needed to prevent its progression. Here, we used an inducible breast cancer transformation system 17 that allows investigation of early molecular changes at high temporal resolution. HER2 18 overexpression to similar levels as those observed in a subtype of HER2 positive breast cancer 19 patients induced transformation of MCF10A cells and resulted in gross morphological changes, 20 increased anchorage-independent growth of cells, and altered transcriptional programme of genes 21 associated with oncogenic transformation. Global phosphoproteomic analysis during the first few 22 hours of HER2 induction predominantly detected an increase in protein phosphorylation. 23 Intriguingly, this correlated with a wave of chromatin opening, as measured by ATAC-seq on acini 24 isolated from 3D cell culture. We observed that HER2 overexpression leads to reprogramming of 25 many distal regulatory regions and promotes reprogramming-associated heterogeneity. We found 26 that a subset of cells acquired a dedifferentiated breast stem-like phenotype, making them likely 27 candidates for malignant transformation. Our data show that this population of cells, which 28 counterintuitively enriches for relatively low HER2 protein abundance and increased chromatin 29 accessibility, possesses transformational drive, resulting in increased anchorage-independent growth 30 in vitro compared to cells not displaying a stem-like phenotype. Our data provide a discovery 31 platform for signalling to chromatin pathways in HER2-driven cancers, offering an opportunity for 32 biomarker discovery and identification of novel drug targets. 33 34 Introduction 35 Metastasis is the main cause of cancer deaths but understanding the root cause of malignant 36 transformation remains poorly understood. Many questions remain unanswered as to what triggers 37 cancer formation beyond DNA mutations in pre-cancerous tissue (Ciccarelli and DeGregori, 2020). 38 Perturbed signalling due to dysregulated phosphorylation of oncogenic proteins is known to alter 39 pathway activity and contributes to cellular transformation (Sever and  To recapitulate the early transformational events and the stochastic nature of early breast cancer 76 development, we generated a controllable in vitro model system by stably transducing a 77 doxycycline-inducible HER2 construct in MCF10A cells (Carter et al., 2017). This model allows for the 78 generation of transformed phenotypes in a synchronised and time-controlled manner and is useful 79 for investigating early transformational events using multi-omic analysis (Fig 1A). To analyse the 80 range of HER2 expression at the protein level, we cultured cells for 24 hours in five different 81 concentrations of doxycycline, using ranges that have been used previously in inducible expression 82 studies with other proteins (Baron et al., 1995;Leitner et al., 2014). In our model, a 24-hour 83 induction with 1µg/ml doxycycline resulted in strong HER2 protein expression (Fig 1B). When grown 84 in three-dimensional cell cultures, control MCF10A cells (MCF10A CTRL ) formed regular, spherical acini, 85 whereas a majority of MCF10A HER2 acini were misshapen, with cells budding into the surrounding 86 matrix (Fig 1C). HER2 overexpression resulted in significantly increased in vitro migratory and 87 invasive potential, as measured by transwell assays (Fig 1D) (Xiang and Muthuswamy, 2006;Paszek 88 and Weaver, 2004). Furthermore, MCF10A HER2 cells displayed a hallmark of in vitro transformation, 89 with increased anchorage-independent growth as compared to control cells (Fig 1E).  (D) Cell migration and invasion was analysed through the 8µm pores of transwell membranes over 16-hour period of chemotactic migration towards full serum media. The ability of cell invasion was measured in collagen or matrigel coated transwells. Migration ability was measured in using uncoated wells. Statistical significance was calculated using student's t-test. Significance is shown as * for p-value < 0.05, ** for p-value < 0.01. N=3.
(E) Colony growth of MCF10A HER2 and control cells in 0.3% ultra-pure agarose over 3 weeks. ImageJ analysis of 5 different size colonies were quantified. Representative microscopic images of colonies stained with crystal violet after three weeks. Statistical significance was calculated using student's t-test). Significance is shown as * for p-value < 0.05, ** for p-value < 0.01, *** for p-value < 0.001.
Images are at 1.6x magnification. Scale bars represent 1000μm. N=3. 8 quite low. Although HER2 protein expression is still low, some of these downstream changes might 132 be present at later timepoints as part of the evolution process. 133 The low levels of HER2 activation at early time points may closely mimic, at least partially, the early 134 signalling changes occurring in HER2 positive breast cancer patients. The signalling changes of low 135 level HER2 induction has not been performed to date. We re-analysed this data by decreasing the 136 significance threshold to log2fold change > 0.5, FDR corrected p-value of < 0.05 for HER2 expression, points analysed, we identified that mitogen-activated protein kinase (MAPK) signalling pathway to 144 be one of the most enriched cascades in our system (Supplementary Fig 1D). The idea that signalling 145 has direct effects on chromatin has already been known, whereby receptor tyrosine kinases can identified in the phosphoproteomic screen are likely to contribute to chromatin changes mediating 153 the transformed phenotypes. Indeed, our phosphoproteomic analysis revealed significant changes in 154 various transcription factors known to affect chromatin dynamics (Supplementary Fig 1E). These In particular, the phosphorylation of JUN at residue S73 could be reconciled by a model in which 157 phosphorylation of JUN triggers dissociation of histone deacetylases (HDACs) and facilitates the 158 rearrangement of chromatin structure (Wolter et al., 2008). Based on these results, we then set out 159 to assess, in an unbiased manner, the effects that signalling changes have on the chromatin 160 organisation. 161 (A) Volcano plots depicting the phosphoproteome upon HER2 protein expression at 0.5-hour, 4 hours, and 7-hour time-points compared to control cells. The statistical significance is shown as (log2 fold change for HER2 > 1.5, p-value < 0.05, and log2 fold change for GFP < 5, p value of > 0.05). The plot contains those phosphosites that are significantly changing upon HER2 protein induction but not significantly changing in the GFP cells at the same time. Those with the highest increase or decrease in fold change are labelled. N=3.
(B) Bar graph depicting the number of detected phosphosites increasing or decreasing in phosphorylation in the phosphoproteomic analysis at the indicated time-points analysed. The statistical significance is shown as (log2 fold change for HER2 > 1.5, p-value < 0.05, and log2 fold change for GFP < 5, p value of > 0.05).
(C) Differential accessibility (log2 fold change > 0.5, FDR corrected p-value of < 0.05) between MCF10A HER2 and control cells, plotted against the mean reads per region. Cells were grown in 3D cell culture from 0-48 hours and ATAC-seq performed on their acini.
Heatmap shows chromatin accessibility across all time points for each replicate in cells expressing HER2 or GFP (controls). N=3.
(D) Fraction of total regions that are differentially accessible (up peaks) or inaccessible (down peaks) in early or late type comparisons.
Log2fold > 2, FDR corrected p-value < 0.05. separated the samples into two groups, "early" (0h, 1h, 4h, and 7h time-points) and "late" (24h and 172 48h time-points) (Supplementary Fig 1F). We selected these conditions with the aim to encompass 173 time-points relevant to both types of analysis. The 0h, 4h and 7h time-points were chosen to 174 characterise early chromatin changes triggered by signalling. The late conditions were selected to 175 detect the resulting delayed chromatin changes occurring later in the process of transformation. We 176 identified 17,868 significant changes between MCF10A HER2 cells relative to control cells (T0 starting 177 population before HER2 protein induction) over the time course, which showed an increase in 178 accessibility in MCF10A HER2 cells relative to controls (Fig 2C & supplementary Fig 2A). We assessed 179 differential accessibility between early and late groups and observed that a much larger fraction of 180 regions, with > 2-fold difference relative to T0, were enriched in the early group compared to in the 181 late group (75% vs 44%, respectively, Fig 2D). Conversely, only ~2.9% of peaks were >4-fold more 182 accessible in the early group and ~6.5% in the late group, which we define as "hyper-accessible" 183 chromatin states ( Supplementary Fig 2A). Even though the numbers of hyper-accessible versus 184 hypo-accessible regions (which lose accessibility > 4-fold) did not show a stark difference, the overall 185 number of accessible regions following HER2 expression outnumbered inaccessible regions. This 186 shows that there is an increase in chromatin accessibility during the early stages of transformation 187 ( Fig 2D). Therefore, this might suggest that the first adaptive response to oncogenic HER2 signalling 188 is altered chromatin accessibility to induce differential gene expression. Subsequently, the changes 189 in chromatin accessibility even out in the later time points, with the number of hypo-accessible regions even exceeding the hyper-accessible ones at late time points, which could indicate that cells 191 have reached an equilibrium. (Supplementary Fig 2A). 192 Next, we performed functional enrichment analyses [Gene Ontology (GO) terms] for upregulated 193 peaks in the early HER2 signature (Fig 2E). We next examined whether peaks were shared between those that were opening (more accessible) 212 and those that were closing (less accessible) between the early and late groups. We found that there 213 was a small overlap between early and late inaccessible peaks but none between the accessible 214 sites with early loss in accessibility relative to T0 could potentially have driving roles in the 216 population drift. We further examined the genomic distribution of the differentially inaccessible 217 chromatin of the overlapping regions, which showed most genomic regions were associated with 218 two nearby genes (Supplementary Fig 2D). Namely, some of the common differential regions 219 correlated with genomic location of FBN2, whose genomic chromosomal coordinates were found to 220 be matching with the promoter region of the FBN2 gene. This gene was found to have aberrant 221 promoter methylation in a number of cancers (Hibi et al., 2012) (Supplementary Fig 2E). Other 222 regions included RIMS2, known to be associated with particularly aggressive breast cancers (Zhang,223 L., Liu and Zhu, 2021) and APIP, which binds HER3 receptor, leading to the heterodimerisation 224 between HER2-HER3 and resulting in sustained activation of downstream signalling (Hong et al., 225 2016). No differentially accessible region was found to be promoter proximal, as all the regions were 226 at least 5 kb upstream of the transcriptional start site (TSS) (Supplementary Fig 2F).  Fig 3A). Seurat 232 clustering found differentially expressed features and separated them into four groups, with cluster 233 0 enriching in the MCF10A CTRL population, wand cluster 1 associating with the highest HER2 234 expression (Supplementary Fig 3B). As expected, we observed a time-dependent increase in HER2 235 gene expression (Supplementary Fig 3C). There is a consensus that high HER2 expression is 236 associated with stem-like phenotype abundance of cancer stem cells) (Fig 3D and Supplementary Fig 3D); LAMB3, which mediates 243 invasive and proliferative behaviours by PI3K-AKT signalling pathway (Zhang, H. et al., 2019), as well 244 as the decrease of genes like MUC1 conversely upon HER2 overexpression, whose downregulation is 245 linked to stem-like phenotype (Stingl, 2009a). While the expression of ID3 is also associated with 246 stemness (Huang et al., 2019) this pattern was not found in our data, suggesting that these 247 processes overlap only partially.  (A) Distance to closest transcriptional start sites (TSSs) of all differentially accessible regions in the early and late cell types. The graphs represents only those regions that are upstream of the TSS. "Open sea" refers to regions that are at least 50 kb or more upstream of the TSS.
(B) Enrichment of transcription factor recognition sequences in differential ATAC-seq peaks comparing MCF10A HER2 and control cells based on HOMER analysis. Down peaks = log2fold < -2, FDR corrected p-value value < 0.05. HER2 cells were able to form cell aggregates, with a > 2-fold increase in colony forming units 264 compared to control cells (Fig 1E). We hypothesised that cells possessing the ability to form colonies 265 under anchorage-independent growth conditions are a selection of aggressive cells out of the total 266 number of cells seeded. Conversely, the proliferative but non-malignant cells that often dominate 267 any heterogeneous parental cell line would be selected against under these conditions. We 268 evaluated whether anchorage-independent growth correlated with reprogramming-associated 269 heterogeneity by testing the expression of proteins found in mammary epithelial stem cell hierarchy 270 by flow cytometry (Stingl, 2009b), in which it has been shown that breast stem cells are 271 characterised by MUC1 -ve , EpCAM low , and CD24 low expression (Fig 4A). We therefore evaluated 272 whether HER2 overexpression could enrich for cells with functional stem-like properties based on 273 these three markers and found that this stem-like phenotype is enriched in MCF10A HER2 cells, as a 274 large proportion of cells lost the expression of MUC1, EpCAM, and CD24 (Fig 4B, Supplementary Fig  275   3E, and supplementary Fig 4). Since our population is heterogeneous due to differing number of 276 copies of the lentiviral HER2 construct, and we have the same amount of doxycycline used to induce 277 the oncogene, the upper threshold of expression of HER2 will depend on the transgene copy 278 number. We therefore hypothesised that stem-like markers would be positively correlated with 279 HER2 levels in our heterogeneous population, i.e., cells having many HER2 copies would also be 280 more likely to express stem-like markers. Surprisingly, we found that cells expressing relatively low 281 HER2 levels had the most pronounced stem-like phenotype compared to other flow sorted 282 populations of cells with increasing levels of HER2 (Fig 4B). We confirmed the different levels of 283 HER2 protein expression after sorting cells into three compartments of low, medium, and high HER2 transformational potential of these cell types by measuring anchorage-independent growth, we flow 286 sorted MCF10A HER2 cells into the three different cell populations and paradoxically found that low 287 HER2 expressing cells had increased transformational potential relative to the other populations of 288 sorted cells (Fig 4D). We thought that high HER2 cells may be undergoing oncogene-induced 289 senescence (OIS), thus resulting in reduced colony formation compared to other cell types. To 290 confirm this, we measured proteins implicated in OIS but found no significant increase in OIS 291 markers in the high HER2 cells compared to other populations, indicating other biological effects 292 being responsible for the lower capacity in anchorage-independent growth of high HER2 expressing 293 cells (Supplementary Fig 3F). It is possible that high oncogene expression induces cancer cells to 294 dormancy that is associated with loss of ability to self-replicate and differentiate (Bellovin,Das and 295 Felsher, 2013). 296 Since we found that chromatin opening was the feature associated with early signalling to chromatin 297 response, we wanted to know if this was reflected in the phenotypic heterogeneity, in particular low 298 versus high HER2 levels. To this end, we used ATAC-seq to determine the genome-wide chromatin 299 accessibility landscape in the five different populations of cells (MCF10A CTRL , low HER2, medium 300 HER2, high HER2 and MCF10A HER2 cells). We analysed these data by comparing each cell type to the 301 control cells (MCF10A CTRL ) and comparing the percentage of differentially accessible regions 302 between the cell types. We found that low HER2 expressing cells exhibited the highest percentage of 303 chromatin opening compared to other cell populations (Fig 4E), confirming that the phenotypes 304 associated with invasiveness and anchorage independent growth were driven by molecular features 305 in stem-like cells and opening of chromatin. To put the magnitude of these chromatin differences in 306 context, i.e., the differential accessibility between low HER2 and high HER2 expressing cells, we 307 found that a dramatic ~95% of peaks were accessible in low HER2 cells. Conversely, only ~42% of the 308 peaks were open (accessible) in the high HER2 cells. Overall, these data indicate that a sharp 309 increase in HER2, which may result in triggering cell intrinsic defensive systems, whereas a low-level sustained presence of HER2 can shift cell identity, via chromatin remodelling, towards tumour-311 promoting phenotypes. 312 We found that a subset of these scRNAseq-unique differentially expressed genes that were either 313 upregulated or downregulated in multiple time-points were also associated with heterogeneity of 314 breast cancer, related to cancer progression and stem cells (Fig 4F). For example, expression of 315 HMGA1, which is known to promote breast cancer angiogenesis through the transcriptional activity 316 of FOXM1 (Zanin et al., 2019), increased in a time dependent manner (Fig 4F). On the other hand, 317 expression of FOS, a pro-proliferative transcription factor, which has been validated in breast 318 tumour samples and is highly expressed in relapse samples and treatment failures (Vendrell et al., 319 2008) was found to be downregulated at all time points (Fig 4F). Intriguingly, high proliferation rates 320 as a result of FOS expression can lead to improved outcomes for patients with breast cancer, as it 321 could lead to higher apoptosis-effector gene expression (Fisler et   (D) HER2 expression was induced for 3 days, and cells were sorted based on HER2 expression into low, medium, and high HER2 expression. 5000 cells from each condition were plated into ultra-pure agarose to investigate there in vitro transformative potential.
Results are plotted as box plots from three biological replicates. Student t-test was performed to compare "HER2 med" and "HER2 high" groups to the "HER2 low" group, p-values are displayed on the graph. One-way anova test was performed to show statistical significance. N=3. independent growth accompanied by the formation of spindle-like conformations in 3D cell culture 338 (Fig 1B and D).  (Stingl, 2009). Our 359 data suggest an alternative model in which dedifferentiation is paradoxically driven by low levels of 360 HER2 protein expression and creates a programme of stem cell marker expression that drives 361 transformational ability (Fig 4B). 362 We observed leucine aminopeptidase 3 (LAP3) to be significantly activated in our phosphoproteomic 363 screen in all of the time-points analysed in MCF10A HER2 compared to MCF10A GFP cells (Fig 2A).

364
LAP3 is known to play a critical role in breast cancer cells by regulating migration, invasion and is 365 associated with metastasis (Fang et al., 2019). In addition, we found that phosphorylation of 366 nucleolar and coiled-body phosphoprotein 1 (NCOL1) at residue S622 was also significantly increased 367 all time-points (Fig 2A). This protein is found to be highly expressed in nasopharyngeal carcinomas 368 (NPC) (Hwang et al., 2009) and in breast cancer cells (Sacco et al., 2016). The consistent and highly 369 stable activation of these two proteins may serve as potential biomarkers for late-stage disease and 370 provide important targets for antimetastatic therapeutic targets. Furthermore, zinc finger protein 371 (ZFP36) is correlated with lower tumour grade breast cancer. Interestingly, we find that ZFP36 (S188) 372 is significantly activated in the 4-and 7-hour time-points but not in the earlier 30-minute time-points 373 (Fig 2A), indicating that low HER2 expressing cells prefer a programme of signalling phosphosites 374 associated with worse patient outcome (Canzoneri et al., 2020). 375 The morphological changes in breast cancer models are often used to indicate the high 376 transformational characteristics of those cells (Petsalaki E. et al., 2021). We found that proteins 377 associated with aggressive basal-like phenotype were found to be increased in our 378 phosphoproteomic screen, which included ADGRA2 (S1079) and DENND4C (S1250). This shows that 379 the morphological changes observed in our system (Fig 1B) were likely due to HER2 induced 380 transformation. 381 downregulation of IFITM family members is associated with resistance maintenance following anti-384 HER2 therapy, trastuzumab . Our single cell RNA-seq data reveal downregulation 385 of IFITM3 within 24 hours of HER2 overexpression, that is maintained until at least 72 hours, which 386 could show that this does not decrease as a result of resistance but may predispose resistance to 387 therapies at the very early stage of disease. Overall, our data show the power of combining genome-388 wide molecular approaches using an in vitro transformation model system to uncover subtle but ATAC-seq library preparation and differential analysis 446 5x10 5 cells were directly recovered from cell culture by trypsin from 2D cell culture or by using the 447 recovery solution (CORNING #354253) for cells grown in either 2D or 3D cell culture. ATAC-seq 448 libraries were generated as described previously (Buenrostro et al., 2015), with minor amendments.
We performed 10 initial PCR amplification cycles followed by direct purification of the transposed 450 DNA, without performing qPCR to calculate the additional numbers of required cycles. Sequencing 451 data was aligned to the human genome (grch38) using bowtie2. Peaks were called on each biological 452 replicate of all ATAC-seq reads using MACS2, and putative copy number and mitochondrial regions 453 were removed. Peak dataset for differential analysis was generated by applying a threshold using a 454 The intensity values were calculated by determining the peak of each individual extracted ion 471 chromatograms (XIC) and plotted as heatmaps. The resulting quantitative data were transferred and 472 visualised in Microsoft Excel. The significance (log2 fold change < -0.5-fold, FDR corrected p-value of 473 0.05 for the upregulated phosphosites) of each phosphosite was annotated by an asterisk; we used 475 the "filter" function on excel to filter out those phosphosites that were not significant. All of the 476 significant MCF10A GFP data was filtered out, whilst simultaneously filtering out non-significant data 477 for the MCF10A HER2 cells, giving us significant changes in MCF10A HER2 cells that are not significantly 478 changing in the MCF10A GFP cells. The number of phosphosites was determined by the number of 479 columns as each column contains one phosphosite, unless overlapping sites were present, in which 480 case they were manually counted. 481

Migration/Invasion assays 482
Chilled matrigel or collagen mixture was directly pipetted on the centre of 8 μm pore size transwell 483 inserts (MILLICELL #MCEP12H48) that were placed into a 12-well plate, and allowed to solidify at 484 37°C. Meanwhile, cells were trypsinised and pipetted onto the transwell inserts -which were either 485 coated with matrix or left uncoated -and cultured for 16 hours. Highly migratory/invasive cells were 486 stained with 0.05% of crystal violet dye. Images of random regions were taken using a standard light 487 microscope and quantified using imageJ. 488 Soft agar colony formation assays 489 A 0.8% base layer was formed in plates using ultra-pure culture grade agarose (THERMOFISHER 490 #16500500) allowed to settle at room temperature. 5000 cells per well were mixed with 0.3% 491 agarose and plated evenly, drop-wise, on top of the base layer. Media was changed every 2 days for 492 3 weeks. Colonies were fixed using 4% paraformaldehyde (PFA) and permeabilised using 100% 493 methanol. Colonies were stained using 0.05% crystal violet dye and images were taken using a 494 dissecting microscope. Binary masks were applied to each of the images and thresholding 495 parameters for the diameter ranging from 10um 100μm were set on ImageJ. Colonies were counted 496 using ImageJ only if they satisfied criteria above the threshold values, and colony counts were then 497 manually checked and adjusted if necessary. Single cell data was run through the 10X cellRanger pipeline to produce counts tables for gene 513 expression counts, HTO counts for sample identification and ADT counts for HER2 expression. Cells 514 were identified and assigned to a timepoint using the HTO counts table and the HTODemux method 515 in Seurat. To exclude cells that did not respond to the doxycycline induction, treated cells with less 516 than 35 counts of the ERBB2/HER2 expression tag were filtered out. The remaining gene expression 517 data was run through Seurat's basic data processing pipeline. The data was normalized, scaled, and 518 the effects of cell cycle were regressed out using Seurat's cell cycle regression strategy. The data was 519 then run through principal component analysis (PCA). The principal components were used to 520 identify clusters and UMAP was run for visualization. Two different differential expression analyses 521 were run using Seurat's FindAllMarkers function, one across the different clusters and one across the 522 different timepoints. For viewing samples on genome bowser or assessing reproducibility and data exploration, all 560 samples were "down sampled" to the same number of reads; setting: samtools view -b -s 561 [downsampling_ratio] in.bam > out.downsampled.bam. Peaks calling was done for each individual 562 non-downsampled file with MACS2 "callpeak"; setting: MACS2 callpeak -t inbamfile -f BAMPE -n 563 in.bamfile -g ce -keep-dup all. These files were then analysed using DiffBind for differential analysis 564 on R. For each sample, a path to the peaks and the bam file were listed in Microsoft Excel and loaded 565 in R; setting: db.object = dba(sampleSheet="name_of_sample_sheet"). Then, the next step was to 566 take the alignment files and compute count information for each of the peaks/regions in the 567 consensus set; setting: db.object = dba.contrast(db.object, categories=DBA_TREATMENT, 568 block=DBA_CONDITION, minMembers = 2); setting: db.object = 569 dba.plotMA(db.object,th="0.05",method=DBA_DESEQ2). Significant changes could then be saved 572 from up or down peaks e.g.; setting: up_peaks_db.object.SigChanges.0.05FDR <-573 db.object.SigChanges."0.05FDR"[db.object.SigChanges.0.05FDR$Fold > 0,] and counted using the 574 command line and can be plotted as percentage in Prism or Microsoft excel in the form of a 575 chart/graph. 576