Change point detection for clustered expression data

Background: To detect changes in biological processes, samples are often studied at several time points. We examined expression data measured at different developmental stages, or more broadly, historical data. Hence, the main assumption of our proposed methodology was the independence between the examined samples over time. In addition, however, the examinations were clustered at each time point by measuring littermates from relatively few mother mice at each developmental stage. As each examination was lethal, we had an independent data structure over the entire history, but a dependent data structure at a particular time point. Over the course of these historical data, we wanted to identify abrupt changes in the parameter of interest - change points. Results: In this study, we demonstrated the application of generalized hypothesis testing using a linear mixed effects model as a possible method to detect change points. The coefficients from the linear mixed model were used in multiple contrast tests and the effect estimates were visualized with their respective simultaneous confidence intervals. The latter were used to determine the change point(s). In small simulation studies, we modelled different courses with abrupt changes and compared the influence of different contrast matrices. We found two contrasts, both capable of answering different research questions in change point detection: The Sequen contrast to detect individual change points and the McDermott contrast to find change points due to overall progression. We provide the R code for direct use with provided examples. The applicability of those tests for real experimental data was shown with in-vivo data from a preclinical study. Conclusion: Simultaneous confidence intervals estimated by multiple contrast tests using the model fit from a linear mixed model were capable to determine change points in clustered expression data. The confidence intervals directly delivered interpretable effect estimates representing the strength of the potential change point. Hence, scientists can define biologically relevant threshold of effect strength depending on their research question. We found two rarely used contrasts best fitted for detection of a possible change point: the Sequen and McDermott contrasts

Adaptation of the Oxygen Sensing System during Lung Development

During gestation, the most drastic change in oxygen supply occurs with the onset of ventilation after birth. As the too early exposure of premature infants to high arterial oxygen pressure leads to characteristic diseases, we studied the adaptation of the oxygen sensing system and its targets, the hypoxia-inducible factor- (HIF-) regulated genes (HRGs) in the developing lung. We draw a detailed picture of the oxygen sensing system by integrating information from qPCR, immunoblotting, in situ hybridization, and single-cell RNA sequencing data in ex vivo and in vivo models. HIF1α protein was completely destabilized with the onset of pulmonary ventilation, but did not coincide with expression changes in bona fide HRGs. We observed a modified composition of the HIF-PHD system from intrauterine to neonatal phases: Phd3 was significantly decreased, while Hif2a showed a strong increase and the Hif3a isoform Ipas exclusively peaked at P0. Colocalization studies point to the Hif1a-Phd1 axis as the main regulator of the HIF-PHD system in mouse lung development, complemented by the Hif3a-Phd3 axis during gestation. Hif3a isoform expression showed a stepwise adaptation during the periods of saccular and alveolar differentiation. With a strong hypoxic stimulus, lung ex vivo organ cultures displayed a functioning HIF system at every developmental stage. Approaches with systemic hypoxia or roxadustat treatment revealed only a limited in vivo response of HRGs. Understanding the interplay of the oxygen sensing system components during the transition from saccular to alveolar phases of lung development might help to counteract prematurity-associated diseases like bronchopulmonary dysplasia

Heterogeneous Effects of Direct Hypoxia Pathway Activation in Kidney Cancer

<div><p>General activation of hypoxia-inducible factor (HIF) pathways is classically associated with adverse prognosis in cancer and has been proposed to contribute to oncogenic drive. In clear cell renal carcinoma (CCRC) HIF pathways are upregulated by inactivation of the von-Hippel-Lindau tumor suppressor. However HIF-1<b>α</b> and HIF-2<b>α</b> have contrasting effects on experimental tumor progression. To better understand this paradox we examined pan-genomic patterns of HIF DNA binding and associated gene expression in response to manipulation of HIF-1<b>α</b> and HIF-2<b>α</b> and related the findings to CCRC prognosis. Our findings reveal distinct pan-genomic organization of canonical and non-canonical HIF isoform-specific DNA binding at thousands of sites. Overall associations were observed between HIF-1<b>α</b>-specific binding, and genes associated with favorable prognosis and between HIF-2<b>α</b>-specific binding and adverse prognosis. However within each isoform-specific set, individual gene associations were heterogeneous in sign and magnitude, suggesting that activation of each HIF-<b>α</b> isoform contributes a highly complex mix of pro- and anti-tumorigenic effects.</p></div

HIF-1α and HIF-2α binding genes confer opposing prognosis in kidney cancer.

<p>(A) The genes nearest to re-expressed HIF-1<b>α</b> (blue bars) and overexpressed HIF-2<b>α</b> (red bars) binding sites were defined and examined for enrichment amongst genes annotated in different cancers using the Human Disease Ontology database (<a href="http://www.disease-ontology.org/" target="_blank">http://www.disease-ontology.org</a>).–log10 Binomial p-values are plotted for each set of HIF-binding genes in each type of cancer. Grey bar denotes p = 0.05 (-log10, 1.3) level of significance. HIF-2<b>α</b> nearest binding genes are consistently more significantly enriched amongst cancer-associated genes than are HIF-1<b>α</b> binding genes. (B) Differential HIF-1<b>α</b> binding genes or (C) differential HIF-2<b>α</b> binding genes were filtered for significant associations with overall survival and used to generate a weighted gene predictor of prognosis for each set of genes. Patients were then divided into those with above or below median values for each gene predictor and subjected to Kaplan-Meier survival analysis. The Cox proportional hazard model indicated a significant protective effect for patients with above median gene predictor values based on the HIF-1<b>α</b> binding genes. Conversely, patients with above median values for the HIF-2<b>α</b> binding gene predictor had a significantly worse prognosis.</p

Preferential distribution of AP-1 binding motifs at HIF-2α versus HIF-1α binding loci.

<p>In addition to a hypoxia response element (HRE) motif, analysis of sites binding endogenous and overexpressed HIF-2<b>α</b> also identified an AP-1 motif. For each site, the maximum normalized log likelihood ratio for the AP-1 motif in red and the HRE motif in blue is plotted on the vertical axis as a bar chart. A smooth spline cubic fit line is overlaid to show the trend. The smoothing parameter is automatically determined using a ‘leave-one-out’ cross validation as implemented by the Smooth.spline function in R. Sites were categorized as binding (A) re-expressed HIF-1<b>α</b>, (B) overexpressed HIF-2<b>α</b> and ranked according to the HIF-1<b>β</b> signal at each site. Spline fit curves are overlaid (solid/dashed lines) to indicate overall trends across both forward and reverse strands. (A) Sites binding re-expressed HIF-1<b>α</b> show specific enrichment (positive score) for the HRE motif that decreases as the HIF-1<b>β</b> signal falls. In contrast, these same sites show depletion of the AP-1 motif. (B) Sites binding overexpressed HIF-2<b>α</b> show enrichment of the HRE motif that declines more steeply as the HIF-1<b>β</b> signal falls. In contrast to sites binding re-expressed HIF-1<b>α</b>, those binding overexpressed HIF-2<b>α</b> show enrichment of the AP-1 motif that increases (and exceeds that seen for the HRE) as the HIF-1<b>β</b> signal falls.</p

HIF-2α overexpression in 786-O cells.

<p>(A) HIF-2<b>α</b> binding sites in the HIF-2<b>α</b> overexpressing cells were identified by peak calling and ranked on the vertical axis according to signal intensity. Heat maps of these sites (±5kb on horizontal axis) showing ChIP-seq read density for the indicated HIF subunits were generated for both the control cells (i, ii) and the cells with HIF-2<b>α</b> overexpressed (iii, iv). In contrast to re-expressed HIF-1<b>α</b>, overexpressed HIF-2<b>α</b> binds to a large number of sites (compare i and iii), without HIF-1<b>β</b> (compare iii and iv) and has little effect on the distribution of HIF-1<b>β</b> (compare ii and iv). (B) Biplot showing Principal Component Analysis (PCA) of ChIP-seq signal intensity (RPKM values) for both individual binding sites (dots) and HIF-subunits (vectors) across all HIF-binding sites identified in control cells and in HIF-2<b>α</b> overexpressing cells. Sites binding endogenous HIF-2<b>α</b> in control cells are shown in blue while sites binding re-expressed HIF-1<b>α</b> are shown in red, sites binding both are colored purple and the remaining sites are colored grey. PCA for HIF subunits shows that HIF-2<b>α</b> and HIF-1<b>β</b> co-vary more closely in the control cells (compare HIF2<b>α</b>(VA) and HIF1<b>β</b>(VA)) than in the overexpressing cells (compare HIF2<b>α</b>(2<b>α</b>OE) and HIF1<b>β</b>(2<b>α</b>OE)). (C) Histogram of the distance to nearest transcription start site (TSS) for HIF-2<b>α</b> binding sites in cells overexpressing HIF-2<b>α</b>. (D) HIF-2<b>α</b> binding sites in the HIF-2<b>α</b> overexpressing cells were categorized according to the class (Ensemble) of the nearest gene. The relative frequency of each class is shown by pie chart. Gene set enrichment analysis (GSEA) for the set of genes nearest to (E) HIF-2<b>α</b> binding sites in the control cells and (F) newly identified HIF-2<b>α</b> binding sites in the overexpressing cells, when genes are ranked according to fold-change and significance in mRNA expression following overexpression of HIF-2<b>α</b> (horizontal axis).</p

HIF-1α re-expression in 786-O cells.

<p>(A) HIF-1<b>α</b> binding sites in the HIF-1<b>α</b> re-expressing cells were identified by peak calling and ranked on the vertical axis according to signal intensity. Heat maps of these sites (±5kb on horizontal axis) showing ChIP-seq read density for the indicated HIF subunits were generated for both the control cells (i, ii) and the HIF-1<b>α</b> re-expressing cells with HIF-1<b>α</b> re-introduced (iii-v). The pattern of HIF-2<b>α</b> binding is minimally affected by the re-expression of full-length HIF-1<b>α</b> (compare i and iii). Sites binding re-expressed HIF-1<b>α</b> are largely co-occupied by HIF-1<b>β</b> (compare iv and v). (B) Biplot showing Principal Component Analysis (PCA) of ChIP-seq signal intensity (RPKM values) for both individual binding sites (dots) and HIF-subunits (vectors) across all HIF-binding sites identified in control cells and in HIF-1<b>α</b> re-expressing cells. Sites binding endogenous HIF-2<b>α</b> in control cells are shown in blue while sites binding re-expressed HIF-1<b>α</b> are shown in red, sites binding both are colored purple and the remaining sites are shown in grey. PCA for each subunit shows high co-variance between HIF-2<b>α</b> binding in the control cells and in the HIF-1<b>α</b> re-expressing cells (compare HIF2<b>α</b>(VA) and (HIF2<b>α</b>(1<b>α</b>RE)). This indicates only minimal change in the HIF-2<b>α</b> binding as a consequence of the HIF-1<b>α</b> re-expression. Conversely, the HIF-1<b>β</b> vector changes dramatically with HIF-1<b>α</b> re-expression (compare HIF1<b>β</b>(VA) with HIF1<b>β</b>(1<b>α</b>RE)) and aligns closely with the vector for re-expressed HIF-1<b>α</b> (HIF1<b>α</b>(1<b>α</b>RE)). The individual binding sites in the control and HIF-1<b>α</b> re-expressing cells (blue and red dots) aligned closely with their respective PCA vectors. (C) Histogram of the distance to nearest transcription start site (TSS) for HIF-1<b>α</b> binding sites in cells re-expressing HIF-1<b>α</b>. (D) HIF-1<b>α</b> binding sites in the re-expressing cells were categorized according to the class (Ensemble) of the nearest gene. The relative frequency of each class is show by pie chart. (E) Gene set enrichment analysis (GSEA) for the set of genes nearest to HIF-1<b>α</b> binding sites when genes are ranked according to fold-change and significance in mRNA expression following re-expression of HIF-1<b>α</b> (horizontal axis).</p