Penalized estimation of the Gaussian graphical model from data with replicates.

Gaussian graphical models are usually estimated from unreplicated data. The data are, however, likely to comprise signal and noise. These two cannot be deconvoluted from unreplicated data. Pragmatically, the noise is then ignored in practice. We point out the consequences of this practice for the reconstruction of the conditional independence graph of the signal. Replicated data allow for the deconvolution of signal and noise and the reconstruction of former's conditional independence graph. Hereto we present a penalized Expectation-Maximization algorithm. The penalty parameter is chosen to maximize the F-fold cross-validated log-likelihood. Sampling schemes of the folds from replicated data are discussed. By simulation we investigate the effect of replicates on the reconstruction of the signal's conditional independence graph. Moreover, we compare the proposed method to several obvious competitors. In an application we use data from oncogenomic studies with replicates to reconstruct the gene-gene interaction networks, operationalized as conditional independence graphs. This yields a realistic portrait of the effect of ignoring other sources but sampling variation. In addition, it bears implications on the reproducibility of inferred gene-gene interaction networks reported in literature.

Mechanical stress regulates bone regulatory gene expression independent of estrogen and vitamin D deficiency in rats.

Mechanical stress determines bone mass and structure. It is not known whether mechanical loading affects expression of bone regulatory genes in a combined deficiency of estrogen and vitamin D. We studied the effect of mechanical loading on the messenger RNA (mRNA) expression of bone regulatory genes during vitamin D and/or estrogen deficiency. We performed a single bout in vivo axial loading with 14 N peak load, 2 Hz frequency and 360 cycles in right ulnae of nineteen weeks old female control Wistar rats with or without ovariectomy (OVX), vitamin D deficiency and the combination of OVX and vitamin D deficiency (N = 10/group). Total bone RNA was isolated 6 hours after loading, and mRNA expression was detected of Mepe, Fgf23, Dmp1, Phex, Sost, Col1a1, Cyp27b1, Vdr, and Esr1. Serum levels of 25(OH)D, 1,25(OH)2 D and estradiol were also measured at this time point. The effect of loading, vitamin D and estrogen deficiency and their interaction on bone gene expression was tested using a mixed effect model analysis. Mechanical loading significantly increased the mRNA expression of Mepe, and Sost, whereas it decreased the mRNA expression of Fgf23 and Esr1. Mechanical loading showed a significant interaction with vitamin D deficiency with regard to mRNA expression of Vdr and Esr1. Mechanical loading affected gene expression of Mepe, Fgf23, Sost, and Esr1 independently of vitamin D or estrogen, indicating that mechanical loading may affect bone turnover even during vitamin D deficiency and after menopause.

Identification of Deregulated Pathways, Key Regulators, and Novel miRNA-mRNA Interactions in HPV-Mediated Transformation.

Next to a persistent infection with high-risk human papillomavirus (HPV), molecular changes are required for the development of cervical cancer. To identify which molecular alterations drive carcinogenesis, we performed a comprehensive and longitudinal molecular characterization of HPV-transformed keratinocyte cell lines. Comparative genomic hybridization, mRNA, and miRNA expression analysis of four HPV-containing keratinocyte cell lines at eight different time points was performed. Data was analyzed using unsupervised hierarchical clustering, integrated longitudinal expression analysis, and pathway enrichment analysis. Biological relevance of identified key regulatory genes was evaluated in vitro and dual-luciferase assays were used to confirm predicted miRNA-mRNA interactions. We show that the acquisition of anchorage independence of HPV-containing keratinocyte cell lines is particularly associated with copy number alterations. Approximately one third of differentially expressed mRNAs and miRNAs was directly attributable to copy number alterations. Focal adhesion, TGF-beta signaling, and mTOR signaling pathways were enriched among these genes. PITX2 was identified as key regulator of TGF-beta signaling and inhibited cell growth in vitro, most likely by inducing cell cycle arrest and apoptosis. Predicted miRNA-mRNA interactions miR-221-3p_BRWD3, miR-221-3p_FOS, and miR-138-5p_PLXNB2 were confirmed in vitro. Integrated longitudinal analysis of our HPV-induced carcinogenesis model pinpointed relevant interconnected molecular changes and crucial signaling pathways in HPV-mediated transformation.

NOX2 Expression Is Increased in Keratinocytes After Burn Injury.

Reepithelialization is crucial for effective wound repair in burn wounds. Reactive oxygen species (ROS) have shown to be important in this. Recent studies suggest that NOX proteins produce ROS in keratinocytes. In the present study, we have studied NOX proteins in burn wounds, including the effect of C1-esterase inhibitor (C1inh) hereon, which is the endogenous inhibitor of complement activity whereof we have shown previously that it also increased the rate of reepithelialization in burn wounds. Skin tissue derived from healthy control Wistar rats (n = 6) were compared with burn-injured rats, with (n = 7) or without C1inh treatment (n = 7). After 14 days, rats were terminated. From the burn-injured rats, the entire wound and nonburned skin from the hind leg, that is, internal control was excised. From the control rats, dorsal skin was excised. In these skin samples, NOX2 and NOX4 were analyzed immunohistochemically. In nonburned rats, NOX2 was found in keratinocytes in both the basal layer and suprabasal layer of the epidermis; and the number of NOX2-positive keratinocytes was 367/mm2 (254-378). In burned rats, the number of NOX2-positive keratinocytes was significantly increased in the newly forming epidermis in the burned area to 1019/mm2 (649-1172), especially in the suprabasal layer, but significantly decreased in remote nonburned skin to 22/mm2 (6-89). C1inh treatment counteracted these changes in epidermal NOX2 expression in burned rats, both in the burned area as in remote nonburned skin. No NOX4 expression was found in the epidermis in none of the groups. NOX2 expression was increased in keratinocytes in newly forming epidermis after burn injury. C1inh, a drug that increases the rate of reepithelialization, counteracted this effect. These results suggest a role for NOX2 in the reepithelialization of burn wounds.

Ridge estimation of network models from time-course omics data.

Time-course omics experiments enable the reconstruction of the dynamics of the cellular regulatory network. Here, we describe the means for this reconstruction and the downstream exploitation of the inferred network. It is assumed that one of the various vector-autoregressive models (VAR) models presented here serves as a reasonably accurate description of the time-course omics data. The models are estimated through ridge penalized likelihood maximization, accompanied by functionality for the determination of optimal penalty paramaters. Prior knowledge on the network topology is accommodated by the estimation procedures. Various routes that translate the fitted models into more tangible implications for the medical researcher are described. The network is inferred from the-nonsparse-ridge estimates through empirical Bayes probabilistic thresholding. The influence of a (trait of a) molecular entity at the current time on those at future time points is assessed by mutual information, impulse response analysis, and path decomposition of the covariance. The presented methodology is applied to the omics data from the p53 signaling pathway during HPV-induced cellular transformation. All methodology is implemented in the ragt2ridges package, freely available from the Comprehensive R Archive Network.

Reconstruction of molecular network evolution from cross-sectional omics data.

Cross-sectional studies may shed light on the evolution of a disease like cancer through the comparison of patient traits among disease stages. This problem is especially challenging when a gene-gene interaction network needs to be reconstructed from omics data, and, in addition, the patients of each stage need not form a homogeneous group. Here, the problem is operationalized as the estimation of stage-wise mixtures of Gaussian graphical models (GGMs) from high-dimensional data. These mixtures are fitted by a (fused) ridge penalized EM algorithm. The fused ridge penalty shrinks GGMs of contiguous stages. The (fused) ridge penalty parameters are chosen through cross-validation. The proposed estimation procedures are shown to be consistent and their performance in other respects is studied in simulation. The down-stream exploitation of the fitted GGMs is outlined. In a data illustration the methodology is employed to identify gene-gene interaction network changes in the transition from normal to cancer prostate tissue.

Prognostic modeling of oral cancer by gene profiles and clinicopathological co-variables.

Accurate staging and outcome prediction is a major problem in clinical management of oral cancer patients, hampering high precision treatment and adjuvant therapy planning. Here, we have built and validated multivariable models that integrate gene signatures with clinical and pathological variables to improve staging and survival prediction of patients with oral squamous cell carcinoma (OSCC). Gene expression profiles from 249 human papillomavirus (HPV)-negative OSCCs were explored to identify a 22-gene lymph node metastasis signature (LNMsig) and a 40-gene overall survival signature (OSsig). To facilitate future clinical implementation and increase performance, these signatures were transferred to quantitative polymerase chain reaction (qPCR) assays and validated in an independent cohort of 125 HPV-negative tumors. When applied in the clinically relevant subgroup of early-stage (cT1-2N0) OSCC, the LNMsig could prevent overtreatment in two-third of the patients. Additionally, the integration of RT-qPCR gene signatures with clinical and pathological variables provided accurate prognostic models for oral cancer, strongly outperforming TNM. Finally, the OSsig gene signature identified a subpopulation of patients, currently considered at low-risk for disease-related survival, who showed an unexpected poor prognosis. These well-validated models will assist in personalizing primary treatment with respect to neck dissection and adjuvant therapies.

Gene Network Reconstruction using Global-Local Shrinkage Priors.

Reconstructing a gene network from high-throughput molecular data is an important but challenging task, as the number of parameters to estimate easily is much larger than the sample size. A conventional remedy is to regularize or penalize the model likelihood. In network models, this is often done locally in the neighbourhood of each node or gene. However, estimation of the many regularization parameters is often difficult and can result in large statistical uncertainties. In this paper we propose to combine local regularization with global shrinkage of the regularization parameters to borrow strength between genes and improve inference. We employ a simple Bayesian model with non-sparse, conjugate priors to facilitate the use of fast variational approximations to posteriors. We discuss empirical Bayes estimation of hyper-parameters of the priors, and propose a novel approach to rank-based posterior thresholding. Using extensive model- and data-based simulations, we demonstrate that the proposed inference strategy outperforms popular (sparse) methods, yields more stable edges, and is more reproducible. The proposed method, termed ShrinkNet, is then applied to Glioblastoma to investigate the interactions between genes associated with patient survival.

Cytotoxic Effects of the Therapeutic Radionuclide Rhenium-188 Combined with Taxanes in Human Prostate Carcinoma Cell Lines.

OBJECTIVE: Rhenium-188-HEDP is an effective radiopharmaceutical for the treatment of painful bone metastases from prostate cancer. The effectiveness of the ß-radiation emitted by 188Re might be enhanced by combination with chemotherapy, using the radiosensitization concept. Therefore, the authors investigated the combined treatment of the taxanes, docetaxel and cabazitaxel, with 188Re in prostate carcinoma cell lines. MATERIALS AND METHODS: The cytotoxic effects of single and combined treatment with taxanes and 188Re were investigated in three human prostate carcinoma cell lines (PC-3, DU 145, and LNCaP), using the colony-forming assay. The half maximal effective concentration (EC50) of all individual agents was determined. The combined treatment was studied at 0.25, 0.5, 1, 2, and 4 times the EC50 of each agent. The interaction was investigated with a regression model. RESULTS: The survival curves showed dose-dependent cell growth inhibition for both the taxanes and 188Re. The regression model showed a good capability of explaining the data. It proved additivity in all combination experiments and confirmed a general trend to a slight subadditive effect. CONCLUSIONS: This proof-of-mechanism study exploring radiosensitization by combining 188Re and taxanes showed no synergism, but significant additivity. This encourages the design of in vivo studies. Future research should explore the potential added value of concomitant treatment of bone metastases with chemotherapy and 188Re-HEDP.

Ridge estimation of the VAR(1) model and its time series chain graph from multivariate time-course omics data.

Omics experiments endowed with a time-course design may enable us to uncover the dynamic interplay among genes of cellular processes. Multivariate techniques (like VAR(1) models describing the temporal and contemporaneous relations among variates) that may facilitate this goal are hampered by the high-dimensionality of the resulting data. This is resolved by the presented ridge regularized maximum likelihood estimation procedure for the VAR(1) model. Information on the absence of temporal and contemporaneous relations may be incorporated in this procedure. Its computational efficient implemention is discussed. The estimation procedure is accompanied with an LOOCV scheme to determine the associated penalty parameters. Downstream exploitation of the estimated VAR(1) model is outlined: an empirical Bayes procedure to identify the interesting temporal and contemporaneous relationships, impulse response analysis, mutual information analysis, and covariance decomposition into the (graphical) relations among variates. In a simulation study the presented ridge estimation procedure outperformed a sparse competitor in terms of Frobenius loss of the estimates, while their selection properties are on par. The proposed machinery is illustrated in the reconstruction of the p53 signaling pathway during HPV-induced cellular transformation. The methodology is implemented in the ragt2ridges R-package available from CRAN.

Aberrant methylation-mediated silencing of microRNAs contributes to HPV-induced anchorage independence.

Cervical cancer and a subset of anogenital and head-and-neck carcinomas are caused by high-risk types of the human papillomavirus (hrHPV). During hrHPV-induced malignant transformation keratinocytes become able to grow anchorage independently, a tumorigenic trait at least partly associated with inactivation of tumor suppressor genes. We used hrHPV-containing keratinocytes to investigate the role of DNA methylation-mediated silencing of microRNAs (miRNAs) in the acquisition of anchorage independence.Anchorage dependent (n=11) and independent passages (n=19) of 4 hrHPV-immortalized keratinocyte cell lines were treated with 2'-deoxy-5-azacytidine (DAC). Genome-wide miRNA expression profiles before and after treatment were compared to identify miRNAs silenced by methylation. Bisulfite sequencing and methylation-specific PCR showed increased methylation of hsa-mir-129-2/-137/-935/-3663/-3665 and -4281 in anchorage independent HPV-transformed keratinocytes and cervical cancer cell lines. Mature miRNAs derived from hsa-mir-129-2/-137/-3663 and -3665 showed functional relevance as they decreased anchorage independence in cervical cancer cell lines. Cervical (pre)cancerous lesions demonstrated increased methylation of hsa-mir-129-2/-935/-3663/-3665 and -4281, underlining the clinical relevance of our findings.In conclusion, methylation-mediated silencing of tumor suppressive miRNAs contributes to acquisition of an anchorage independent phenotype. This study further substantiates the importance of miRNAs during early stages of carcinogenesis and underlines their potential as both disease markers and therapeutic targets.

Better prediction by use of co-data: adaptive group-regularized ridge regression.

For many high-dimensional studies, additional information on the variables, like (genomic) annotation or external p-values, is available. In the context of binary and continuous prediction, we develop a method for adaptive group-regularized (logistic) ridge regression, which makes structural use of such 'co-data'. Here, 'groups' refer to a partition of the variables according to the co-data. We derive empirical Bayes estimates of group-specific penalties, which possess several nice properties: (i) They are analytical. (ii) They adapt to the informativeness of the co-data for the data at hand. (iii) Only one global penalty parameter requires tuning by cross-validation. In addition, the method allows use of multiple types of co-data at little extra computational effort. We show that the group-specific penalties may lead to a larger distinction between 'near-zero' and relatively large regression parameters, which facilitates post hoc variable selection. The method, termed GRridge, is implemented in an easy-to-use R-package. It is demonstrated on two cancer genomics studies, which both concern the discrimination of precancerous cervical lesions from normal cervix tissues using methylation microarray data. For both examples, GRridge clearly improves the predictive performances of ordinary logistic ridge regression and the group lasso. In addition, we show that for the second study, the relatively good predictive performance is maintained when selecting only 42 variables.

Transcriptomic Heterogeneity in Cancer as a Consequence of Dysregulation of the Gene-Gene Interaction Network.

Many pathways are dysregulated in cancer. Dysregulation of the regulatory network results in less control of transcript levels in the cell. Hence, dysregulation is reflected in the heterogeneity of the transcriptome: the more dysregulated the pathway, the more the transcriptomic heterogeneity. We identify four scenarios for a transcriptomic heterogeneity increase (i.e., pathway dysregulation) in cancer: (1) activation of a molecular switch, (2) a structural change in a regulator, (3) a temporal change in a regulator, and (4) weakening of gene-gene interactions. These mechanisms are statistically motivated, explored in silico, and their plausibility to occur in vivo illustrated by means of oncogenomics data of breast cancer studies.

tigaR: integrative significance analysis of temporal differential gene expression induced by genomic abnormalities.

BACKGROUND: To determine which changes in the host cell genome are crucial for cervical carcinogenesis, a longitudinal in vitro model system of HPV-transformed keratinocytes was profiled in a genome-wide manner. Four cell lines affected with either HPV16 or HPV18 were assayed at 8 sequential time points for gene expression (mRNA) and gene copy number (DNA) using high-resolution microarrays. Available methods for temporal differential expression analysis are not designed for integrative genomic studies. RESULTS: Here, we present a method that allows for the identification of differential gene expression associated with DNA copy number changes over time. The temporal variation in gene expression is described by a generalized linear mixed model employing low-rank thin-plate splines. Model parameters are estimated with an empirical Bayes procedure, which exploits integrated nested Laplace approximation for fast computation. Iteratively, posteriors of hyperparameters and model parameters are estimated. The empirical Bayes procedure shrinks multiple dispersion-related parameters. Shrinkage leads to more stable estimates of the model parameters, better control of false positives and improvement of reproducibility. In addition, to make estimates of the DNA copy number more stable, model parameters are also estimated in a multivariate way using triplets of features, imposing a spatial prior for the copy number effect. CONCLUSION: With the proposed method for analysis of time-course multilevel molecular data, more profound insight may be gained through the identification of temporal differential expression induced by DNA copy number abnormalities. In particular, in the analysis of an integrative oncogenomics study with a time-course set-up our method finds genes previously reported to be involved in cervical carcinogenesis. Furthermore, the proposed method yields improvements in sensitivity, specificity and reproducibility compared to existing methods. Finally, the proposed method is able to handle count (RNAseq) data from time course experiments as is shown on a real data set.

A test for comparing two groups of samples when analyzing multiple omics profiles.

BACKGROUND: A number of statistical models has been proposed for studying the association between gene expression and copy number data in integrated analysis. The next step is to compare association patterns between different groups of samples. RESULTS: We propose a method, named dSIM, to find differences in association between copy number and gene expression, when comparing two groups of samples. Firstly, we use ridge regression to correct for the baseline associations between copy number and gene expression. Secondly, the global test is applied to the corrected data in order to find differences in association patterns between two groups of samples. We show that dSIM detects differences even in small genomic regions in a simulation study. We also apply dSIM to two publicly available breast cancer datasets and identify chromosome arms where copy number led gene expression regulation differs between positive and negative estrogen receptor samples. In spite of differing genomic coverage, some selected arms are identified in both datasets. CONCLUSION: We developed a flexible and robust method for studying association differences between two groups of samples while integrating genomic data from different platforms. dSIM can be used with most types of microarray/sequencing data, including methylation and microRNA expression. The method is implemented in R and will be made part of the BioConductor package SIM.

Penalized differential pathway analysis of integrative oncogenomics studies.

Through integration of genomic data from multiple sources, we may obtain a more accurate and complete picture of the molecular mechanisms underlying tumorigenesis. We discuss the integration of DNA copy number and mRNA gene expression data from an observational integrative genomics study involving cancer patients. The two molecular levels involved are linked through the central dogma of molecular biology. DNA copy number aberrations abound in the cancer cell. Here we investigate how these aberrations affect gene expression levels within a pathway using observational integrative genomics data of cancer patients. In particular, we aim to identify differential edges between regulatory networks of two groups involving these molecular levels. Motivated by the rate equations, the regulatory mechanism between DNA copy number aberrations and gene expression levels within a pathway is modeled by a simultaneous-equations model, for the one- and two-group case. The latter facilitates the identification of differential interactions between the two groups. Model parameters are estimated by penalized least squares using the lasso (L1) penalty to obtain a sparse pathway topology. Simulations show that the inclusion of DNA copy number data benefits the discovery of gene-gene interactions. In addition, the simulations reveal that cis-effects tend to be over-estimated in a univariate (single gene) analysis. In the application to real data from integrative oncogenomic studies we show that inclusion of prior information on the regulatory network architecture benefits the reproducibility of all edges. Furthermore, analyses of the TP53 and TGFb signaling pathways between ER+ and ER- samples from an integrative genomics breast cancer study identify reproducible differential regulatory patterns that corroborate with existing literature.

Identification of lethal microRNAs specific for head and neck cancer.

PURPOSE: The prognosis of head and neck squamous cell carcinomas (HNSCC) remains disappointing and the development of novel anti-cancer agents is urgently awaited. We identified by a functional genetic screen microRNAs that are selectively lethal for head and neck cancer cells but not for normal cells. We further investigated the genes targeted by these microRNAs. EXPERIMENTAL DESIGN: A retroviral expression library of human microRNAs was introduced in HNSCC cell lines and normal oropharyngeal keratinocytes to identify tumor-selective lethal microRNAs. Potential downstream gene targets of these microRNAs were identified by gene expression profiling and validated by functional assays. RESULTS: We identified six microRNAs that selectively inhibit proliferation of head and neck cancer cells. By gene expression profiling and 3'-untranslated region (UTR) luciferase reporter assays, we showed that the ataxia telangiectasia mutated (ATM) gene is a common target for at least two and likely three of these microRNAs. Specific inhibition of ATM resulted in a similar tumor-specific lethal effect, whereas the phenotype was reverted in rescue experiments. CONCLUSIONS: These six microRNAs might be developed as novel anti-cancer agents and highlight ATM as an interesting novel therapeutic target for head and neck cancer.

Modeling the DNA copy number aberration patterns in observational high-throughput cancer data.

The process of occurrence of genomic aberrations over time in the genetic material of cancer cells reflects the progression of the cancer. Modern technologies like aCGH (array Comparative Genomic Hybridization) and MPS (Massive Parallel Sequencing) provide high-resolution measurements of DNA copy number aberrations, that reveal the full scale of genomic aberrations. A continuous time Markov chain model is proposed to describe the accumulation of aberrations over time. Time however is a latent variable (with the number of aberrations as a proxy). Integrating out time, yields the distribution of the observed DNA copy number data. The model parameters are estimated from high-dimensional DNA copy number data by means of penalized maximum pseudo- and likelihood and method of moments procedures. Having fitted the model, posterior time estimates of the advancement of each sample's cancer are obtained and the most likely locations of a sample's aberrations are predicted. The three estimation methods are compared in a simulation study. The paper closes with an application of the proposed methodology on cancer data.

CD98 marks a subpopulation of head and neck squamous cell carcinoma cells with stem cell properties.

Patients with advanced head and neck squamous cell carcinomas (HNSCCs) are often treated with concomitant chemotherapy and radiotherapy, but only 50% is cured. A possible explanation for treatment failure is therapy resistance of the cancer stem cells (CSCs). The application of compounds specifically targeting these CSCs, in addition to routinely used therapeutics, would likely improve clinical outcome. We demonstrate that the previously described monoclonal antibody K984 recognizes the CD98 cell surface protein, which is specifically expressed by cells forming the squamous basal cell layer, the region where the squamous stem cells reside. Moreover, CD98 is highly resistant to the proteolytic enzymes required for CSC enrichment procedures. We show that CD98(high) cells, in contrast to CD98(low) cells, are able to generate tumors in immunodeficient mice, indicating that CD98(high) cells have stem cell characteristics. Furthermore, the CD98(high) subpopulation expresses high levels of cell cycle control and DNA repair genes, while the CD98(low) fraction shows expression patterns that represent the more differentiated cells forming the bulk of the tumor. CD98 is a promising CSC enrichment marker in HNSCC. Our data support the CSC concept in head and neck cancer and the potential relevance of these cells for treatment outcome.

Genomic aberrations relate early and advanced stage ovarian cancer.

BACKGROUND: Because of the distinct clinical presentation of early and advanced stage ovarian cancer, we aim to clarify whether these disease entities are solely separated by time of diagnosis or whether they arise from distinct molecular events. METHODS: Sixteen early and sixteen advanced stage ovarian carcinomas, matched for histological subtype and differentiation grade, were included. Genomic aberrations were compared for each early and advanced stage ovarian cancer by array comparative genomic hybridization. To study how the aberrations correlate to the clinical characteristics of the tumors we clustered tumors based on the genomic aberrations. RESULTS: The genomic aberration patterns in advanced stage cancer equalled those in early stage, but were more frequent in advanced stage (p = 0.012). Unsupervised clustering based on genomic aberrations yielded two clusters that significantly discriminated early from advanced stage (p = 0.001), and that did differ significantly in survival (p = 0.002). These clusters however did give a more accurate prognosis than histological subtype or differentiation grade. CONCLUSION: This study indicates that advanced stage ovarian cancer either progresses from early stage or from a common precursor lesion but that they do not arise from distinct carcinogenic molecular events. Furthermore, we show that array comparative genomic hybridization has the potential to identify clinically distinct patients.