<p>Figure S1. Experimental Figure Scheme. Figure S2. ChIP-Seq quality assessment for different samples and histone marks. Figure S3. Comparable ChIP-Seq and ChIP-qRT-PCR detection of histone enrichment for different histone marks near three reference genes. Figure S4. Comparable expression of reference genes. Figure S5. H3K9ac histone mark enrichment distribution near transcription start sites. Figure S6. H3K27ac histone mark enrichment distribution near transcription start sites. Figure S7. H3K9me3 histone mark enrichment distribution near transcription start sites. Figure S8. Correlation between tissue specific H3K9ac histone peak enrichment and expression of nearby genes. Figure S9. Correlation between tissue specific H3K9me3 histone peak enrichment and expression of nearby genes. Figure S10. Correlation between tissue specific H3K4me3 histone peak enrichment and expression of nearby genes for six study samples and two PDX-parental tissues. Figure S11. Correlation between tissue specific H3K27ac histone peak enrichment and expression of nearby genes for six study samples and two PDX-parental tissues. Figure S12. Correlation between tissue specific H3K9ac histone peak enrichment and expression of nearby genes for six study samples and two PDX-parental tissues. Figure S13. Correlation between tissue specific H3K9me3 histone peak enrichment and expression of nearby genes for six study samples and two PDX-parental tissues. Figure S14. Scheme for integration of ChIP-Seq Specific Histone Peaks with Expression Variation Analysis (EVA). Figure S15. Correlation of H3K27ac-enriched with genes differentially regulated in mesenchymal and classical HPV+ HNSCC subtypes.</p>

<p>Table S1. Clinical data for patient-derived PDX1, PDX2, UPPP1, and UPPP2 samples used for ChIP-based analysis. Table S2. ChIP-DNA qRT-PCR primers-probe sequence Table S3. Sample-specific enrichment of H3K4me3 histone mark at 5''UTR of individual genes. Table S4. Sample-specific enrichment of H3K27ac histone mark at 5''UTR of individual genes. Table S5. Sample-specific enrichment of H3K9ac histone mark at 5''UTR of individual genes. Table S6. Sample-specific enrichment of H3K9me3 histone mark at 5''UTR of individual genes. Table S7. Gene set enrichment analysis of genes linked to H3K27ac-enrichment specific for tumor samples. Table S8. Gene set enrichment analysis of genes linked to H3K27ac-enrichment specific for normal samples. Table S9. Gene set enrichment analysis of genes linked to H3K4me3-enrichment specific for tumor samples. Table S10. Gene set enrichment analysis of genes linked to H3K4me3-enrichment specific for normal samples. Table S11. Correlation of H3K27ac-enriched with genes differentially regulated in HPV-KRT and HPV-IMU.</p>

The script for ChIP-Seq and RNA-Seq data integration

Current literature suggests that epigenetically regulated super-enhancers (SEs) are drivers of aberrant gene expression in cancers. Many tumor types are still missing chromatin data to define cancer-specific SEs and their role in carcinogenesis. In this work, we develop a simple pipeline, which can utilize chromatin data from etiologically similar tumors to discover tissue-specific SEs and their target genes using gene expression and DNA methylation data. As an example, we applied our pipeline to human papillomavirus-related oropharyngeal squamous cell carcinoma (HPV + OPSCC). This tumor type is characterized by abundant gene expression changes, which cannot be explained by genetic alterations alone. Chromatin data are still limited for this disease, so we used 3627 SE elements from public domain data for closely related tissues, including normal and tumor lung, and cervical cancer cell lines. We integrated the available DNA methylation and gene expression data for HPV + OPSCC samples to filter the candidate SEs to identify functional SEs and their affected targets, which are essential for cancer development. Overall, we found 159 differentially methylated SEs, including 87 SEs that actively regulate expression of 150 nearby genes (211 SE-gene pairs) in HPV + OPSCC. Of these, 132 SE-gene pairs were validated in a related TCGA cohort. Pathway analysis revealed that the SE-regulated genes were associated with pathways known to regulate nasopharyngeal, breast, melanoma, and bladder carcinogenesis and are regulated by the epigenetic landscape in those cancers. Thus, we propose that gene expression in HPV + OPSCC may be controlled by epigenetic alterations in SE elements, which are common between related tissues. Our pipeline can utilize a diversity of data inputs and can be further adapted to SE analysis of diseased and non-diseased tissues from different organisms.

Head and neck squamous cell carcinomas (HNSCC) are the sixth leading cause of cancer worldwide with different incidences, mortalities, and prognosis for different subsites. Infection by Human Papilloma Viral (HPV) can cause the development of HPV+ HNSCC, the most rapidly growing population of cancer patients. The molecular biology of HPV+ HNSCC is related to abnormal transcriptional regulation. Super-enhancers (SEs) were recently identified as critical regulators of gene expression during cell differentiation and disease development. SEs are long clusters of enhancers that recruit master transcription factors (TFs) and coactivators to regulate the expression of target genes. We have recently developed a new whole-genome analytical pipeline to optimize ChIP-Seq procedures on primary surgical samples, and this protocol helped us to define the whole-genome distribution of H3K27ac, the hallmark of SEs, in HPV+ HNSCC samples. To detect SEs, we used a novel LILY algorithm that corrects the ChIP-Seq signal for copy number variation and CpG density. The analysis revealed that HPV+ HNSCC-specific SEs regulate genes critical for head and neck cancer development, such as EGFR, TP63, JAG1, IRF6, SMAD3, MAP4K2, SNAI1, TNFAIP3, ANO2, and TNFRSF1A. These genes are located in the vicinity of HNSCC-specific SEs, and the expression of those genes was altered by JQ1, the inhibitor of BRD4, the main component of SE machinery, supporting the role of SE in their regulation. The following Cistrome analysis allowed us to elucidate the TF composition of the SE machinery. Thus, we have found out that P63, P53, SOX2, FOSL1, SMAD3, and SNAI2 are TFs that are the most overrepresented in the HPV+ HNSCC-specific SEs. These comprehensive analyses have revealed novel insight into the HPV+ HNSCC biology and paved the basis for the implications for epigenetic therapeutics for this tumor types, and potentially other virus-related malignancies. Citation Format: Fernando T. Zamuner, Ilya Vorontsov, Emily Flam, Vera Mukhina, Ludmila Danilova, Tingting Ou, Elena Stavrovskaya, Theresa Guo, Eric Windsor, Dylan Z. Kelley, Michael Parfenov, Michael Considine, Elana J. Fertig, Alexander Favorov, Daria A. Gaykalova. Role of super-enhancers in HPV+ head and neck squamous cell carcinoma [abstract]. In: Proceedings of the American Association for Cancer Research Annual Meeting 2019; 2019 Mar 29-Apr 3; Atlanta, GA. Philadelphia (PA): AACR; Cancer Res 2019;79(13 Suppl):Abstract nr 5206.

Abstract Human papillomavirus-related head and neck squamous cell carcinoma (HPV+ HNSCC) is characterized by abundant gene expression changes, which cannot be explained by limited genetic alterations alone. We hypothesized that epigenetically regulated super-enhancers (SEs) are the drivers of aberrant gene expression in HPV+ HNSCC. We used integrated analysis of DNA methylation and gene expression data for cancer and noncancer oropharyngeal tissues to discover actionable SEs that regulate expression of target genes in HPV+ HNSCC. Since no SEs are currently defined for any HNSCC samples; we investigated 6196 SE elements that were obtained from the public domain for closely related tissues, including normal and tumor lung, and cervical cancer cell lines. We analyzed the methylation of these elements and gene expression of their nearby genes for 47 HPV+ HNSCC and their 25 noncancer controls. Overall, we found 122 differentially methylated SEs that had putative target (nearby) genes whose expression was correlated with enhancer methylation in HPV+ HNSCC. Of these, 107 were hypermethylated, and 15 were hypomethylated in tumors relative to normal samples. The pathway analysis revealed that the inferred SE-regulated genes were associated with pathways known to regulate carcinogenesis. Our data demonstrate that gene expression in HPV+ HNSCC may be regulated by epigenetic alterations in SE elements of related tissues. Citation Format: Emily L. Flam, Dylan Z. Kelley, Elena Stavrovskaya, Ludmila Danilova, Theresa Guo, Michael Considine, Jiang Qian, Joseph A. Califano, Alexander V. Favorov, Elana J. Fertig, Daria A. Gaykalova. Differentially methylated super-enhancers regulate gene expression in human papillomavirus-related head and neck squamous cell carcinoma [abstract]. In: Proceedings of the American Association for Cancer Research Annual Meeting 2018; 2018 Apr 14-18; Chicago, IL. Philadelphia (PA): AACR; Cancer Res 2018;78(13 Suppl):Abstract nr 364.

Abstract This project develops a novel experimental technique to perform ChIP-Seq (chromatin immunoprecipitation with massively parallel DNA sequencing) analysis of chromatin structure in primary tumor tissues from high risk HPV-related head and neck squamous cell carcinomas (HPV+ HNSCC). Recent data suggest that chromatin structure is the central regulator and predictor of cancer-specific expression and mutagenesis landscape of diseased cells. Genome-wide gene expression dysregulation in many tumors, including HPV+ HNSCC, are incompletely described by current knowledge. Methods for study of chromatin structure in primary tumor tissue are needed to better understand the role global epigenetic changes may play in these tumors. However, ChIP-Seq, which is the state-of-the-art method of elucidating chromatin structure, until now, has not been reliably performed on any HNSCC samples. Because chromatin structure is disrupted at room temperature, ChIP-Seq is especially complicated for primary patient tissues, which are primarily obtained as surgical waste after pathology review. Snap freezing of leftover waste surgical tissues and further tissue thawing for the analysis decreases chromatin structure integrity necessary for highly sensitive ChIP-Seq methodology, especially for tumor samples with chromatin structure deformed during carcinogenesis. To improve the chromatin structure integrity in tumor sample we added a xenografting step and minimized the exposure of cancer tissue to room temperature conditions after mouse surgery. We also minimized patient non-cancer tissue preservation at ambient temperature after patient surgery. We successfully performed ChIP-Seq for H3K4me3, H3K9me3, and H3K9ac on frozen uvulopalatopharyngoplasty (UPPP) primary tissues, frozen patient derived xenograft tissues, and freshly-cultured head and neck squamous cell carcinoma cell lines, revealing comparable success rates between tissue type and sample preservation techniques. ChIP-Seq techniques were performed and cross validated using tried and true qRT-PCR methods to demonstrate data reproducibility. The biological relevance of the ChIP-Seq data was confirmed through massive RNA-Seq analysis of 47 HPV+ HNSCC samples and 25 non-cancer controls. Analysis revealed that most H3K9ac and H3K9me3 enrichment is similar in primary tissues, regardless of disease status. Only small portion of them showed differential histone enrichment, which correlated with differential expression of corresponding genes. On the other hand, H3K4me3 showed strong tissue specificity and were found differentially enriched especially in tumor samples. The proposed experimental pipeline demonstrates high reproducibility between biological replicates, diversity of tissue models, and low dependence of ChIP-Seq analysis on tissue preservation techniques. Citation Format: Dylan Z. Kelley, Emily L. Flam, Hildegard A. Wulf, Theresa Guo, Evgeny Izumchenko, Dzov A. Singman, Ludmila V. Danilova, Elena D. Stavrovskaya, Michael Considine, Justin A. Bishop, William H. Westra, Zubair Khan, Wayne M. Koch, David Sidransky, Sarah Wheelan, Joseph A. Califano, Alexander V. Favorov, Elana J. Fertig, Daria A. Gaykalova. The in-parallel whole-genome ChIP-Seq analysis of primary tissues, patient derived xenografts, and cancer cell lines from HPV-relative HNSCC samples [abstract]. In: Proceedings of the American Association for Cancer Research Annual Meeting 2017; 2017 Apr 1-5; Washington, DC. Philadelphia (PA): AACR; Cancer Res 2017;77(13 Suppl):Abstract nr 2424. doi:10.1158/1538-7445.AM2017-2424

Abstract Chromatin alterations mediate mutations and gene expression changes in cancer. Chromatin immunoprecipitation followed by sequencing (ChIP-Seq) has been utilized to study genome-wide chromatin structure in human cancer cell lines, yet numerous technical challenges limit comparable analyses in primary tumors. Here we have developed a new whole-genome analytic pipeline to optimize ChIP-Seq protocols on patient-derived xenografts from human papillomavirus–related (HPV+) head and neck squamous cell carcinoma (HNSCC) samples. We further associated chromatin aberrations with gene expression changes from a larger cohort of the tumor and normal samples with RNA-Seq data. We detect differential histone enrichment associated with tumor-specific gene expression variation, sites of HPV integration in the human genome, and HPV-associated histone enrichment sites upstream of cancer driver genes, which play central roles in cancer-associated pathways. These comprehensive analyses enable unprecedented characterization of the complex network of molecular changes resulting from chromatin alterations that drive HPV-related tumorigenesis. Cancer Res; 77(23); 6538–50. ©2017 AACR.

Motivation: Genomics features with similar genome-wide distributions are generally hypothesized to be functionally related, for example, colocalization of histones and transcription start sites indicate chromatin regulation of transcription factor activity. Therefore, statistical algorithms to perform spatial, genome-wide correlation among genomic features are required. Results: Here, we propose a method, StereoGene, that rapidly estimates genome-wide correlation among pairs of genomic features. These features may represent high-throughput data mapped to reference genome or sets of genomic annotations in that reference genome. StereoGene enables correlation of continuous data directly, avoiding the data binarization and subsequent data loss. Correlations are computed among neighboring genomic positions using kernel correlation. Representing the correlation as a function of the genome position, StereoGene outputs the local correlation track as part of the analysis. StereoGene also accounts for confounders such as input DNA by partial correlation. We apply our method to numerous comparisons of ChIP-Seq datasets from the Human Epigenome Atlas and FANTOM CAGE to demonstrate its wide applicability. We observe the changes in the correlation between epigenomic features across developmental trajectories of several tissue types consistent with known biology and find a novel spatial correlation of CAGE clusters with donor splice sites and with poly(A) sites. These analyses provide examples for the broad applicability of StereoGene for regulatory genomics.

The modern high-throughput sequencing methods provide massive amounts of genome-focused, DNA-positioned data. This data is often represented as a function of the DNA coordinate (e.g. coverage). The genomeor chromosome-wide correlations between data from different sources may provide information about functional biological interrelation of the investigated features, e.g., trancription and histone modification. The task to compute the correlation was already successfully solved for interval annotations ([1]) as well as for coverage (functional) data ([2], [3], [4]). The key idea of the correlation studies is that two features that are similarly distributed along a chromosome may be functionally related. The point we are addressing here is a that peaks of dependent functional features can be located in a similar, although somewhat different, way. To account for these similarities, we propose here a fast method for calculation of kerneled correlation between two numeric annotations of the genome. The kernel represents the mutual position of related features; e.g., a Gaussian shape corresponds to ’somewhere around’, etc.

The modern high-throughput sequencing methods provide massive amounts of genomefocused, DNA-positioned data. This data is often represented as a function (e.g. coverage) onf the DNA coordinate. The genomeor chromosome-wide correlations between data from different sources may provide information about functional biological interrelation of the investigated features, e.g., about the trancription and histone modification. The task to compute the correlation was already successfully solved for interval annotations [1] as well as for coverage (functional) data ([2], [3], [4]). The key idea of the correlation studies is that two features that are similarly distributed along a chromosome may be functionally related. The point we are addressing here is a that peaks of dependent functional features can be located in a similar, although somewhat different, way. To account for these similarities, we propose here a fast method for calculation of kerneled correlation between two numeric annotations of the genome. The kernel represents the mutual position of related features; e.g., a Gaussian shape corresponds to 'somewhere around', etc. The approach is implemented as a computer program using C++ language. It allows counting of correlation not only for single features, but also for their combinations.

Identification of genes regulated by the same transcription factor (TF) is a major problem in analysis of regulation. The key step in detection of a group of co-regulated genes (regulon) is prediction of TF binding sites (TFBS). This is what positional weight matrix (PWM) is for. This matrix is applied to upstream region of a gene, and high-scoring sites are considered as putative TFBSs. Choice of threshold for the scoring function is a separate complicated problem. Usually, the threshold is chosen manually. Some methods for automated threshold detection exist, but they are based on selection of threshold for different functions. In this paper, we present an approach for regulon prediction based on a probabilistic method of threshold detection. The optimal probability computed by this method can be used to estimate the quality of the PWM itself. It can be useful when the matrix is a result of regulatory motif prediction program.

BACKGROUND:Genome-scale prediction of gene regulation and reconstruction of transcriptional regulatory networks in bacteria is one of the critical tasks of modern genomics. The Shewanella genus is comprised of metabolically versatile gamma-proteobacteria, whose lifestyles and natural environments are substantially different from Escherichia coli and other model bacterial species. The comparative genomics approaches and computational identification of regulatory sites are useful for the in silico reconstruction of transcriptional regulatory networks in bacteria.RESULTS:To explore conservation and variations in the Shewanella transcriptional networks we analyzed the repertoire of transcription factors and performed genomics-based reconstruction and comparative analysis of regulons in 16 Shewanella genomes. The inferred regulatory network includes 82 transcription factors and their DNA binding sites, 8 riboswitches and 6 translational attenuators. Forty five regulons were newly inferred from the genome context analysis, whereas others were propagated from previously characterized regulons in the Enterobacteria and Pseudomonas spp.. Multiple variations in regulatory strategies between the Shewanella spp. and E. coli include regulon contraction and expansion (as in the case of PdhR, HexR, FadR), numerous cases of recruiting non-orthologous regulators to control equivalent pathways (e.g. PsrA for fatty acid degradation) and, conversely, orthologous regulators to control distinct pathways (e.g. TyrR, ArgR, Crp).CONCLUSIONS:We tentatively defined the first reference collection of ~100 transcriptional regulons in 16 Shewanella genomes. The resulting regulatory network contains ~600 regulated genes per genome that are mostly involved in metabolism of carbohydrates, amino acids, fatty acids, vitamins, metals, and stress responses. Several reconstructed regulons including NagR for N-acetylglucosamine catabolism were experimentally validated in S. oneidensis MR-1. Analysis of correlations in gene expression patterns helps to interpret the reconstructed regulatory network. The inferred regulatory interactions will provide an additional regulatory constrains for an integrated model of metabolism and regulation in S. oneidensis MR-1.

RegPredict web server is designed to provide comparative genomics tools for reconstruction and analysis of microbial regulons using comparative genomics approach. The server allows the user to rapidly generate reference sets of regulons and regulatory motif profiles in a group of prokaryotic genomes. The new concept of a cluster of co-regulated orthologous operons allows the user to distribute the analysis of large regulons and to perform the comparative analysis of multiple clusters independently. Two major workflows currently implemented in RegPredict are: (i) regulon reconstruction for a known regulatory motif and (ii) ab initio inference of a novel regulon using several scenarios for the generation of starting gene sets. RegPredict provides a comprehensive collection of manually curated positional weight matrices of regulatory motifs. It is based on genomic sequences, ortholog and operon predictions from the MicrobesOnline. An interactive web interface of RegPredict integrates and presents diverse genomic and functional information about the candidate regulon members from several web resources. RegPredict is freely accessible at http://regpredict.lbl.gov.

Author(s): Novichkov, Pavel S.; Stavrovskaya, Elena D.; Gelfand, Mikhail s.; Mironov, Andrey A.; Dubchak, Inna; Rodionov, Dmitry A. | Abstract: One of the major challenges for the bioinformatics community in view of constantly growing number of complete genomes is providing effective tools to enable high-quality reconstruction of transcriptional regulatory networks (TRN). Definition of a particular TRN includes specification of which transcription factors (TF) bind to TF-binding sites (TFBS) in the promoter regions of which genes and what is the integrated effect of all these TFs on the expression of al these genes. Reconstruction of TRNs helps to better understand the metabolism and functions of bacteria. Among different approaches that are used for TRN reconstruction are an expression data-driven approach, and comparative genomic approaches that are either computing-driven, or subsystem (pathway) -driven. DNA microarrays, reporting gene expression, continue to be an important tool for high-throughput measurements on transcriptional levels, and machine-learning approaches were used to identify TRN (without a TFBS component) from a compendium of microarray expression profiles . However, in many cases the complexity of the interactions between regulons makes it difficult to distinguish between direct and indirect effects on transcription. Availability of a large number of complete genomes opens an opportunity to apply modern approaches of comparative genomics to expand the known regulons to yet uncharacterized organisms and to predict and describe new regulons with high precision.

There exist numerous algorithms for identification of regulatory signals in unaligned DNA fragments. Here we present two genetic algorithms for signal identification and describe their implementation and testing on simulated and real data. The first algorithm selects the start position of the signal in a given fragment. The second one builds a "universal" word that is recognized by the transcription factor. We compare these approaches and study the behavior of the genetic algorithm.

Author(s): Stavrovskaya, Elena D.; Rodionov, Dmitry A.; Mironov, Andrey A.; Dubchak, Inna; Novichkov, Pavel S. | Abstract: Reconstruction of transcriptional regulatory networks is one of the major challenges facing the bioinformatics community in view of constantly growing number of complete genomes. The comparative genomics approach has been successfully used for the analysis of the transcriptional regulation of many metabolic systems in various bacterial taxa. The key step in this approach is, given a position weight matrix, find an optimal threshold for the search of potential binding sites in genomes. Here we demonstrate that this problem is tightly bound to a problem of discovering the optimal content of regulon and suggest an approach to solve both problems simultaneously