dataesr/scientific-paragraphs-categorization huggingface dataset:
Annotation of potential enhancers The active ChromHMM states 1, 8, 9, and 11 comprised 238,820 regions (368 Mb) and represent loci that potentially contain enhancers. This set overlapped with the corresponding B-and T-cell chromatin states of the Roadmap consortium to 89. 8 and 72. 2%, respectively  (Fig EV4A, Dataset EV9). The majority of previously reported enhancer loci from ENCODE and FANTOM projects were found in our state "Active 2" (state 9, H3K4me1, H3K27ac) and with a lower enrichment in states 1, 8,10, and 11 with the bivalent chromatin state 7 being the most enriched repressive state (Fig  EV4B  and C ). Furthermore, binding sites of MEF2, IRF4, and FOXM1 identified by ENCODE were enriched in these states  (Fig EV4B). The most prominent chromatin state transitions from NBCs to CLL cells were from "Quiescent" (state 12) to "Poised" (state 8, H3K4me1), from "Transcribed" (state 2, H3K36me3 in the body of transcribed genes) to "Active 3" (state 1, H3K4me1, H3K27ac, and H3K36me3, predictive for active intragenic enhancers) and from "Poised" to "Active 2" (state 9, H3K4me1 and H3K27ac, predictive of intergenic enhancers; Fig  4A  and B ).
Genome-wide, the highest overlap with bidirectional expression of RNA was observed with "Active 1" (state 10, predictive for TSS) and "Active 2" (state 9, predictive for active intergenic enhancer; Fig  4A  and C ). Bidirectional RNA expression at enhancers could also discriminate between CLL and NBCs, identifying a set of disease-specific active enhancers  (Fig 4D). Next, super-enhancers were annotated by clusters of open chromatin regions co-occurring with large regions of H3K27ac domains in the range of 10 kb and above by using the ROSE software tool  (Whyte et al, 2013; Dataset EV9). In total, 310 out of a union set of 1,289 SEs showed differential activity (false discovery rate (FDR) < 0. 01). There were 219 up-and 91 downregulated SEs in CLL as compared to  NBCs (Fig 4E). Exemplary predicted SEs that became active and upregulated the associated genes in CLL are shown for TCF4, CREB3L2  ,  and FMOD (Figs 1A,  and EV4D  and E ).
NBC samples were in the left cluster with high NFKB2 expression, which was reduced in the right cluster containing the CLL samples. Bottom: A switch of interactions between the NFKB2 promoter and two different enhancers in CLL (red line) vs. NBCs (gray line) was observed. Based on the CTCF ChIP-seq analysis, both intronic enhancers at the NFKB2 and the FBXL15 gene show constitutively bound CTCF nearby, even though the targets of the two enhancers were switched. For the further analysis, we compiled lists of "active enhancers in CLL" and "active enhancers in NBCs". They were derived from the consensus of individual ATAC-seq peak regions flanked by active chromatin states 1 and 9 (Dataset EV9). The two lists comprised 10,145 loci in CLL and 7,312 in NBCs of which 4,771 were shared between the two groups after merging.
Mapping promoter-enhancer interactions from single-cell ATAC-seq data Next, we identified accessible chromatin loci in single cells via single-cell ATAC-seq (scATAC-seq), which can resolve the heterogeneity between individual cells (Appendix Fig  S6A-C ). We analyzed 494 malignant cells from three CLL patients and 343 NBCs from two healthy donors. Distinct signatures of co-occurring open loci were found across different cells with a high overall overlap of open chromatin regions identified by bulk ATAC-seq and scATACseq. The heterogeneity of TF binding between single cells was assessed from the integration frequency of ATAC at known TF binding motifs  (Appendix Fig S6B). Several motifs like NRF1, NFYA/B, and JUN/FOS (AP-1) identified above from the analysis of deregulated CLL chromatin features at promoters and enhancers displayed also the highest heterogeneous accessibility between individual cells in CLL.
Single-cell ATAC-seq can also reveal co-regulation of genomic loci by correlating concomitant accessibility of loci in heterogeneous cell populations. Based on the scATAC-seq data, we computed correlation matrices that contain information about pairs of genomic regions that are simultaneously open in the same cell  (Fig 5F). Correlations between two regions can arise from spatial contacts between them or from other types of co-regulation, e. g. , binding of a common regulator to both loci in the same cell. Using the annotated "active enhancers in CLL/NBC" lists (Dataset EV9), about 800 accessible enhancers per cell were detected on an average. With a threshold derived from the correlation coefficients obtained in a permuted data set (Appendix Fig  S6D  and F ), we then identified enhancer-promoter pairs within 100-kb windows (Dataset EV10).
We found that 68% of these pairs were also listed as spatial contacts in the 4D genome database  (Teng et al, 2015) , suggesting that 2/3 of the promoter-enhancer pairs involve physical contacts. In total, 3,955 promoter-enhancer pairs were identified, with most promoters being connected to one enhancer at mean and median distances of 32 and 20 kb (CLL) and 23 and 10 kb (NBCs), respectively  (Appendix Fig S6E and G). A total of 1,612 pairs were specific for CLL, 1,932 pairs were specific for NBCs, and 411 pairs were present in both CLL and NBCs. In addition, a set of 205 rewired promoters was defined as promoters that do not share any enhancer between CLL and NBCs but were correlated with at least one enhancer in each of the cell types. Within this set of rewired promoters, 70 and 49 were linked to genes down-and upregulated in CLL, respectively.
Assignment of TFs to aberrant CLL chromatin features Aberrant regulatory epigenetic signals were detected at 81% of the transcriptionally deregulated genes in CLL  (Fig 6A, Datasets EV11 and EV12). Based on our TF motif analysis at promoters and enhancers, we selected the most relevant TFs that displayed differential binding with respect to two or more of our chromatin readouts as described in Materials and Methods (Appendix Table  S2 ). These TFs were assigned to four different deregulated CLL chromatin features at promoters or enhancers  (Fig 6B). A gene set enrichment analysis of targets of these TFs retrieved pathways that are highly relevant for the CLL pathophysiology (BCR,  and MAPK signaling; Fig 6C, Appendix Fig S7). Relevant links in this context were, for example, from SYK (Benkisser-Petersen et al, 2016) to NFATC1, BCL2  (Roberts et al, 2016)  to TCF4, and NOTCH1  (Fabbri et al, 2017; Ryan et al, 2017)  to NF-jB.
Integrative gene regulatory enhancer network analysis The pair-wise correlations between enhancers and promoters computed from the scATAC-seq data was combined with the B-cell GRN to derive a "gene regulatory enhancer containing network" (GREN; Datasets EV13 and EV14). From the complete GREN, a CLLspecific GREN was extracted that includes the connected network derived from the CLL-specific TF list, their target genes, as well as linked chromatin modifiers that affect the aberrant chromatin features in CLL. As an additional requirement for this selected network part, all included factors needed to be deregulated in their activity/expression between CLL and NBCs. A part of this network that is centered around TCF4 and EBF1 is shown in Interestingly, genes that in CLL lost EBF1 binding at their enhancers could be clustered into the functional groups of immune response and cell activation, two pathways of central relevance in the pathophysiology of CLL  (Appendix Fig S8E).
Interestingly, boundaries of PMDs were enriched for occupied CTCF sites, which might indicate a role of CTCF to limit the expansion of PMDs. This view is supported by the finding that CTCF binding sites can act as bifurcation points for differential DNA methylation and that CTCF binding can be modulated by DNA demethylation activity  (Teif et al, 2014; Wiehle et al, 2019). Our integrative analysis of a large set of readouts revealed altered chromatin features at promoter/enhancer elements for 81% of the differentially regulated genes in CLL cells (Fig  6A ). The underlying CLL-specific TF network that we derived from aberrant chromatin features displayed strong connectivity with epigenetic readers and writers such as SIN3 and the NuRD and SWI/SNF chromatin remodeling complexes that were deregulated in CLL  (Appendix Fig S7, Dataset EV14).
Furthermore, we find a number of additional links in our core TF network to modifiers of H3K4me3 that included, for example, KMT2E, KDM5A, SETD7  (Appendix Fig S7, Dataset EV14). Extensive changes of enhancer activity were observed that allowed us to clearly discriminate between CLL and NBCs, with a pronounced loss of enhancer activity in CLL. It could be related to an elevated HDAC activity measured in CLL samples and was reverted in part by panobinostat treatment that led to increased H3K27ac at enhancers. Our scATAC-seq analysis revealed correlations between bona fide active sites and led us to propose that enhancer-mediated changes of transcription activity occurred via rewiring to different target promoters in CLL cells. A similar data set does not exist for CLL, as only direct physical contact have been mapped at low resolution by Hi-C  (Beekman et al, 2018)  or for selected interactions by 4C  (Ott et al, 2018).
From our chromatin feature maps, we derived a highly connected CLL-specific network centered around the TFs targeting 17 central binding motifs and an enrichment of BCR signaling genes  (Fig 6B and C, Appendix Fig S7 and Table S2, Dataset EV14). These central motifs include gained motifs in CLL for NFAT, TCF4, and LEF1 and lost motifs for EBF1 and AP-1, which have similarly been reported in two other studies  (Oakes et al, 2016; Beekman et al, 2018). The integrated view on the interplay of TFs, chromatin modifiers, and associated target genes derived here provides a rich resource to generate hypotheses for novel molecular links to the CLL pathophenotype. For example, the regulation of BCL2 via TCF4 predicted in our GREN (Appendix Fig  S7 ) is of special interest since it is the target of venetoclax, a highly effective treatment option in CLL  (Roberts et al, 2016).
Normal B-cell programming might account for some of the methylation changes from the comparison of CLL to NBCs  (Oakes et al, 2016). To account for this, DMRs were filtered out if they show "similar" methylation changes (< 0. 2 b-value) also in the comparison of high-maturity memory B cells to naive B cells as described previously  (Oakes et al, 2016). Methylation profiles were employed to assign CLL samples into three distinct CLL subtypes (low, intermediate, and high) programmed using a reference 450K methylation data from 185 CLL cases. Raw 450K data were processed using the Minfi Bioconductor package (v1. 16. 1; Aryee et al, 2014). Phylogenetic analysis was carried out as previously reported  (Brocks et al, 2014; Oakes et al, 2016)  on the 450K data of CLL cases and normal B-cell subsets together with the WGBS data of CLL cases from our study.
ChIP-seq data analysis After mapping with bowtie to the hg19 genome assembly with decoy regions (options -best -strata -v 1 -m 1), peaks were called with MACS for H3K4me1, H3K4me3, H3K9ac, and H3K27ac using H3 as a background control. For H3K9me3, H3K27me3, and H3K36me3, SICER was used to determine enriched regions with histone H3 serving again as a control. In order to computationally validate the quality of the ChIP-seq experiments, we clustered the different IP samples using their Pearson correlation. Additionally, quality scores suggested by the ENCODE project were used to assess the overall quality of each ChIP-seq data set (Dataset EV4). Peaks were then used for chromatin state analysis with ChromHMM as stated below. Topology and size distributions of modified chromatin regions were calculated with MCORE as described in  Molitor et al (2017)  with the modification that histone H3 was used as a background file for normalization.
CTCF and EBF1 ChIP-seq sequencing libraries were mapped as described above, and peaks were called with MACS2 using standard parameters for narrow peak detection. Differential histone modification and TF binding analysis from ChIP-seq data were performed using the peaks called for each sample and the corresponding reads counts within the alignment files using DiffBind (Dataset EV8). Read counts were scaled to the control, and normalization was done relative to the number of reads within peak regions. Differential enhancers were called using edgeR with an FDR threshold of 0. 01 within DiffBind. For the 24-h mock treatment and panobinostat samples, paired analysis was performed using DESeq2 with an FDR threshold of 0. 01 within DiffBind. The panobinostat differential H3K27ac regions were annotated as promoter overlapping with respect to any TSS from Gencode v17 and annotated as enhancer overlapping with respect to the enhancers of combined CLL occurring in at least three untreated or two 24-h-treated samples.
Genome segmentation and chromatin annotation The chromatin was segmented and annotated using a multivariate Hidden Markov model with the ChromHMM software  (Ernst & Kellis, 2012). The model was trained using seven histone modification marks (H3K4me1, H3K4me3, H3K27ac, H3K27me3, H3K9ac, H3K9me3, and H3K36me3) across all samples where histone modification ChIP-seq data were of sufficient quality for all histone marks (samples H1, H3-5, H8, CLL1-3, CLL5-11). The model was learned by using the merged peak called with MACS or SICER from the ChIP-seq data as binarized input for ChromHMM and allowing for a maximum of 400 iterations. We generated and compared models with 7-50 states and decided to use 12 states (Datasets EV6 and EV7). This model captured all combinations of biologically interpretable states with median correlation of 0.
ATAC-seq data analysis Both scATAC-seq and bulk ATAC-seq were processed in a similar manner. Reads were trimmed using Trimmomatic version 0. 36  (Bolger et al, 2014)  and mapped to the human genome (hs37d5) using bowtie2  (Langmead et al, 2009)   genome were also removed. For further analysis and calculation of coverage tracks, all reads were shifted to be centered on the cutting position of the transposase and extended to a total of 29 bp to represent the region occupied by the transposase  (Adey et al, 2010). For bulk ATAC-seq, peak calling was done with MACS2 version 2. 1. 1  (Zhang et al, 2008)  with the parameters "-nomodel -shift -10extsize 28 -broad". Differential ATAC-seq signals were calculated with DiffBind (Ross-Innes et al, 2012), within peak regions called in at least four replicates (Dataset EV10.
1). Fold change thresholds were determined as the intersections of three Gaussian functions fitted to the distribution of log2 fold changes from DiffBind (Fig  EV1E ). The FDR threshold used for the differential ATAC-seq signal analysis was set to the x-value of the inflection point of the number of differentially accessible regions as function of FDR threshold. The analysis was also conducted with only fixed frozen NBC samples to ensure that the results were not affected by differences between fixed and viably frozen control samples (Dataset EV10. 2, used for TF motif analysis in Appendix Fig  S5A  and B ). For scATAC-seq, typically 40,000-100,000 unique integrations per cell were obtained. Cells with very few (less than 0. 2-times the average) or very frequent (more than three standard deviations above the average) integrations were removed from the analysis.
Variations in accessibility for different TF binding motifs were determined with chromVAR  (Schep et al, 2017). Correlation coefficients for simultaneous insertions across all cells at promoters and enhancers were computed with the R-package RWire, which has been deposited at https://github. com/FabianErdel/RWire. Target promoters of enhancers within TAD-sized regions were identified based on the respective correlation values using a threshold (0. 22) derived from the spurious correlations obtained after randomly permuting insertions (see Appendix Fig S6 for details). The resulting set of promoter-enhancer pairs was used for downstream analysis and integrated into the gene regulatory network.
A AE 1 kb region around the ATAC/Bidi signal center was intersected with the predicted intergenic and genic active enhancer states 1 and 9. The size of AE 1 kb was selected based on the average extension of key enhancer marks (H3K27ac, H3K4me1, p300, Bidi; Chen et al, 2018). Subsequently, overlapping regions were merged. Consensus active enhancer lists for CLL and NBCs were generated from individual patient or NBC control lists by selecting regions which were identified in at least three of the available CLL patient sample and one (ATAC) or two (Bidi) of the NBC controls, for which both the ChromHMM annotation and ATAC-seq/Bidi data were available. (iii) SEs were predicted with the ROSE tool  (Whyte et al, 2013 ) by using the H3K27ac read signal within the overlap region of ATAC and H3K27ac peaks.
The Gencode v17 genes were defined as a background gene model for filtering transcription start sites within ROSE. Differential analysis of H3K27ac at enhancers was performed using DiffBind on the enhancer regions occurring in at least three samples with the original H3K27ac data for the individual samples. Differential loci were called using DiffBind with an FDR threshold of 0. 01. Functional assignment of differential SEs was done with  GREAT (McLean et al, 2010)  with modified parameters of 10 kb upstream, 5 kb downstream, and 2 Mb distal which annotates cis-regulatory regions.
TF motif analysis To identify TF motifs in regions of interest with gained or lost occupancy, we used the HOMER package with appropriate background controls for each individual data set  (Heinz et al, 2010). All TF motif analysis was conducted within the ATAC peak regions. For the TF footprinting analysis, loci with a certain motif were retrieved with HOMER's "annotatePeaks. pl" and "known motifs" analysis. Coverage tracks AE 100 bp around the center of the motifs were generated with ngs. plot  (Shen et al, 2014)  using combined reads from all CLL patients or NBC controls collapsed to the center of the transposase binding event. In Appendix Table  S2 , the TF sequence-binding motifs and selected factors are listed that were associated with at least one additional aberrant chromatin feature in CLL in addition to chromatin accessibility changes as detected by ATAC-seq.
RNA-seq data analysis Sequence reads were mapped to the human genome (hs37d5,) using STAR version 2. 3. 0  (Dobin et al, 2013)  with the parameters: -genome-Load LoadAndRemove -alignIntronMax 500000 -alignMatesGapMax 500000 -outSAMunmapped Within -outFilterMultimapNmax 1 -outFilterMismatchNmax 3 -outFilterMismatchNoverLmax 0. 3 -sjdbOverhang 50 -chimSegmentMin 15 -chimScoreMin 1 -chimS-coreJunctionNonGTAG 0 -chimJunctionOverhangMin 15. RPKM and TPM values were computed for the long and the total RNA. The mapped reads were quantified using HTSeq-count from the HTSeq framework version 0. 6. 0 using the intersection_nonempty option in a stranded fashion according to the gencode17 annotation  (Anders et al, 2015). Normalization and differential expression analysis were performed using DESeq2 (Love et al, 2014; Dataset EV11).
The activity values of regulators were calculated based on our RNA-seq data set by using the VIPER algorithm  (Alvarez et al, 2016). For this, the raw RNA-seq counts were normalized by variance-stabilizing transformation using DESeq2  (Love et al, 2014). Activities for 2,804 of the 3,862 regulators could be computed. To identify regulators which showed a significantly different activity between the CLL samples and the NBC samples, a two-sided Student's t-test was performed and the P-values were adjusted for multiple testing  (Benjamini & Hochberg, 1995)  yielding 1,378 differentially active regulators with an adjusted P-value < 0. 05. Differential expression of regulators target genes of the regulators was determined with DESeq2. Gene set enrichment analysis was performed using the R-package "gPro-fileR"  (Reimand et al, 2016)  with the B-cell network as universe to identify significantly enriched KEGG pathways  (Kanehisa & Goto, 2000)  for the identified core TFs and their target genes.
The network was visualized using Cytoscape version 3. 5. 1  (Shannon et al, 2003). Regulators with a differential activity change between the CLL and the NBC samples (P < 0. 05) and all differential expressed non-regulators (P < 0. 01) with a log fold change of À1. 7 < LFC > 1. 7 were selected for further analysis. To identify aberrant chromatin modifier in CLL vs. NBCs, the list of deregulated regulators/genes was intersected with the EpiFactors database  (Medvedeva et al, 2015). This list was then reduced to those genes that were linked to the deregulated chromatin features identified in our study (meC loss, H3K27me3 gain, H3K27ac loss/gain, and H3K4me3 loss/gain). All enhancer-promoter pairs with a correlation of 0. 22 or higher in the scATAC-seq analysis  (Appendix Fig S6D and F)  were used from the CLL and the NBC samples.
To remove duplicates, overlapping regions were put together and the gene promoters were mapped to gene symbols. The enhancer-promoter pairs were grouped into CLL only (C), NBCs only (H), and found in both groups (CH). For all these cases, the enhancers were numbered consecutively (C1-C1208 in case of the CLL-specific enhancers) and the information about the genomic positions was added as a node attribute. For all the genes in the B-cell network, the enhancer-gene interactions were integrated. For the CLL-specific network, the interactions between the core TFs (Appendix Table  S2 ), their deregulated target genes, and the selected chromatin modifier as well as the enhancers were extracted from the B-cell network (Datasets EV13 and EV14). Exemplarily, the part of the network around TCF4 and EBF1 is shown in Fig 6D.
Acknowledgements We thank the members of our laboratories and the  CancerEpiSys and PRECiSe consortia  for help and discussion. The main funding of this project was from the  CancerEpiSys  ( 0316049 ) and  PRECiSe  ( 031L0076 ) grants of the  German Federal Ministry of Education and Research (BMBF)  in the e: Bio program. The  DKFZ-Heidelberg Center for Personalized Oncology (DKFZ-HIPO)  provided additional sequencing ( HIPO_H038 ). Further support was via  BMBF  project  CancerTelSys  ( 01ZX1302B ,  01ZX1602B ) to KR and RK,  Wellcome Trust  grant  200733/Z/16/Z  to VBT and  DFG SFB1074  (subprojects  B1  and  B2 ) to SS and DM, and the NCT 3. 0 project ENHANCE to PL and KR. We thank  Andre Holzer ,  Martina Seiffert , and  Yevhen Vainshtein  for help and discussions;  Gideon Zipprich ,  Jules Kerssemakers , and  Nick Kepper  for data management; and the  DKFZ Genomics and Proteomics Core Facility  for sequencing services.
Data availability The data and computer code produced in this study are available from the following sources: All original sequencing data have been deposited at European Genome-phenome Archive under the accession number EGAS00001002518 (http://www. ebi. ac. uk/ega/). Processed WGBS, ChIP-seq, (sc)ATAC-seq, and RNA-seq data are available under GSE113336 (https://www. ncbi. nlm. nih. gov/ge o/query/acc. cgi?acc=GSE113336) at Gene Expression Omnibus as full bigBed/BigWig files. Metadata from the comprehensive analyses are provided as Datasets EV1-EV14 with the manuscript. They are described in.