{'fraction_non_alphanumeric': 0.039521644755702995, 'fraction_numerical': 0.028607117991563256, 'mean_word_length': 4.821313624985099, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 6, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Enhancer–silencer transitions in the human genome

Dual-function regulatory elements (REs), acting as enhancers in some cellular contexts and as silencers in others, have been reported to facilitate the precise gene regulatory response to developmental signals in Drosophila melanogaster. However, with few isolated examples detected, dual-function REs in mammals have yet to be systematically studied. We herein investigated this class of REs in the human genome and profiled their activity across multiple cell types. Focusing on enhancer–silencer transitions specific to the development of T cells, we built an accurate deep learning classifier of REs and identified about 12,000 silencers active in primary peripheral blood T cells that act as enhancers in embryonic stem cells. Compared with regular silencers, these dual-function REs are evolving under stronger purifying selection and are enriched for mutations associated with disease phenotypes and altered gene expression. In addition, they are enriched in the loci of transcriptional regulators, such as transcription factors (TFs) and chromatin remodeling genes. Dual-function REs consist of two intertwined but largely distinct sets of binding sites bound by either activating or repressing TFs, depending on the type of RE function in a given cell line. This indicates the recruitment of different TFs for different regulatory modes and a complex DNA sequence composition of these REs with dual activating and repressive encoding. With an estimated >6% of cell type–specific human silencers acting as dual-function REs, this overlooked class of REs requires a specific investigation on how their inherent functional plasticity might be a contributing factor to human diseases.

Transcriptional silencers, exerting negative regulatory impact and counterbalancing positive regulatory elements (REs, such as enhancers), have long been recognized to be essential for fine-tuning gene regulation and precisely responding to cellular signals and environmental stimuli in metazoans (Jacob and Monod 1961;Brand et al. 1985;Johnson et al. 2015;Rojano et al. 2019;Halfon 2020). However, because of the inherent difficulty in assaying and identifying silencers, the research focused on silencers has been greatly overshadowed by studies targeting enhancers. With the scarcity of known examples, knowledge about silencers has lagged behind that of enhancers, promoters, and even that of insulators (Ngan et al. 2020). Recently, the advance of massively parallel reporter assays (MPRAs) has greatly facilitated the large-scale detection of silencers in human cells and ignited interest in silencers (Della Rosa and Spivakov 2020;Doni Jayavelu et al. 2020;Pang and Snyder 2020).
Despite recent progress, large-scale silencer maps are still restricted to very few cell lines in humans, specifically HepG2 liver carcinoma and K562 chronic myelogenous leukemia cell lines (Doni Jayavelu et al. 2020;Pang and Snyder 2020). This only offers a limited characterization of the effects of silencers in various biological contexts. Furthermore, although histone H3 lysine 27 acetylation (H3K27ac) and histone H3 lysine 4 mono-methylation (H3K4me1) have been widely used for enhancer identification, the most promising mark for silencer identification, histone H3 lysine 27 trimethylation (H3K27me3), is not specific enough to accurately separate silencers from enhancers (Della Rosa and Spivakov 2020;Gisselbrecht et al. 2020). In human embryonic stem cells (hESCs), almost two-thirds of H3K27me3 chromatin immunoprecipitation sequencing (ChIP-seq) peaks carry signals of the histone marks associated with transcriptional activation. In primary T cells, a differentiated cell type, half of the H3K27me3 ChIP-seq peaks show the same property. Furthermore, there is a lack of consistency among silencer maps reported by different groups. In the HepG2 cell line, for example, two independent sets of silencers, which were identified using massively parallel platforms and contain about 4000 elements each, have only four silencers in common (Doni Jayavelu et al. 2020;Pang and Snyder 2020). This small overlap could be partially explained by different silencer pooling strategies and suggests an abundance of silencers active in a human cell type. Inspired by the success of sequence-based deep learning models in the prediction of chromatin states, TF binding occupancies, and enhancers (Alipanahi et al. 2015;Zhou and Troyanskaya 2015;Thibodeau et al. 2018), as well as our prior work on silencer detection (Huang et al. 2019), we trained a sequence-based deep learning model to identify silencers and showed that our model greatly outperforms other methods for silencer detection.
After generating accurate genome-wide maps of cell-specific silencers, we decided to focus on dual-function REs (DFREs), which act either as silencers or enhancers in different cell lines and thus represent the versatile component of gene regulatory machinery. Limited DFREs have been reported in Drosophila melanogaster (Stathopoulos and Levine 2005) and mice (Bessis et al. 1997;Kehayova et al. 2011). More recently, DFRE abundance in D. melanogaster has been reported (Erceg et al. 2017;Gisselbrecht et al. 2020). Almost 25% of nearly 1000 Pho Recessive Complex (PhoRC)-bound silencers in mesodermal cells function as enhancers at other developmental stages (Erceg et al. 2017), whereas 22 out of 29 (76%) identified mesodermal silencers were shown to act as DFREs in large enhancer screens (Gallo et al. 2011;Bonn et al. 2012) and were later augmented by additional six DFREs picked up by a targeted enhancer analysis (thus, bringing the total to 28/29 [Gisselbrecht et al. 2020]). Although the estimates of the fraction of silencers acting as DFREs differed between these two studies, likely because different silencer subclasses were investigated, both studies clearly showed that DFREs are widespread, at least in D. melanogaster embryonic cells. In this study, we focused on human DFREs, investigating distinct functional characteristics of these elements in comparison to regular silencers/enhancers and molecular mechanisms enabling their functional transition across cellular contexts.

Convolutional neural network model accurately predicts silencers
We built a multiclass convolutional neural network (CNN) model with three nodes in the output layer representing silencers, enhanc-ers, and regulatory neutral DNA sequences (see Methods) ( Fig. 1A). Our CNN model achieved a respectable prediction performance for silencer and enhancer identification across multiple human cell types. The area under the curve of the receiver operating characteristic (AUC-ROC) ranged from 0.84 to 0.94 for enhancer detection and from 0.74 to 0.90 for silencer detection across six cell types (H1 hESCs, hematopoietic progenitor cells [HPCs], HepG2s, K562s, monocytes, and T cells) ( Fig. 1B; Supplemental Fig. S1).
To perform an independent assessment of the accuracy of silencer predictions, we applied our CNN model to experimentally identified silencers. On two sets of silencers active in K562 cells (see Methods), our CNN model achieved the performance of AUC-ROC ≥0.78 and AUC-precision recall curve (AUC-PRC) A B D C E Figure 1. A deep learning model to predict silencers. (A) A schematic of the deep learning model used to predict cell type-specific silencers and enhancers. The number of kernels or neurons in a layer is listed in parentheses. The input of the model consists of 1-kb genomic sequences, and the output is a set of three probabilities of the input sequence being a silencer (ys), an enhancer (ye), or a nonfunctional sequence (yn). (B) ROCs and PRCs of the prediction models in three human cell types. The results in the additional three cell lines are presented in Supplemental Figure S1. (C) ROC and PRC classification accuracy in prediction of ReSE (Pang and Snyder 2020) and MPRA (Doni Jayavelu et al. 2020) experimentally characterized silencers. (D) Sharpr-MPRA MaxPos scores in K562 cells. EN and SL are the predicted enhancers and silencers, respectively. NonEN represents the DNase-seq peaks that overlap with H3K27ac ChIP-seq peaks but have not been predicted as enhancers in this study. NonSL represents the H3K27me3 ChIP-seq peaks not predicted as silencers, and the DHS column corresponds to all DNase-seq peaks. The count under a group label (x-axis) is the size of the corresponding sequence set. (E) The density of GWAS SNPs in the predicted silencers (SL), the sequences with H3K27me3 peaks but no H3K27ac peaks (H3K27me3/-H3K27ac), and the predicted enhancers (EN). The background consists of the genomic sequences having the GC content and repeat density matching to T cell SLs. ( * * ) P < 10 −10 . White asterisks are the significance enrichments of GWAS SNPs compared with H3K27me3/-H3K27ac, and black asterisks are the significant enrichments in enhancers compared with silencers.
≥0.47 (Fig. 1C). It outperforms the linear support vector machine (SVM) model we developed previously (AUC-ROC ≥0.66 and AUC-PRC = 0.32) (Huang et al. 2019). The improved accuracy of the CNN model in comparison to the SVM model can be attributed to the superiority of CNNs in capturing spatial patterns in sequences and retrieving nonlinear connections among these patterns (Alipanahi et al. 2015;Zhou and Troyanskaya 2015;Koo and Ploenzke 2020). It is important to note that our method accurately predicts experimentally identified ReSE (Pang and Snyder 2020) and MPRA (Doni Jayavelu et al. 2020) silencers (AUC ROC = 0.81 and 0.78 for ReSE and MPRA silencers, respectively) ( Fig.  1C). Despite a small overlap between these two data sets (<0.2% of matching sequences), the silencers from these data sets share common DNA sequence motifs, such as the enrichment in TFBSs of the repressive TFs REST and SETDB1 (Supplemental Fig. S2).
Silencers predicted by the CNN model also feature similar TFBS enrichment (Supplemental Fig. S2), indicative of this model capturing the motifs shared across different silencer sets and showing a reliable performance on two otherwise distinct silencer sets. There are also differences in the TFBS enrichment profiles of these three silencer sets, which could reflect different types of silencers within them.
To profile human silencers, we applied the trained CNN model to all DNA sequences carrying open chromatin signals (i.e., DNase-ChIP peaks) or repressive histone marks (i.e., H3K27me3 ChIP-seq peaks) and labeled the sequences having a silencer prediction score large enough to correspond a less-than-0.1 false-positive rate (FPR) on test samples as silencers (see Methods). We also applied this scheme to identify enhancers. To this end, we predicted approximately 130,000 silencers and 120,000 enhancers per cell line (Supplemental Table S1; Supplemental Fig. S3). Among the false-positive silencer predictions, that is, those that carry no H3K27me3 signals but are predicted as silencers, 7.7% correspond to enhancer candidates, that is, the DNase-seq peaks having H3K27ac but no H3K27me3 signals in the corresponding cell type. The remaining 92.3% of false-positive silencer predictions correspond to the background, that is, genomic regions with no H3K27ac, H3K27me3, or DNase signals in the corresponding cell type. In comparison to the test sample sets in which 5.6% are enhancer candidates (Supplemental Fig. S4), these results suggest a bias of false-positive silencer predictions to enhancer candidates. This likely reflects shared sequence motifs of enhancers and silencers active in a particular cell type that correspond to the binding sites of multifunctional TFs (Berest et al. 2019). For example, the NFKB2 complex functions as both activator and repressor in T cells (Senftleben et al. 2001;Grinberg-Bleyer et al. 2018), and its TFBSs are enriched in both T cell silencers and enhancers (see the section "Distinct sequence syntax of the DFREs").
In a systematic high-resolution activation and repression profiling with MPRA (Sharpr-MPRA), the activating/repressive capabilities of individual nucleotides within REs have been measured by MaxPos scores (Ernst et al. 2016). A low/high MaxPos score indicates a strong repressive/activating regulatory effect. In K562 cells, predicted silencers show the lowest average MaxPos score among all elements, which is 4.7 times lower than the average MaxPos score of H3K27me3 ChIP-seq peaks not predicted as silencers (dubbed nonSLs, Student's t-test P = 4 × 10 −6 ) (Fig. 1D). Similarly, predicted enhancers display a higher average MaxPos score than other elements, including H3K27ac ChIP-seq peaks not predicted as enhancers (dubbed nonEN, Student's t-test P = 10 −22 ) (Fig. 1D). These experimental results provide an independent support of our silencer and enhancer predictions.
In the five tested cell lines (monocytes were not included in this analysis owing to the lack of gene expression data for this cell type), the predicted silencers were enriched in the proximity of low-expression genes (Supplemental Fig. S5). Additionally, the density of single-nucleotide polymorphisms (SNPs) detected in genome-wide association studies (GWASs) is greater in predicted silencers and enhancers than in the sequences carrying H3K27me3 but no H3K27ac ChIP-seq peaks for all examined cell lines but hESCs (binomial test, P < 10 −10 ) (Fig. 1E). Furthermore, we observed that silencers show higher genomic mappability than H3K27me3 regions overall (Supplemental Fig. S6). To exclude a potential ascertainment bias arising from mappability differences between different sequence sets, we restricted the GWAS SNP analysis to the elements with at least 50% of certainly mappable nucleotides (i.e., mappability score = 1) (Supplemental Fig. S7). As there is no qualitative change in the results, this further confirms the functional importance of the predicted silencers, at least in comparison to H3K27me3 regions and randomly selected background sequences with matching length, GC content, and repeat density. For simplicity, we refer to these predicted silencers/enhancers as silencers/enhancers for the remainder of the manuscript.

DFREs are associated with TFs and chromatin remodeling genes
In our search for REs capable of alternating between enhancer and silencer function, we focused on H1 hESC enhancers that transition to silencers in primary T cells, a type of differentiated lymphoid cells. Based on overlap with H1 hESC enhancers, the silencers in primary T cells have been split into DFRE and non-DFRE silencers (see Methods). The latter were termed regular silencers (SLrs). To address the function of DFREs, we compared DFREs with SLrs and hESC enhancers that do not function as enhancers or silencers in T cells (named as ENrs below).
In primary T cells, 11,888 silencers (6% of all silencers) are DFREs ( Fig. 2A). Although DFREs and SLrs show similar H3K27me3 intensities (Supplemental Fig. S8), 8% of these DFRE sequences are conserved across placental species (Siepel et al. 2005), which is significantly higher than the 4.1% of SLrs and 3.6% of the background (binomial test P < 10 −10 ) but lower than the 10.6% of ENrs (P < 10 −10 ) (Fig. 2B). Also, the DFREs harbor 5.06 common SNPs per kilobase, which is significantly lower than the 5.44 common SNPs per kilobase in SLrs (P < 10 −10 ) and is similar to that of ENrs (Fig. 2B). In addition, 7.96% of DFRE SNPs have a derived allele frequency (DAF) of greater than 0.9, which is significantly lower than the 9.23% of the SLr SNPs (P < 10 −10 ) (Supplemental Fig. S9). These results suggest that DFREs have been and are still evolving under strong purifying selection, which is indicative of their functional importance.
Next, we used genes flanking DFREs to assess their potential function. Compared with all silencers, the DFREs show a strong preference for the loci of genes participating in transcriptional regulation (such as mRNA transcription, hypergeometric test P = 2 × 10 −8 ) (Supplemental Fig. S10) and/or genes associated with binding activity (such as T cell receptor binding, P = 3 × 10 −5 ) (Fig. 2C). Overall, 19.2% of DFREs are located in the loci of genes encoding TFs and/or chromatin modifiers (CMs), which is significantly higher than the 16.7% of SLrs (binomial test P = 10 −11 ; see Methods) ( Fig. 2D). This trend strengthens in the loci of the top 1000 high-expression TFs and/or CMs by 1.3 times (P = 6 × 10 −5 , high-expression vs. all TFs/CMs). Furthermore, compared with ENrs, DFREs show a similar enrichment in the loci of all TFs/ CMs (P = 0.03) but a contrasting strong preference to TFs/CMs highly expressed in T cells (binomial test P = 10 −7 ) (Fig. 2D), suggesting a distinct association of DFREs with the regulation of the TFs/CMs specific to T cells. The TFs are the core of gene regulation networks, and CM enzymes, which reshape the chromatin topologies globally or locally to facilitate or impede the binding of TFs, are crucial for regulating the responses to cellular and external signals, especially during cell differentiation and in the context of immunity (Dixon et al. 2015;Tartey and Takeuchi 2015;Chen et al. 2020). The strong association of DFREs with these genes suggests a primary contribution of DFREs to governing the activity and identity of T cells.
Next, we used RNA-seq gene expression (Roadmap Epigenomics Consortium et al. 2015) and Hi-C contact data in hESCs and T cells (Jung et al. 2019;Yang et al. 2020 in DFRE function might result in change of the DFRE target gene. Among intergenic DFREs that contact only one of their flanking genes (named A) and not the other flanking gene (named X) in hESCs, we identified DFREs that also feature an expression level of A at least fivefold greater than the corresponding X and the average gene expression level in hESCs. This set of DFRESs thus constituted the most pronounced DFRE enhancer effects on a specific target gene in hESCs. We then observed that 84.6% of these A genes in T cells are expressed at a less than twofold greater level than both the corresponding X genes and the average gene expression in T cells, and 100% of these A genes are expressed at a less than fivefold greater level than the corresponding X genes and the average gene expression in T cells. This indicates a pronounced drop in the level of A gene expression upon DFRE transition consistent with the enhancer-silencer DFRE change with no change in the target gene. Therefore, these results strongly suggest that the target gene of DFREs remains largely unchanged upon enhancer-tosilencer transition. In addition, using the Hi-C data ( The lack of change in chromatin contacts to proximal genes further highlights that the function transition of DFREs is pivotal for precisely regulating the transcription of these T cell differentiation genes during the development of T cells.

DFREs are enriched in the T cell-related genetic variants
To evaluate the influence of DFREs, we next explored GWAS SNPs and their linkage-disequilibrium (LD) counterparts (r 2 > 0.8; see Methods). The DFREs and SLrs harbor 0.93 and 0.91 GWAS SNP per kilobase, respectively, which represents a significant enrichment over the 0.74 GWAS SNP per kilobase in the background sequences (i.e., randomly sampled genomic sequences with the length, GC content, and repetitive element density matching T cell silencers; binormal test P < 10 −10 ) (Fig. 2E). Furthermore, 10.6% of GWAS SNPs in DFREs are associated with immunity-related traits and 9.8% of SLr SNPs are immunity associated. Both percentages are higher than the 8.5% of ENr SNPs and 9.4% of background SNPs (P < 10 −4 ) (Fig. 2E), highlighting the contribution of these silencers, especially DFREs (P = 0.0002, DFREs vs. SLrs), to the regulation of the immune system genes. Analyzing individual GWAS traits, we observed that T cell DFREs are specifically enriched in SNPs linked to immune system disorders, such as polyangiitis with granulomatosis, Sjogren syndrome, and type I diabetes mellitus (binormal test P < 0.05) ( Fig. 2F; Supplemental Table  S2). Eighteen out of the top 50 DFRE-associated GWAS traits (36%) are immunity related (Fig. 2F, azure diamonds), which is significant given that there are ∼8% of immunity-related traits among all GWAS traits (hypergeometric test P = 10 −7 ). In addition, T cell DFREs host 0.98 whole-blood expression quantitative trait loci (eQTLs) per kilobase, which is significantly higher than the 0.91 in the SLrs, 0.89 in ENrs, and 0.8 in the background sequences (binormal test P < 10 −5 ) (Fig. 2E). Furthermore, the density of T cell eQTLs in the DFREs is 0.012 per kilobase, which is 1.5 times that in SLrs (P = 2 × 10 −6 ) and 1.7 times that in ENrs (P = 2 × 10 −8 ) and the background sequences (P = 9 × 10 −9 ). Compared with whole-blood eQTLs, T cell eQTLs show elevated enrichment in the T cell silencers, especially the DFREs (chisquared test P < 10 −5 ) (Fig. 2E), suggesting strong tissue specificity and pronounced regulatory impact of these silencers on target gene expression in T cells. Also, the DFREs harbor 0.009 eQTL detected in induced pluripotent stem cells (iPSCs) per kilobase, which is similar to that in ENrs (binomial test P = 0.13) but is 1.7 times that in SLrs (P = 2 × 10 −7 ) and 1.4 times that in the background (P = 0.002). Combined, the enrichment in immunity-associated SNPs and stem-cell-associated SNPs supports the functional duality of the DRFEs, a distinct feature among all tested elements.
DFRE mutations are more likely to negatively affect silencer activity than SLr mutations To further characterize the effects of silencer mutations, we compared the output of our CNN models for wild-type (WT) alleles to those for mutant alleles, defining the CNN-based silencing alteration score (CNN-SAS; see Methods) (Supplemental Fig. S12). A large positive value of CNN-SAS suggests that the corresponding mutation causes a strong decrease in repressive activity. To evaluate CNN-SASs, we first turned to the experimental results in reporter assay quantitative trait locus (raQTL) studies where the regulatory effect of a SNP mutation was quantified by the expression change of the reporter gene (van Arensbergen et al. 2019). Positive raQTL scores represent mutations causing a decrease in reporter gene expression. Using the HepG2 CNN model, we calculated the CNN-SASs of HepG2 raQTL mutations and observed that the absolute values of the CNN-SASs of raQTLs were significantly higher than those of non-raQTLs (Wilcoxon rank-sum test P = 4 × 10 −164 ; see Methods) ( Fig. 3A; Supplemental Fig. S13). We also observed a gradual decrease in raQTL scores with an increase in CNN-SASs (Fig. 3B). Among silencer raQTLs with the 5% highest CNN-SASs, 94% had negative raQTL scores, which is a significant enrichment given that 54% of all tested raQTLs had negative raQTL scores (binomial test P = 3 × 10 −25 ) (Fig. 3B). The fraction of negative raQTLs dwindled with the decrease of CNN-SASs, ending at 12% among the raQTLs with the 5% lowest CNN-SASs (P = 10 −25 compared with the expected on all tested raQTLs). The strong negative correlation between CNN-SAS and raQTL scores validates the ability of our method to quantify the strength of silencer-disrupting mutations.
In T cells, for which raQTL data are not available, we used eQTL data to evaluate CNN-SASs. T cell CNN-SASs are correlated with T cell eQTLs (Spearman's correlation r = 0.064, P = 0.008) (Fig. 3C), potentially confirming the ability of CNN-SASs to predict the regulatory impact of mutations in T cells. The weak correlation between CNN-SASs and whole-blood eQTLs (r = 0.003, P = 0.35) (Supplemental Fig. S14) could be attributed to the indirect measurement of regulatory activity using eQTLs, as eQTLs reflect association and not causality of noncoding mutations with gene expression owing to (1) LD in loci containing regulatory variants and (2) a cumulative effect of multiple REs on target gene expression (as opposed to a direct readout of single mutation effects in raQTL experiments).
We next used the T cell CNN-SASs to assess the impact of silencer mutations. A CNN-SAS was considered significant when its absolute value was larger than 1% of those of all possible single-nucleotide mutations in the silencers (Methods) (Supplemental Fig. S15). Compared with the ENrs and background sequences, DFREs and SLrs host more significant CNN-SASs among GWAS SNPs and eQTLs (binomial test P ≤ 0.002), reflecting the functional importance of the DFREs and SLrs to the regulation of T cell phenotypes and gene expression. In addition, compared with SLrs, the DFREs are enriched for significant CNN-SASs (i.e., potential damaging mutations) across a panel of mutation sets (binomial test P < 10 −5 ) (Fig. 3D). For example, in DFREs, 3.5% of T cell eQTLs SNPs correspond to a significant CNN-SAS (that is 3.5-fold greater than expected by chance; P < 10 −10 ). These results validate the utility of significant CNN-SASs in detecting most likely silencer-damaging mutations.
An example DFRE mutation having a significant CNN-SAS is rs12631656. Its T-to-C (i.e., major to minor allele) mutation corresponds to a CNN-SAS of 0.0816 (P = 0.003) (Fig. 3E). This SNP is located within a tight LD block (r 2 > 0.96) of rs6445975 and rs4681851, the two SNPs that have been reported to be significantly associated with systemic lupus erythematosus and systemic sclerosis (Harley et al. 2008;Mayes et al. 2014). This LD block hosts five DFRE SNPs, four SLr SNPs, and one enhancer SNP in T cells. Among these SNPs, rs12631656 is the only one with a significant CNN-SAS (Fig. 3E). The TF motif mapping shows that the mutation at rs12631656 potentially weakens the binding affinity of ARID5B and SOX13, two known repressors in T cells (Lefebvre 2010;Wang et al. 2020).

Distinct sequence syntax of the DFREs
To address function encryption in DFRE sequences, we considered two primary hypotheses: (1) there is a single set of TF binding sites (TFBSs) within each DFRE bound by a stable set of TFs with alternating repressor and activator functions, and (2) there are two dis-tinct sets of TFBSs within DFREs-one encoding silencer function and another encoding enhancer function. To test these hypotheses, we mapped the TFBSs in TF ChIP-seq peaks (see "TFBS prediction in TF ChIP-seq peaks" in the Supplemental Notes). With these TFBSs, we first characterized TFBS signatures of the DFREs for their activating and repressive functions. In T cells, although inactive ENrs are depleted of ChIP-seq TFBSs of all TFs, silencers (DFREs and SLrs) are depleted of ChIP-seq TFBSs of activators (such as JUND, IRF3, and IRF4) but are enriched in ChIP-seq TFBSs of ubiquitous repressors, including REST and repressors acting predominantly in blood cells, such as EBF1 ( Fig. 4A; Györy et al. 2012), which is in line with their repressive function. Functioning as enhancers in hESCs, DFREs and ENrs show a similar TFBS profile. They are depleted of ChIP-seq TFBSs of REST (P ≤ 0.02) but are enriched in ChIP-seq TFBSs of hESC-specific TFs, such as POU5F1 and NANOG (Fig. 4B). These results suggest that DFREs recruit a distinct set of TFs for activating and repressive functions. Furthermore, 70% of silencer TFBSs are silencer specific, and 45% of enhancer TFBSs are enhancer specific (Fig. 4C), indicating a presence of distinct DFRE fragments that establish activating and repressive functions.
The enrichment of CTCF ChIP-seq TFBSs in DFREs prompted us to compare DFREs and insulators. In hESCs, where the topologically associated domains (TADs) have been reported , 0.84% of the DFREs are located at the TAD boundaries, which are known to be enriched for insulators (Ong and Corces 2014). This fraction is comparable to 0.94% of all DNase-seq peaks (P = 0.25) and is significantly lower than the 1.12% of CTCF ChIP- seq TFBSs (P = 10 −11 ) (Supplemental Fig. S16A). These trends were also observed in T cells, where 0.95% of DFREs, 0.93% of all DNaseseq peaks (P = 0.34 vs. DFREs), and 1.31% of CTCF ChIP-seq TFBSs (P = 10 −11 vs. DFREs) (Supplemental Fig. S16B) are located within the T cell TAD boundaries detected in Yang's study . Also, in T cells, DFREs are enriched in the loci of the 1000 lowest expressed genes (9.1% of DFREs vs. 6.3% of the whole genome, binomial test P = 5 × 10 −20 ) (Supplemental Fig. S16C) and are depleted in the loci of the 1000 highest expressed genes (P = 2 × 10 −5 ). A reverse trend was observed for DFREs in hESCs (P ≤ 0.001) (Supplemental Fig. S16D). Similarly, the DFREs frequently target the lowly expressed genes in T cells (P = 10 −6 ) (Supplemental Fig. S16E) but highly expressed genes in hESCs (P = 4 × 10 −15 ) (Supplemental Fig. S16F). These results are consistent with the silencer function of DFREs in T cells and the enhancer function of DFREs in hESCs. On the other hand, CTCF ChIP-seq TFBSs show no consistent correlation with either highly or lowly expressed genes in both cell types (Supplemental Fig. S16), supporting a functional distinction between DFREs and CTCF-defined insulators. We also observed that the majority of DFRE TFBSs, specifically 60% of enhancer TFBSs and 62% of silencer TFBSs, are located within the central component (±200 bp from the midpoint) of the DFREs, which is significantly >43% of randomly scattered TFBSs (binomial test P < 10 −10 ) (Fig. 4D). Furthermore, 37% of the silencer TFBSs are located within 50 bp to their nearest enhancer TFBSs, which is much higher than 24% of randomly scattered TFBSs in DFREs (binomial test P = 10 −160 ) (Fig. 4E). Combined, these results advocate for intertwined and spatially proximal distribution of distinct activating and repressive TFBSs within DFREs.
Next, to account for TFs not yet profiled in ChIP-seq experiments, we used CNN-SASs to predict TFBSs as short segments overrepresented in high CNN-SASs (see "TFBS prediction using a CNN model" in the Supplemental Notes). More than 43% of the CNNpredicted TFBSs overlap the reported TF ChIP-seq peaks in both T cells and hESCs, which is more than 1.4 times that of randomly scattered TFBSs (Supplemental Fig. S17A) and, therefore, justifies the use of this predictive approach for additional TFBS discovery. Predicted TFBSs in T cell (silencer) DFREs are enriched with the binding motifs of repressors, such as SNAI1/2, REST, and LMO2 (Supplemental Fig.  S17B). Predicted TFBSs in hESC (enhancer) DFREs host binding motifs of developmental activators, including POU5F1, NANOG, and SOX4 (Supplemental Fig. S17C). The CNN-predicted TFBSs confirm and further strengthen the trends observed using ChIP-seq TF data. Fifty-nine percent of these predicted TFBSs reside within the central region (±200 bp from the midpoint) of DFREs, which is significantly more than the expected 45% of randomly scattered TFBSs (binomial test P = 10 −323 ). Also, 43% of silencer TFBSs are located within 50 bp of enhancer TFBSs, which is significantly more than the expected 32% of randomly scattered TFBSs (P = 10 −323 ) (Supplemental Fig.  S17). Again, these results suggest a tight intertwining between TFBSs active for opposite functions within DFREs as well as a "coencryption" within the central regions of DFREs.

Developmental dynamics of DFRE formation
To understand the progression in DFRE formation during development, we applied our analysis pipelines for identifying DFREs in primary HPCs, which represent an intermediate step in

Enhancer-silencer transitions
Genome Research 443 www.genome.org differentiation from embryonic stem cells to T cells. Twenty-seven percent of HPC silencers are DFREs, acting as enhancers in hESCs (Fig. 5A). Compared with HPC SLrs, these DFREs are more conserved across placental species (10.2% of DFREs vs. 7.4% of SLrs, binomial test P < 10 −10 ) and feature a lower density of common SNPs (4.93 SNPs/kb in DFREs vs. 5.45 SNPs/kb in the SLrs, P < 10 −10 ) (Fig.  5B). HPC DFREs also harbor 0.98 GWAS SNPs and 1.3 whole-blood eQTLs per kilobase, which is significantly higher than the 0.87 GWAS SNPs and 0.94 whole-blood eQTLs in the SLrs, 0.91 GWAS SNPs and 0.96 whole-blood eQTLs in the ENrs (i.e., the hESC enhancers acting neither as silencers nor as enhancers in HPCs), and 0.74 GWAS SNPs and 0.8 whole-blood eQTLs in the background sequences (P < 10 −5 ) (Fig. 5C). HPC DFREs are more highly enriched in GWAS and eQTL SNPs than the SLrs, ENrs, and background in both tested cell lines, suggesting their functional importance across cellular contexts. Overall, HPC DFREs feature functional characteristics very similar to T cell DFREs, with the main exception of approximately four times as many HPC silencers being DFREs compared with T cell silencers (binomial test P < 10 −5 ). Furthermore, among 18,965 HPC DFREs, 11.6% (2201) are the DFREs in T cells, which is significantly lower than the 18.6% of HPC SLrs and 46.5% of HPC enhancers shared by T cells (binomial test P < 10 −5 ) (Fig. 5D). The shrinking fraction of DFREs during differentiation and the strong cell specificity of DFREs in a mature cell suggest the loss of functional plasticity in REs upon differentiation and a relatively small fraction of REs preserving multifunctional activity after the differentiation.

Discussion
Advances in high-throughput sequencing and MPRAs vastly increase the knowledge of gene regulation in human cells. However, the detection and functional characterization of active silencers remain challenging (Della Rosa and Spivakov 2020;Halfon 2020). In this study, we developed a deep learning silencer model that accurately detects experimentally identified silencers and quantifies the impact of silencer mutations. Identified silencers predominantly populate the loci of low-expression genes and are enriched in GWAS SNPs associated with the cell type matching disease etiology. In principle, silencers exert the repressive impact mainly through reducing the activity of basal promoters or smothering the activity of proximal enhancers (Gisselbrecht et al. 2020). H3K27me3 modification has been evidenced to take part in both repressive mechanisms (Ogiyama et al. 2018;Cai et al. 2021). As such, the CNN models, which were built with H3K27me3 ChIPseq peaks as training silencer samples, predict both types of silencers. These predicted silencers can be further categorized by using chromatin interaction maps with the knowledge that the silencers that reduce promoter activity interact with promoters more frequently than the other silencers. Also, our models, centered on the H3K27me3-associated silencers, might be less accurate in identifying silencers not marked by H3K27me3. Further studies are needed for investigate the differences and similarities between H3K27me-based and H3K27me3-independent silencers. By combining silencer detection with enhancer detection, we observed enhancer-silencer transitions and found that 6% of the T cell silencers and 28% of HPC silencers are DFREs, functioning as enhancers in hESCs and changing their regulatory function during differentiation. Compared with regular silencers, DFREs feature greater evolutionary sequence conservation and are enriched in GWAS SNPs and eQTLs. Moreover, DFREs are preferentially distributed in the proximity of genes governing transcriptional regulation. DFRE mutations are more than 1.5 times as likely to significantly damage silencer activity (as measured using CNN-SASs) as regular silencer mutations. These results support the essentiality and a rather overlooked important role of DFREs in regulating the development and maintaining the fate of blood cell types (and, likely, other differentiated cell types not included in this study). Earlier studies in D. melanogaster embryonic cells, coupled with the scope of our study targeting a single differentiation path, suggest that a large fraction of developmental enhancers might be wired to act as silencers upon and after differentiation (Erceg et al. 2017;Gisselbrecht et al. 2020).
To gain insights into the mechanisms of DFRE functional duality, we compared the distribution of active TFBSs corresponding to silencer and enhancer functions of DFRE sequences. We observed largely distinct but proximal DFRE regions occupied by binding sites of opposing functions. Our findings suggest a "double" genomic regulatory encryption within DFRE sequences with independent active regulatory binding sites, reflective of the use of different TFs for establishing activating and repressive functions after effectively rejecting an alternative hypothesis of the same pool of TFs bound to DFREs and switching their function based on cellular context.
Silencers are an integral part of the regulatory machinery in metazoans shaping and fine-tuning gene regulatory programs. The difference in the function of DFREs reported in D. melanogaster was primarily attributed to different cell types. In our study, we focus on the functional switch from enhancer to silencer activity, which is associated with the development and maturation of T cells. Therefore, it is logical to assume that the presented DFREs are more likely to be associated with developmental abnormalities and carcinogenesis. It is also possible that a subset of reported DFREs might underlie the environmental response, and their mutations are more likely to have a phenotypic effect owing to their role in development. Here we report thousands of DFREs during development of T cells and uncover unique features and profound phenotypic/pathogenic contributions of these elements. These findings together suggest the wide spread of DFREs in various biological processes across species. Follow-up studies extending to other biological contexts and additional chromatin data will further characterize human DFREs and help understanding how the same parts of the genome are reused for multiple regulatory functions at different developmental time points and in different cellular contexts. This is a crucial step toward revealing how gene regulation swiftly and precisely adapts to changing of cellular contexts in various biological processes.

Building CNN models to predict silencers
We designed a deep learning model composed sequentially of five convolutional layers and two fully connected network layers (Fig.  1A). Each CNN layer is followed by a max-pooling and a dropout layer. The input of this model is 1000-bp-long DNA sequences, whereas the output layer contains three nodes, each representing one of the sequence classes: silencer, enhancer, and background. Each output node predicts the probability of a given sequence belonging to the corresponding sequence classes. We used the Python library Keras version 2.4.0 (https://github.com/kerasteam/keras) to implement our model. To deal with this multiclass classification task, categorical cross-entropy was used as the cost function. By minimizing this cost function, the parameters of the models were adjusted by using the gradient-based algorithm Adam and then Adagrad as implemented in Keras (see "CNN model" in the Supplemental Notes; Supplemental Fig. S18).
We built a CNN model for each of six cell types, including H1 hESCs, the K562 cell line, the HepG2 cell line, primary HPCs, pri-mary monocytes, and primary T cells from peripheral blood. During training, we set aside Chromosome 6 for validation and Chromosomes 7 and 8 for testing. All other autosomes and Chromosome X were used for training.
To identify silencers, we calculated the cutoff in the silencer prediction scores (silencer scores) generated by the CNN model, which corresponds to the FPR of 0.1 in test samples, among which background samples were randomly selected so that the ratio of background samples to enhancers/silencers was 9:1 in each tested cell type. Sequences that carry a DNase-seq peak or a H3K27me3 ChIP-seq peak and have a silencer score greater than this cutoff were labeled as silencers. A similar approach with the FPR of 0.1 in CNN model enhancer prediction scores (enhancer scores) was used for enhancer identification. DNaseseq peaks that carry a H3K27ac signal and have an enhancer score greater than this cutoff established the set of predicted enhancers. After applying this scheme to all tested cell types, we identified approximately 130,000 silencers and 120,000 enhancers per cell type (Supplemental Table S1). The size of these enhancer sets is comparable to those identified by ChromHMM (Ernst and Kellis 2017) as well as those collected in the Benton and Gao databases (Benton et al. 2019;Gao and Qian 2020).

Data for training and validating CNN models
We trained a CNN enhancer/silencer model independently for each cell type. Cell type-specific DNase-seq and histone modification ChIP-seq data sets have been used for training (all ChIP-seq data used in this study were downloaded from the NIH Roadmap Epigenomics Project) (Roadmap Epigenomics Consortium et al. 2015). H3K27me3 ChIP-seq peaks with neither H3K27ac nor H3K4me1/ 3 ChIP-seq peaks overlapping their 400-bp central region constituted the set of training silencers. DNase-seq peaks with an overlapping H3K27ac peak but no overlapping H3K27me3 ChIP-seq peaks within their 400-bp central region constituted the set of training enhancers. Background sequences were randomly selected from DNase-seq peaks and H3K4me1, H3K4me3, H3K27ac, and H3K27me3 histone modification ChIP-seq peaks detected in any cell line other than the tested cell line. CNN models were trained on 1-kb regions centered on enhancer, silencer, and background sequences. For T cells, we collected 416,773 DNase-seq peaks and 392,442 H3K27me3 and 129,500 H3K27ac "broad" peaks from the NIH Roadmap Epigenomics Project (Roadmap Epigenomics Consortium et al. 2015). This translated into 239,925 and 158,887 nonredundant putative silencer and enhancer sequences, respectively. The CNN models identified 219,971 T cell silencers and 156,451 T cell enhancers from these sequences.
In the K562 cell line, experimentally validated silencers were used to evaluate the accuracy of the silencer prediction model. Those silencers were acquired from two resources: 3796 K562 silencers identified with self-transcribing active regulatory region sequencing platform (STARR-seq) (Doni Jayavelu et al. 2020) and 3909 K562 silencers reported using the repressive ability of silencer element assay (ReSE) (Pang and Snyder 2020). After redundancy exclusion, we assembled 7701 distinct silencers experimentally validated in K562 cells. To apply the CNN model in which the input is 1000-bp sequences on these silencers, we extended these silencers into 1000-bp genomic sequences centered at the midpoint of tested segments. To prevent a bias from prior knowledge on the estimate of the model accuracy, we excluded all 1406 experimentally validated silencers that overlap training or validation sequences (either for building the SVM or the CNN model) from the computational model evaluation.
In addition, the results of Sharpr-MPRA were used to evaluate silencer/enhancer predictions (Ernst et al. 2016). A low Sharpr-MPRA MaxPos score represents a negative RE. In this study, a genomic sequence was assigned with a MaxPos score only when this sequence hosted the whole corresponding Sharpr-MPRA sequence.

Enhancer-silencer transitions during development
To identify enhancer-silencer transitions during development, we compared silencer profiles in mature cell types with the enhancer map in H1 hESCs. A T cell silencer overlapping an H1 hESC enhancer by >600 bp was considered a DFRE. During embryonic development, the differentiation of T cells is a stepwise process from ESCs to mesodermal cells, HPCs, common lymphoid progenitors, and finally T cells (Kumar et al. 2018). In adulthood, T cells are differentiated from HPCs.

Evaluating the repressive effect of mutations
To evaluate the impact of silencer mutations, we focused on ys and ye in a built CNN model (Fig. 1A) as they predict the capability of genomic sequences being silencers and activators, respectively. Given the WT and mutant allele (MU) of a mutation, the difference of ys and ye between these alleles directly measures the alternation of activating and repressing power of the host sequence. That is, the CNN-SAS of a given mutation is defined as where ys WT and ys MU are the probabilities of silencer activity in the WT and the MU sequences, respectively. ye WT and ye MU are the probabilities of enhancer activities of the corresponding sequences. In addition, we measured the silencing alteration caused by a mutation with the odd ratios of (ys WT − ye WT ) to (ys MU − ye MU ) after deriving the probability function of (ys − ye) and noticed the similar performance between CNN − SAS scores and odd-ratio scores (see "CNN-based silencing odds ratio of mutations" in the Supplemental Notes; Supplemental Fig. S12).
The CNN − SAS distribution of all possible single-nucleotide silencer mutations forms a T distribution with 69% of mutations having an absolute CNN − SAS value smaller than 0.01 (Supplemental Fig. S15). Using this distribution as a background, we evaluated the significance of observed CNN − SAS scores. A low value of the probability represents a mutation having a strong impact on the activity of the host silencer. A CNN − SAS score was considered significant if this probability was <1%.

Data and tools used for analysis
We downloaded the GWAS SNPs curated in the National Human Genome Research Institute (NHGRI) catalog as of January 2019 (Buniello et al. 2019). There were 70,578 unique SNPs in the catalog at that time. To account for the fact that GWAS SNPs might not be directly responsible for phenotypic alterations but might be in LD with untested causal genetic variants (Cantor et al. 2010), we further expanded the SNP set by adding SNPs in a strong LD (r 2 > 0.8) block with GWAS SNPs in at least one population from the 1000 Genomes Project Consortium. To that end, we acquired 1,409,462 GWAS SNPs associated with 2722 traits. SNPs showing a strong phenotypic influence are likely to be associated with multiple traits according to independent studies (Zhou et al. 2018). To account for this, GWAS SNPs were linearly weighted by the number of the traits linked to them when we evaluated the GWAS SNP density. The GWAS traits were identified as immunity related when they contain keywords such as allergy, arthritis, asthma, Graves disease, infection, inflammation, lupus, multiple sclerosis, Sjogren syndrome, type I diabetes mellitus, etc.
Whole-blood eQTLs used in this study were downloaded from the Genotype-Tissue Expression Project, GTEx database version 8 (The GTEx Consortium 2015). To further evaluate the tissue specificity of the predicted DFREs, we collected eQTLs detected in T cells from the Blueprint epigenome project (Chen et al. 2016). EQTLs in iPSC were downloaded from the study of DeBoever et al. (2017).
We downloaded 14,183 HepG2 raQTLs from the study (van Arensbergen et al. 2019). These SNPs are associated with a significant change in reporter gene activity in HepG2 cells. We also downloaded additional 14,183 SNPs that correspond to an insignificant change in the expression of reporter genes (i.e., non-raQTL mutations) as analysis background.
We used the Genomic Regions Enrichment of Annotation Tool (GREAT) (McLean et al. 2010) to evaluate Gene Ontology (GO) biological processes and molecular functions associated with different RE groups. The element-gene association setting was selected as "two nearest genes." We explored conserved elements predicted using the 46-way placental phastCons scores (Siepel et al. 2005). A notable overlap of a genome sequence with conserved elements is indicative of an evolutionary constraint imposed on that sequence. The density of SNPs in genomic sequences was used to assess the evolutionary pressure acting on genomic segments in the human lineage. A low SNP density corresponds to strong selective pressure. The SNPs reported by The 1000 Genomes Project Consortium et al. (2015) were used for this analysis.
HESC TAD boundaries were downloaded from the Topologically Associating Domain Knowledge Base (TADKB) database in which the boundaries, in the resolution of 10 kb, were detected by using the Gaussian mixture model and proportion test ). The chromatin contacts in hESCs were downloaded from the database assembled in the study of Jung et al. (2019). The TADs and genomic contacts in T cells were downloaded from the report of Yang et al. (2020). We used the detections for unstimulated T cells. The TAD boundaries were defined as the 10-kb-long genomic regions centering at all the ends of TADs.

Software availability
Custom Python scripts (training CNN models and predicting enhancers/silencers with a built CNN model) are available on GitHub (https://github.com/ncbi/SilencerEnhancerPredict) and as Supplemental Code.