Detection of new pioneer transcription factors as cell-type specific nucleosome binders

Reviewer #1 (Public Review):

Peng et al develop a computational method to predict/rank transcription factors (TFs) according to their likelihood of being pioneer transcription factors--factors that are capable of binding nucleosomes--using ChIP-seq for 225 human transcription factors, MNase-seq and DNase-seq data from five cell lines. The authors developed relatively straightforward, easy to interpret computational methods that leverage the potential for MNase-seq to enable relatively precise identification of the nucleosome dyad. Using an established smoothing approach and local peak identification methods to estimate positions together with identification of ChIP-seq peaks and motifs within those peaks which they referred to as "ChIP-seq motifs", they were able to quantify "motif profiles" and their density in nucleosome regions (NRs) and nucleosome depleted regions (NDRs) relative to their estimated nucleosome dyad positions. Using these profiles, they arrived at an odd-ratio based motif enrichment score along with a Fisher's exact test to assess the odds and significance that a given transcription factor's ChIP-seq motifs are enriched in NRs compared to NDRs, hence, its potential to be a pioneer transcription factor. They showed that known pioneer transcription factors had among the highest enrichment scores, and they could identify a number of relatively novel pioneer TFs with high enrichment scores and relatively high expression in their corresponding cell line. They used multiple validation approaches including (1) calculating the ROC-AUC and Matthews correlation coefficient (MCC) and generating ROC and precision-recall curves associated with their enrichment score based on 32 known pioneer TFs among their 225 TFs which they used as positives and the remaining TFs (among the 225) as negatives; (2) use of the literature to note that known pioneer TFs that acted as key regulators of embryonic stem cell differentiation had a highest enrichment scores; (3) comparison of their enrichment scores to three classes of TFs defined by protein microarray and electromobility shift assays (1. strong binder to free and nucleosomal DNA, 2. weak binder to free and nucleosomal DNA, 3. strong binding to free but not nucleosomal DNA); and (4) correlation between their calculated TF motif nucleosome end/dyad binding ratio and relevant data from an NCAP-SELEX experiment. They also characterize the spatial distribution of TF motif binding relative to the dyad by (1) correlating TF motif density and nucleosome occupancy and (2) clustering TF motif binding profiles relative to their distance from the dyad and identifying 6 clusters.

The strengths of this paper are the use of MNase-seq data to define relatively precise dyad positions and ChIP-seq data together with motif analysis to arrive at relatively accurate TF binding profiles relative to dyad positions in NRs as well as in NDRs. This allowed them to use a relatively simple odds ratio based enrichment score which performs well in identifying known pioneer TFs. Moreover, their validation approaches either produced highly significant or reasonable, trending results.

The weaknesses of the paper are relatively minor, and the authors do a good job describing the limitations of the data and approach.

Reviewer #2 (Public Review):

In this study, the authors utilize a compendium of public genomic data to identify transcription factors (TF) that can identify their DNA binding motifs in the presence of nuclosome-wrapped chromatin and convert the chromatin to open chromatin. This class of TFs are termed Pioneer TFs (PTFs). A major strength of the study is the concept, whose premise is that motifs bound by PTFs (assessed by ChIP-seq for the respective TFs) should be present in both "closed" nucleosome wrapped DNA regions (measured by MNase-seq) as well as open regions (measured by DNAseI-seq) because the PTFs are able to open the chromatin. Use of multiple ENCODE cell lines, including the H1 stem cell line, enabled the authors to assess if binding at motifs changes from closed to open. Typical, non-PTF TFs are expected to only bind motifs in open chromatin regions (measured by DNaseI-seq) and not in regions closed in any cell type. This study contributes to the field a validation of PTFs that are already known to have pioneering activity and presents an interesting approach to quantify PTF activity.

For this reviewer, there were a few notable limitations. One was the uncertainty regarding whether expression of the respective TFs across cell types was taken into account. This would help inform if a TF would be able to open chromatin. Another limitation was the cell types used. While understandable that these cell types were used, because of their deep epigenetic phenotyping and public availability, they are mostly transformed and do not bear close similarity to lineages in a healthy organism. Next, the methods used to identify PTFs were not made available in an easy-to-use tool for other researchers who may seek to identify PTFs in their cell type(s) of interest. Lastly, some terms used were not define explicitly (e.g., meaning of dyads) and the language in the manuscript was often difficult to follow and contained improper English grammar.

Reviewer #3 (Public Review):

Peng et al. designed a computational framework for identifying pioneer factors using epigenomic data from five cell types. The identification of pioneer factors is important for our understanding of the epigenetic and transcriptional regulation of cells. A computational approach toward this goal can significantly reduce the burden of labor-intensive experimental validation.

The authors have addressed my previous comments.

The main issue identified in this re-review is based on the authors' additional experiments to investigate the reproducibility of the pioneer factors identified in the previously analysis that anchored on H1 ESCs.

The additional analysis that uses the other four cell types (HepG2, HeLa-S3, MCF-7, and K562) as anchors reveals the low reproducibility/concordance and high dependence on the selection of anchor cell type in the computational framework. In particular, now several stem cell related TFs (e.g. ESRRB, POU5F1) are ranked markedly higher when H1 ESC is not used as the anchor cell type as shown in Supplementary Figure 5.

Of note, the authors have now removed the shape labels that denote Yamanaka factors in Figure 2c (revised manuscript) that was presented in the main Figure 2a in the initial submission. The NFYs and ESRRB labels in Supplementary 4a are also removed and the boxplot comparing NFYs and ESRRB with other TF are also removed in this figure. Removing these results effectively hides the issues of the computational framework we identified in this revision. Please justify why this was done.

In summary, these new results reveal significant limitations of the proposed computational framework for identifying pioneer factors. The current identifications appear to be highly dependent on the choice of cell types.

Author Response

The following is the authors’ response to the original reviews.

Public Reviews:

Reviewer #1 (Public Review):

“Peng et al develop a computational method to predict/rank transcription factors (TFs) according to their likelihood of being pioneer transcription factors--factors that are capable of binding nucleosomes--using ChIP-seq for 225 human transcription factors, MNase-seq and DNase-seq data from five cell lines. The authors developed relatively straightforward, easy to interpret computational methods that leverage the potential for MNase-seq to enable relatively precise identification of the nucleosome dyad. Using an established smoothing approach and local peak identification methods to estimate positions together with identification of ChIP-seq peaks and motifs within those peaks which they referred to as "ChIP-seq motifs", they were able to quantify "motif profiles" and their density in nucleosome regions (NRs) and nucleosome free regions (NFRs) relative to their estimated nucleosome dyad positions. Using these profiles, they arrived at an odd-ratio based motif enrichment score along with a Fisher's exact test to assess the odds and significance that a given transcription factor's ChIP-seq motifs are enriched in NRs compared to NFRs, hence, its potential to be a pioneer transcription factor. They showed that known pioneer transcription factors had among the highest enrichment scores, and they could identify 32 relatively novel pioneer TFs with high enrichment scores and relatively high expression in their corresponding cell line. They used multiple validation approaches including (1) calculating the ROC-AUC associated with their enrichment score based on 16 known pioneer TFs among their 225 TFs which they used as positives and the remaining TFs (among the 225) as negatives; (2) use of the literature to note that known pioneer TFs that acted as key regulators of embryonic stem cell differentiation had a highest enrichment scores; (3) comparison of their enrichments scores to three classes of TFs defined by protein microarray and electromobility shift assays (1. strong binder to free and nucleosomal DNA, 2. weak binder to free and nucleosomal DNA, 3. strong binding to free but not nucleosomal DNA); and (4) correlation between their calculated TF motif nucleosome end/dyad binding ratio and relevant data from an NCAP-SELEX experiment. They also characterize the spatial distribution of TF motif binding relative to the dyad by (1) correlating TF motif density and nucleosome occupancy and (2) clustering TF motif binding profiles relative to their distance from the dyad and identifying 6 clusters.

The strengths of this paper are the use of MNase-seq data to define relatively precise dyad positions and ChIP-seq data together with motif analysis to arrive at relatively accurate TF binding profiles relative to dyad positions in NRs as well as in NFRs. This allowed them to use a relatively simple odds ratio based enrichment score which performs well in identifying known pioneer TFs. Moreover, their validation approaches either produced highly significant or reasonable, trending results.

The weaknesses of the paper are relatively minor. The most significant one is that they used ROC-AUC to assess the prediction accuracy of their enrichment score on a highly imbalanced dataset with 16 positives and 209 negatives. ROC-AUC is known to be a misleading prediction measure on highly imbalanced data. This is mitigated by the fact that they find an AUC = 0.94 for their best case. Thus, they're likely to find good results using a more appropriate performance measure for imbalanced data. Another minor point is that they did not associate their enrichment score (focus of Figure 2) with their correlation coefficients of TF motif density and nucleosome occupancy (focus of Figure 3). Finally, while the manuscript was clearly written, some parts of the Methods section could have been made more clear so that their approaches could be reproduced. The description of the NCAP-SELEX method could have also been more clear for a reader not familiar with this approach.”

Reviewer #2 (Public Review):

“In this study, the authors utilize a compendium of public genomic data to identify transcription factors (TF) that can identify their DNA binding motifs in the presence of nuclosome-wrapped chromatin and convert the chromatin to open chromatin. This class of TFs are termed Pioneer TFs (PTFs). A major strength of the study is the concept, whose premise is that motifs bound by PTFs (assessed by ChIP-seq for the respective TFs) should be present in both "closed" nucleosome wrapped DNA regions (measured by MNase-seq) as well as open regions (measured by DNAseI-seq) because the PTFs are able to open the chromatin. Use of multiple ENCODE cell lines, including the H1 stem cell line, enabled the authors to assess if binding at motifs changes from closed to open. Typical, non-PTF TFs are expected to only bind motifs in open chromatin regions (measured by DNaseI-seq) and not in regions closed in any cell type. This study contributes to the field a validation of PTFs that are already known to have pioneering activity and presents an interesting approach to quantify PTF activity.

For this reviewer, there were a few notable limitations. One was the uncertainty regarding whether expression of the respective TFs across cell types was taken into account. This would help inform if a TF would be able to open chromatin. Another limitation was the cell types used. While understandable that these cell types were used, because of their deep epigenetic phenotyping and public availability, they are mostly transformed and do not bear close similarity to lineages in a healthy organism. Next, the methods used to identify PTFs were not made available in an easy-to-use tool for other researchers who may seek to identify PTFs in their cell type(s) of interest. Lastly, some terms used were not defined explicitly (e.g., meaning of dyads) and the language in the manuscript was often difficult to follow and contained improper English grammar.”

Reviewer #3 (Public Review):

Peng et al. designed a computational framework for identifying pioneer factors using epigenomic data from five cell types. The identification of pioneer factors is important for our understanding of the epigenetic and transcriptional regulation of cells. A computational approach toward this goal can significantly reduce the burden of labor-intensive experimental validation. Nevertheless, there are several caveats in the current analysis which may require some modification of the computational methods and additional analysis to maximize the confidence of the pioneer factor prediction results.

A key consideration that arises during this review is that the current analysis anchors on H1 ESC and therefore may have biased the results toward the identification of pioneer factors that are relevant to the four other differentiated cell types. The low ranking of Yamanaka factors and known pioneer factors of NFYs and ESRRB may be due to the setup of the computational framework. Analysis should be repeated by using each of every cell type as an anchor for validating the reproducibility of the pioneer factors found so far and also to investigate whether TFs related to ESC identity (e.g. Yamanaka factors, NFYs and ESRRB) would show significant changes in their ranking. Given the potential cell type specificity of the pioneer factors, the extension to more cell types appears to be important for further demonstrating the utility of the computational framework.

Author Response: We thank all reviewers for their thoughtful and constructive comments and suggestions, which helped us to strengthen our paper. Following the suggestions, we have performed additional analysis to address the reviewer’s comments and the detailed responses are itemized below.

Reviewer #1 (Recommendations For The Authors):

1. The authors should generate precision-recall curves in addition to (or replacing) the ROC-AUC curves shown Figure 2c. They should also calculate the precision-recall AUC and use that as their measure of enrichment score predication accuracy. Precision-recall curves and AUC are more appropriate for imbalanced positive-negative data as is the case in this study.

Response: Following the reviewer’s suggestion, we have performed precision-recall analysis and calculated Matthews correlation coefficients (MCC) (Figure 2). We have further expanded our validation set to 32 known pioneer transcription factors (Supplementary Table 5) and compared the performance of enrichment score using different test sets (Supplementary Table 10). We have attained the highest ROC = 0.71, pr-ROC-AUC = 0.37 and MCC = 0.31 for Test set1 and ROC = 0.92, pr-ROC-AUC = 0.45 and MCC=0.49 for Test set2 (Supplementary Table 11).

1. The authors should generate scatter plots of their TF enrichment scores (focus of Figure 2) and motif-density nucleosome occupancy Pearson correlation coefficients (focus of Figure 3) and calculate the corresponding correlation coefficient and p-value.

Response: We observed a weak but statistically significant correlation between the enrichment scores and the correlation coefficient values (R=0.32 and p-value=1e-9)).

1. The authors should write their computational methods in the Methods section in such a way that a skilled bioinformatician could reproduce their results. This does not require a major rewrite. They are very close. One example of this is that a minimum distance between neighboring local maxima of the smoothed dyad counts was set to 150 bps. How was this algorithmically done? Suppress/ignore weaker local maxima that are within 150bp of other stronger local maxima?

Response: We have revised the Methods section to make it easier to follow and to reproduce the results. For identifying the local maxima, we have used the bwtool with the parameters ‘‘find local-extrema -maxima -min-sep=150’’ so that local maxima located within 150 bp of another neighboring maxima was ignored to avoid local clusters of extrema.

1. Describe the NCAP-SELEX method more clearly so that a reader not familiar with this approach doesn't have to look it up. This can be brief.

Response: Following the reviewer’s suggestion, we have added a detailed description of the NCAP-SELEX method.

Reviewer #2 (Recommendations For The Authors):

To improve the manuscript:

1. The grammar in the manuscript should be read for accuracy to improve readability and clarify the exact meaning.

Response: We have improved the grammar and have clarified the meaning of terms.

1. The exact meaning of dyads needs to be defined up front. In some places seems to mean pairs of reads and others seems to refer to nucleosome positioning.

Response: The meaning of “dyads” has been clarified. The dyad positions were determined by the midpoints of the mapped reads in MNase-seq data and refer to the center of the nucleosomal DNA.

1. Meaning of NCAP-SELEX needs to be defined before use of acronym.

Response: We have defined it in the manuscript.

Reviewer #3 (Recommendations For The Authors):

1. The authors found that Yamanaka factors and several other known pioneer factors (e.g. NFY-A, NFY-B, and ESRRB) are lowly ranked in their pioneer factor analysis. Since the analysis was performed by anchoring on H1 ESCs and comparing them to the other four cell lines, the results may only be relevant to differentiated cell types. It is therefore not unexpected that the Yamanaka factors which are important for iPSC reprogramming and the NFYs which have been experimentally shown to replace nucleosomes for maintaining ESC identity from differentiation (PMID: 25132174; PMID: 31296853) would not be enriched in the analysis. I suggest the authors repeat their analysis by anchoring on differentiated cell types and validate the reproducibility of the pioneer factors found so far and also investigate whether TFs related to ESC identity (e.g. Yamanaka factors, NFYs, and ESRRB) would show significant changes in their ranking as pioneer factors.

Response: Following reviewer’s suggestions, we have repeated the enrichment analysis by redefining differentially open regions as those closed in differentiated cell lines (HepG2, HeLa-S3, MCF-7 and K562) and open in H1 embryonic cell line (Supplementary Figure 6). The results indicate that most known PTFs still showed significantly higher enrichment scores compared with other TFs especially for FOXA, GATA and CEBPB families. Interestingly, ESSRB and Yamanaka pioneer factor POU5F1 (OCT4) have also shown significantly high enrichment scores in this analysis (Supplementary Figure 6). This could be explained by the roles of Yamanaka factors in cellular reprogramming – they reprogram somatic differentiated cells into induced pluripotent stem cells.

1. The authors mentioned the cell-type-specificity of TFs been pioneer factors and the example of CTCF was given. This point relates closely to above point 1 and, in particular, the correlation analysis of Yamanaka factors and NFYs supports their binding to nucleosomes. Together, these results highlight potential caveats of the current analysis in that the analysis is likely to be limited to the available cell types and may be affected by which cell type was used as the anchor cell type.

Response: Differentiated and embryonic cell lines were used to ask specific question about the functional roles of PTFs for cell differentiation and stem cell reprogramming. In the revised manuscript, we have clarified this point and separated our data set into three different sets of PTFs with different functions (Supplementary table 10). We agree with the reviewer, it would be nice to have more data from other cell lines but unfortunately the matching between different Chip-seq, DNAase-seq and Mnase-seq data sets imposes strict limitations.

1. The differential and conserved open chromatin regions are defined based on overlaps found between five cell types using their DNase-seq mapping profiles. The limitation of this definition is its lack of quantitativeness. For example, a chromatin region can have more than 80% overlaps between H1 and another cell type but the level of accessibility (e.g. number of reads mapped to this region) can be quite different between cell types. In such a case, I think it is still more appropriate to define such a region as a differential open chromatin region. The author should explore whether using a more quantitative definition would improve the identification and categorization of differential and conserved open chromatin regions.

Response: we thank the reviewer for these suggestions. In the revised version, we have clarified the definition and further explored different thresholds in defining the differentially and conserved open chromatin regions in enrichment analysis (Supplementary Figure 8). Our results were not significantly affected when different thresholds are applied.

1. While it is mentioned that H3K27ac and H3K4me1 ChIP-seq data from the five human cell lines were used in the study, the information on how enhancers are mapped/defined in these cell types is lacking.

Response: We have clarified the definition in the text. The enhancer regions were identified as the open chromatin regions overlapped with both H3K27ac and H3K4me1 ChIP-seq narrow peaks. We have elucidated the how enhancers are defined in the methods sections. In addition, we have performed additional enrichment analysis using NRs located on differentially active enhancer regions and NDRs located on conserved active enhancer regions (Supplementary Figure 7) between H1 embryonic cell line and any other differentiated cell lines and the performance of enrichment scores in PTF classification was slightly worse compared with those calculated from differentially and conserved open chromatin regions

1. The description of "genome-wide mapping of transcription factor binding sites" is unclear. For example, what does it mean by "In total, ChIP-seq data for 225 transcription factors could be matched with MNase-seq data" and why is this step needed? I would assume that a typical approach for mapping TF binding sites in the five cell types is to obtain the ChIP-seq data for each TF in each cell type and perform sequence alignment to the reference genome. The procedure described by the authors needs a clearer motivation and justification.

Responses: This sentence refers to matching between the ChIP-seq and MNase-seq data from the same cell type. We explain in detail how ChIP-seq data is processed. We have clarified this in the paper.

1. I also suggest the authors clearly justify the use of ROC analyses given that only a ground truth of positive (e.g. 16 known pioneer factors) is available and the "other transcription factors" considered as negative in the analysis in fact are expected to contain unknown pioneer factors and their identification should not be minimized (which lead to the maximization of ROC) by the analysis procedure.

Responses: (This is also pointed by review 1). The fact that unknown transcription factors are treated as negatives actually leads to the lower reported ROC scores (more hits considered to be false positives), not to their maximization. That is the reason we mentioned in the paper that the obtained ROC scores can be considered as lower bound estimates. In addition, we have expanded our validation sets to 32 known pioneer factors and compiled three sets of PTFs for validations. Following the reviewers’ suggestions, we have further performed precision-recall (PR) analysis and calculated the Matthews correlation coefficient (MCC) using three sets of PTFs for validation (Supplementary Table 11 and Supplementary Figure 2).

1. The analysis of pioneer transcription factor binding sites lacks insight. What can we learn these this analysis other than TFs from the same families are likely to be clustered in the same group?

Responses: We thank the review for pointing out it and have added a more detailed discussion of these results in the revised manuscript. Very few PTF-nucleosome structural complexes have currently been solved so far and the binding modes of majority of PTFs with nucleosomes still remain unknow. Our analysis has identified six distinctive clusters of TF binding profiles with nucleosomal DNA, which could provide insight into the binding modes of PTFs with nucleosome. These clusters point to the diversity of binding motifs where transcription factors belonging to the same cluster may also exhibit potential competitive binding.