Scientist Stories: Jeff Vierstra, Footprinting the Human Genome

Hi, I'm Jeff Earstra, and today I'm going to share with you a little bit of our work to build nucleotide resolution maps of transcription factor occupancy across the human genome. Eukaryotic genome is hierarchically packaged into chromatin, the fundamental unit being the nucleosome, which is 150 base pairs of DNA wrapped around a histone octamer. Long arrays of nucleosomes form chromatin fibers and then higher-order structures on top of that. Chromatin is punctuated by cis-regulatory elements, such as promoters, enhancers, silencers, that control the expression of their cognate genes. Within cis-regulatory DNA, transcription factors cooperatively bind in the place of the canonical nucleosome. These TFs recruit complexes that directly control transcription. Thus, transcription factors form the basic building blocks of regulatory DNA, and I would argue that a full mechanistic understanding of Gene regulation requires detailed maps of transcription factor occupancy across the whole genome. So how do we do this? To do this, we employ a strategy called DNase I mapping. Now, over 30 years ago, it was discovered that active regulatory DNA is exquisitely sensitive to cleavage and digestion by the nonspecific endonuclease DNase I. The digestion of chromatin with DNase I generates small DNA fragments that can be purified and and end-sequenced using massively parallel sequencing. Mapping these small fragments back to the genome enables the genome-wide detection of active regulatory DNA. Now, with sufficient sequencing depth, one can begin to visualize the activity of DNase I at single nucleotide resolution. So if we sequence these CD8 T cells to 200 million tags, you can see that DNase cleavage is not uniform with one of these DHSs, but the signal is attenuated by the binding and occupancy of individual transcription factors. For example, NRF1 here, as I'm showing. Now, DNase I footprinting capitalized on the fact that 50% of the total cleavages occur within accessible elements representing 1% of the human genome, highlighting the incredible signal-to-noise ratio inherent to DNase I mapping. Now, intuitively, DNase I footprinting reflects the outcome of a competition between transcription factors and DNase I for access to DNA. Footprinting works because the affinity of sequence-specific transcription factors for their specific binding sites is much greater than the affinity of DNase I for those particular sites. So highly occupied sites will result in marked protection or a DNase I footprint, while lowly occupied sites reflect the intrinsic sequence preference of DNase I itself. So in order to truly fulfill the potential of DNase I footprinting as a method for the de novo genome-wide detection of transcription factor occupancy, we want to do this in a manner that doesn't require any prior knowledge of putative binding sites, only examining the cleavage profiles of DNA-seq cells. And the trick here is to use an algorithm to identify footprints in the data itself. To do this is very conceptually simple. We're just going to find contiguous regions of the genome with cleavage imbalances. So one way you can intuitively think about this is sort of running a window across the genome consisting of a footprint core and the flanking regions. And anytime there's fewer cleavages in the core than the flank, you could call those putative footprints. However, this approach has a number of challenges that may confound footprint detection, namely the variability in cleavage rates at adjacent bases, as well as the variation accessibility amongst different DHSs in general. A major contributor to the variation in DNase I cleavage is that DNA nucleases in general have intrinsic sequence specificity. In other words, DNase I, for example, doesn't cut uniformly on naked DNA. What we know now is that it senses the minor groove width. And in a collaboration with Harmen Bustemak and Remo Roos, we show that the cleavage preference can be effectively explained by a 6-mer sequence model. Now, we took advantage of this and we developed a computational approach that incorporates both chromatin architecture and the empirical DNase I sequence preference to determine the expected nucleotide cleavage rate across the genome. And then for each dataset, we derived a statistical model for testing whether the observed cleavage rate at individual nucleotides deviated significantly from the expectation. So what you can do here is then you can do that statistical test and choose an appropriate cutoff and call footprints. Now we performed de novo footprint discovery independently on all of the data— all the datasets in this particular study, which is about 243, and we detected on average about 650,000 of footprints per dataset. Now, it should be noted that the number of footprints that you can detect is very dependent on sequencing depth. So the deeper you sequence, the more footprints you're going to find. Now, this is great, but we wanted to see if we could do even better. Specifically, we wanted to leverage the fact that there now exist hundreds of regulatory DNA maps encompassing diverse human tissues and cell types. Now, if we use all these cell types, comparative footprinting across these has the potential to really illuminate both the structure and function of regulatory DNA. However, here, systematic approach for the joint analysis of these genomic footprinting data really has been lacking. Now, as I mentioned, given the scale and the diversity of these cell types, we really wanted to develop a framework that could integrate hundreds of available, even thousands of available datasets to increase the precision and resolution of footprint detection. And also, we wanted to build a scaffold to build a common reference index of transcription factor occupancy across the genome, essentially to build a composite view of what individual regulatory DNA looks like. [PAUSE] So to do this, we implemented an empirical-based framework that estimates the posterior probability that a given nucleotide is footprinted. And we did this by incorporating a prior on the presence of a footprint in a likelihood model of cleavage at both occupied and unoccupied sites. Now, this worked amazingly well. As you can see here, if we just plot individual nucleotides across all these data— so the heat map here, I'm showing you the posterior value of that empirical Bayes approach. If we plot the individual nucleotides scaled by their footprint prevalence across all the samples, doing so precisely resolves the core recognition sequences for all these diverse TFs in the bottom. Now, what we wanted to do is we wanted to build this reference set of TF-occupied DNA. So what we did here is we applied the same consensus approach described by Vautrin that he used to build the DHS index. And what we did here is we correlated the overlapping footprint regions across individual datasets into a nucleotide-resolution consensus footprint map. We applied this approach to all DHSs detected in one or more of the 243 datasets in this study. And this collectively delineated approximately 4.6 million consensus footprints, so individual distinct footprints present in one or more cell types. And these footprints were populated within 1.6 million of the 3.3 million DHSs indexed. So slightly about 50% of the DHSs we discovered footprints in. Footprint occupancy across all data sets showed marked enrichment for the recognition sequences of the master regulatory TFs and their corresponding lineages.. For example, you see GATA1 footprints in erythroblast cells, HNF-alpha footprints in fetal intestine, and PAX6 footprints in fetal eye. Now we find enrichment in a cell type for virtually all major classes of DNA binding domain families, suggesting that very few families, or if any at all, are refractory towards DNase I footprints. For degenerate motifs, the same sequence is recognized by many distinct TFs, We observed highly cell-specific occupancy patterns that could be further decomposed into coherent groups that corresponded to cell type and function, as you can see here with two different E-box families. Now, because TF engagement creates alterations in DNA shape and protects underlying phosphate bonds from nuclease attack, we wondered to what extent fluctuations in the DNase I cleavage rates reflected the topology of the transcription factor DNA interface in mass. Here we focused on CTCF. CTCF is a well-known poly zinc finger, as you know, which has two clusters of zinc fingers that bind DNA separated by a hinge. And that hinge region is thought to mediate DNA bending. So what we did is we transposed the overall per nucleotide cleavage propensity, or the average cleavage propensity, onto the co-crystal of CTCF with DNA. And what you can see here is that this accurately traced all the features of the protein-DNA interaction interface known, including the focal hypersensitivity limited to a couple of nucleotides within the hinge region as thought to mediate bending and likely modifying the minor groove with accentuating DNase I cleavage. Now furthermore, we can plot the corrected cleavage counts for all CTCF sites genome-wide in T regulatory cells. And this reveals that these topological features are immediately evident even at the level of individual footprints on the genome. So we're looking at the structure of individual binding sites. We examined this further for a number of TFs and found the average footprint width for diverse TFs tightly tracked the width of their respective recognition sequences. So taken together, we think that demonstrates that the extended profile of the perinucleotide DNA cleavage is really reminiscent or reflective of the 3D structure of the regulatory DNA element. We're definitely interested in looking at this in the future. Another extremely powerful feature of genomic footprinting is its ability to reveal the logic of individual cis-regulatory elements. For example, in non-nervous cell types, occupancy of the repressor NRSF in the SCAMP5 promoter silences transcription. However, nervous tissues like bipolar neurons, NRSF is not expressed, not bound, and transcription of SCAMP5 occurs. To detect differential occupancy and regulation in an unbiased manner, we devised a perinucleotide differential occupancy test similar to those employed for differential gene expression analysis. Here in the SCAMP5 promoter, our tests identified the precise nucleotides differentially occupied in nervous tissues and cell types. This revealed the fundamental logic of the SCAMP5 promoter. In nervous tissues, Two activators, ZFX and TFAP2, replace REST and drive expression, whereas strong REST occupancy in non-nervous tissue silences transcription. While the vast majority of disease and trait-associated variation is non-coding, identifying the specific genetic variants that are likely to affect regulatory functions remained a significant challenge. Deep sequence coverage of individual DHS enables the de novo genotyping of regulatory variation and the simultaneous characterization of their functional effect. And we do this by quantifying and comparing the cleavage at each allele at heterozygous sites. Now our data collectively encompass 243 individual data sets, and these data sets are derived from 147 individuals. And de novo genotyping across all these individuals revealed 4 million variants. Of these 4 million, 1.65 million were heterozygous and had the power to accurately quantify allelic imbalance. Across individuals, we conservatively identified 117,000 variants that altered DNA accessibility on an individual allele. Notably, within DHSs, single nucleotide polymorphisms that were allelically imbalanced were markedly enriched in the core consensus footprints. Highlighting that footprints are pinpointing the functional sites on the genome. Genomic footprinting data provides a unique nucleotide-resolution view and interpretation of the impact of individual DNA variants on chromatin structure. As an example, I'm showing you a variant within a DHS identified within the intron of the gene EGHD1. Here the variant is a C to a G. Allelic resolution of the DNase I cleavages in heterozygous individuals reveals that the C allele results in no footprint, where the G allele results in a strong footprint. We can confirm this in homozygous individuals by looking at individuals that contain either the C or the G. And you can see here that in these particular cell lines, in retina and HD29 cells, that the C allele also has no footprint while the G allele contains a strong footprint. We can confirm this in homozygous individuals by performing our differential footprint test across the two alleles. This reveals the precise nucleotides that are affected by this variant and shows that the G allele is creating a strong NF-γX binding site. So this means that the derived allele is creating a gain-of-function binding on the genome. While this might seem surprising, this actually occurs about half of the time. You can find many more examples of this across the genome with variants affecting NFI-X occupancy resulting either the loss of binding or the gain of binding. Here in this plot, I'm showing you the relative DNase I protection homozygotes on the y-axis versus the proportion of cleavages on the reference allele or allelic imbalance in heterozygotes on the x-axis. And you can see here there's a strong consistency between the variant and these two different configurations. Consistent with this, we can also look at the variant effect on the recognition sequence, or what is the predicted energetic effect on the sequence alone. And what you can see here is that as variants that result in the loss of binding as measured by footprinting the alternate allele creates a weaker motif than the reference allele. In contrast, at the sites that result in the gain of binding, or the variants that result in the gain of binding, the alternative allele creates a stronger motif than the reference allele. So really, taken all together here, this shows that genomic footprinting provides an ultra-high-resolution view of regulatory variation, its impact on transcription factor occupancy. Given that genetic variation affecting chromatin accessibility is enriched within footprints and that trait-associated genetic variation localizes within DHSs in general, we wondered whether disease and trait-associated genetic variation would preferentially localize in genomic footprints versus non-footprinted regulatory DNA. To do this, we performed enrichment of the GWAS catalog SNPs after LD expansion in DHSs but outside footprints and then and increasing stringencies of consensus footprints. What we find here is a strong increase in enrichment in the strongest footprints, while no enrichment within DHS outside footprints, suggesting that the vast majority of trait-associated variation is mediated by variation within footprints. To gain a more accurate view of the enrichment of trait-associated variants in footprints, we compared the SNP-based trait heritability of individual traits. We did this using summary statistic data from individual GWAS studies from the UK Biobank. We applied partitioned LD score regression to compute the relative heritability of these variants within all DHSs versus footprints. What we found here after doing this analysis was a striking enrichment of variants that account for trait heritability in footprints versus DHSs, and most prominently in footprints for their corresponding cell types. For example, here, if you look at red blood cell counts in erythroid footprints, we see an enrichment of approximately 45-fold. So taken together, we conclude that the genetic signals from disease and trait-associated variants are emanating primarily from TF footprints within DHSs, and that the variants within footprints seem to be the major contributors of trait-based heritability in regulatory DNA. To wrap up, digital genomic footprinting is a structural readout of regulatory DNA topology and an incredibly powerful approach to map transcription factor occupancy. We've made some computational advances to integrate hundreds to thousands of different cell types and to build nucleotide-resolution consensus maps of transcription factor occupancy within millions of cis-regulatory elements. Digital genomic footprinting also enables the unbiased and novel discovery of regulatory logic and structure of individual loci. We can measure the effects using this approach of genetic variation on TF occupancy at nucleotide resolution. And using this, we were able to show that disease and trait-associated variation is specifically and preferentially enriched within footprints versus non-footprinted regulatory DNA. Now, there's many more findings in the paper, and I hope you guys give it a read and take a look. I would also like to say that the digital genomic footprinting data is now available in the UCSC browser through a track hub. So you can go to the UCSC browser. You can load up a public track hub. It's called Digital Genomic Footprinting from 243 Cell and Tissue Types. And when you load that up, you'll get something that looks like this, where you can see here there's footprinting data from individual samples. So this includes the per nucleotide cleavage counts, both the observed and the expected tracks, footprint calls at different FDR thresholds in individual samples. A consensus footprint track, and motif matches from the overlapping footprints. So give it a look. Finally, I would like to make some acknowledgments to the people that supported and helped in this project. First, I'd like to acknowledge 3 talented computational biologists that aided in the data analysis. John Lazar was critical in designing some of the statistical frameworks used for footprint detection. Shane Neff and Eric Haugen helped in data processing. I'd like to also specifically acknowledge the ENCODE data production team and the institute for generating such fabulous DNA-SWAN data. My two colleagues at the institute, John Stam and Valter Muleman. Also like to reiterate the data is available at both Zenodo, my personal website, and there's extensive code and documentation at GitHub. Finally, I'd like to acknowledge our funding sources, the NHGRI ENCODE project and a charitable donation from GlaxoSmithKline. Thank you.