Dynamic reprogramming of the kinome in response to targeted MEK inhibition in triple-negative breast cancer.

Kinase inhibitors have limited success in cancer treatment because tumors circumvent their action. Using a quantitative proteomics approach, we assessed kinome activity in response to MEK inhibition in triple-negative breast cancer (TNBC) cells and genetically engineered mice (GEMMs). MEK inhibition caused acute ERK activity loss, resulting in rapid c-Myc degradation that induced expression and activation of several receptor tyrosine kinases (RTKs). RNAi knockdown of ERK or c-Myc mimicked RTK induction by MEK inhibitors, and prevention of proteasomal c-Myc degradation blocked kinome reprogramming. MEK inhibitor-induced RTK stimulation overcame MEK2 inhibition, but not MEK1 inhibition, reactivating ERK and producing drug resistance. The C3Tag GEMM for TNBC similarly induced RTKs in response to MEK inhibition. The inhibitor-induced RTK profile suggested a kinase inhibitor combination therapy that produced GEMM tumor apoptosis and regression where single agents were ineffective. This approach defines mechanisms of drug resistance, allowing rational design of combination therapies for cancer.

A systematic assessment of cell type deconvolution algorithms for DNA methylation data.

We performed systematic assessment of computational deconvolution methods that play an important role in the estimation of cell type proportions from bulk methylation data. The proposed framework methylDeConv (available as an R package) integrates several deconvolution methods for methylation profiles (Illumina HumanMethylation450 and MethylationEPIC arrays) and offers different cell-type-specific CpG selection to construct the extended reference library which incorporates the main immune cell subsets, epithelial cells and cell-free DNAs. We compared the performance of different deconvolution algorithms via simulations and benchmark datasets and further investigated the associations of the estimated cell type proportions to cancer therapy in breast cancer and subtypes in melanoma methylation case studies. Our results indicated that the deconvolution based on the extended reference library is critical to obtain accurate estimates of cell proportions in non-blood tissues.

Efficient gradient boosting for prognostic biomarker discovery.

MOTIVATION: A gradient boosting decision tree (GBDT) is a powerful ensemble machine-learning method that has the potential to accelerate biomarker discovery from high-dimensional molecular data. Recent algorithmic advances, such as extreme gradient boosting (XGB) and light gradient boosting (LGB), have rendered the GBDT training more efficient, scalable and accurate. However, these modern techniques have not yet been widely adopted in discovering biomarkers for censored survival outcomes, which are key clinical outcomes or endpoints in cancer studies. RESULTS: In this paper, we present a new R package 'Xsurv' as an integrated solution that applies two modern GBDT training frameworks namely, XGB and LGB, for the modeling of right-censored survival outcomes. Based on our simulations, we benchmark the new approaches against traditional methods including the stepwise Cox regression model and the original gradient boosting function implemented in the package 'gbm'. We also demonstrate the application of Xsurv in analyzing a melanoma methylation dataset. Together, these results suggest that Xsurv is a useful and computationally viable tool for screening a large number of prognostic candidate biomarkers, which may facilitate future translational and clinical research. AVAILABILITY AND IMPLEMENTATION: 'Xsurv' is freely available as an R package at: https://github.com/topycyao/Xsurv. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

A network approach to compute hypervolume under receiver operating characteristic manifold for multi-class biomarkers.

Computation of hypervolume under ROC manifold (HUM) is necessary to evaluate biomarkers for their capability to discriminate among multiple disease types or diagnostic groups. However the original definition of HUM involves multiple integration and thus a medical investigation for multi-class receiver operating characteristic (ROC) analysis could suffer from huge computational cost when the formula is implemented naively. We introduce a novel graph-based approach to compute HUM efficiently in this article. The computational method avoids the time-consuming multiple summation when sample size or the number of categories is large. We conduct extensive simulation studies to demonstrate the improvement of our method over existing R packages. We apply our method to two real biomedical data sets to illustrate its application.

Evidence synthesis analysis with prioritized benefit outcomes in oncology clinical trials.

Overall survival, progression-free survival, objective response/complete response, and duration of (complete) response are frequently used as the primary and secondary efficacy endpoints for designs and analyses of oncology clinical trials. However, these endpoints are typically analyzed separately. In this article, we introduce an evidence synthesis approach to prioritize the benefit outcomes by applying the generalized pairwise comparisons (GPC) method, and use win statistics (win ratio, win odds and net benefit) to quantify treatment benefit. Under the framework of GPC, the main advantage of this evidence synthesis approach is the ability to combine relevant outcomes of various types into a single summary statistic without relying on any parametric assumptions. It is particularly relevant since health authorities and the pharmaceutical industry are increasingly incorporating structured quantitative methodologies in their benefit-risk assessment. We apply this evidence synthesis approach to an oncology phase 3 study in first-line renal cell carcinoma to assess the overall effect of an investigational treatment by ranking the most clinically relevant endpoints in cancer drug development. This application and a simulation study demonstrate that the proposed approach can synthesize the evidence of treatment effect from multiple prioritized benefit outcomes, and has substantial advantage over conventional methods that analyze each individual endpoint separately. We also introduce a newly developed R package WINS for statistical inference based on win statistics.

Polygenic association of glomerular filtration rate decline in world trade center responders.

BACKGROUND: The factors associated with estimated glomerular filtrate rate (eGFR) decline in low risk adults remain relatively unknown. We hypothesized that a polygenic risk score (PRS) will be associated with eGFR decline. METHODS: We analyzed genetic data from 1,601 adult participants with European ancestry in the World Trade Center Health Program (baseline age 49.68 ± 8.79 years, 93% male, 23% hypertensive, 7% diabetic and 1% with cardiovascular disease) with ≥ three serial measures of serum creatinine. PRSs were calculated from an aggregation of single nucleotide polymorphisms (SNPs) from a recent, large-scale genome-wide association study (GWAS) of rapid eGFR decline. Generalized linear models were used to evaluate the association of PRS with renal outcomes: baseline eGFR and CKD stage, rate of change in eGFR, stable versus declining eGFR over a 3-5-year observation period. eGFR decline was defined in separate analyses as "clinical" (> -1.0 ml/min/1.73 m2/year) or "empirical" (lower most quartile of eGFR slopes). RESULTS: The mean baseline eGFR was ~ 86 ml/min/1.73 m2. Subjects with decline in eGFR were more likely to be diabetic. PRS was significantly associated with lower baseline eGFR (B = -0.96, p = 0.002), higher CKD stage (OR = 1.17, p = 0.010), decline in eGFR (OR = 1.14, p = 0.036) relative to stable eGFR, and the lower quartile of eGFR slopes (OR = 1.21, p = 0.008), after adjusting for established risk factors for CKD. CONCLUSION: Common genetic variants are associated with eGFR decline in middle-aged adults with relatively low comorbidity burdens.

The Association of Posttraumatic Stress Disorder With Longitudinal Change in Glomerular Filtration Rate in World Trade Center Responders.

OBJECTIVE: High levels of psychological distress increase the risk of a wide range of medical diseases. In this study, we investigated the association between posttraumatic stress disorder (PTSD) and kidney disease. METHODS: World Trade Center (WTC) responders were included if they had two or more measures of estimated glomerular filtration rate (eGFR). The PTSD Checklist (PCL) was used to define no PTSD (PCL < 40), "mild" PTSD (40 ≤ PCL <50), and "severe" PTSD (PCL ≥50). Subtypes of PTSD by symptom clusters were analyzed. Multinomial logistic regression was used to estimate the association of PTSD with two GFR change outcomes (decline or increase) compared with the stable GFR outcome. RESULTS: In 2266 participants, the mean age was 53.1 years, 8.2% were female, and 89.1% were White. Individuals with PTSD (n = 373; 16.5%) did not differ in mean baseline GFR from individuals without PTSD (89.73 versus 90.56 mL min-1 1.73 m-2; p = .29). During a 2.01-year mean follow-up, a mean GFR decline of -1.51 mL min-1 1.73 m-2 per year was noted. In multivariable-adjusted models, PTSD was associated with GFR decline (adjusted relative risk [aRR] = 1.74 [1.32-2.30], p < .001) compared with stable GFR, with "hyperarousal" symptoms showing the strongest association (aRR =2.11 [1.40-3.19]; p < .001). Dose-response effects were evident when comparing mild with severe PTSD and comparing PTSD with versus without depression. PTSD was also associated with GFR rise (aRR = 1.47 [1.10-1.97], p < .009). The association between PTSD and GFR change was stronger in participants older than 50 years. CONCLUSIONS: PTSD may be a novel risk factor for exaggerated longitudinal GFR change in young, healthy adults. These findings need to be validated in other cohorts.

Negative binomial additive model for RNA-Seq data analysis.

BACKGROUND: High-throughput sequencing experiments followed by differential expression analysis is a widely used approach for detecting genomic biomarkers. A fundamental step in differential expression analysis is to model the association between gene counts and covariates of interest. Existing models assume linear effect of covariates, which is restrictive and may not be sufficient for certain phenotypes. RESULTS: We introduce NBAMSeq, a flexible statistical model based on the generalized additive model and allows for information sharing across genes in variance estimation. Specifically, we model the logarithm of mean gene counts as sums of smooth functions with the smoothing parameters and coefficients estimated simultaneously within a nested iterative method. The variance is estimated by the Bayesian shrinkage approach to fully exploit the information across all genes. CONCLUSIONS: Based on extensive simulations and case studies of RNA-Seq data, we show that NBAMSeq offers improved performance in detecting nonlinear effect and maintains equivalent performance in detecting linear effect compared to existing methods. The vignette and source code of NBAMSeq are available at http://bioconductor.org/packages/release/bioc/html/NBAMSeq.html.

methylGSA: a Bioconductor package and Shiny app for DNA methylation data length bias adjustment in gene set testing.

MOTIVATION: An important downstream analysis following differential expression from RNA sequencing (RNA-Seq) or DNA methylation analysis is the gene set testing to relate significant genes or CpGs to known biological properties. However, the traditional gene set testing approaches result in biased P-values due to the difference in gene length. Existing methods accounting for length bias were primarily developed for RNA-Seq data. For DNA methylation data profiled using the Illumina arrays, separate methods adjusting for the number of CpGs instead of gene length are necessary. RESULTS: We developed methylGSA, a Bioconductor package for gene set testing in DNA methylation data. Our accompanying Shiny app provides an interactive way of accessing functions and visualizing the results in methylGSA package. AVAILABILITY AND IMPLEMENTATION: methylGSA is available at Bioconductor repository: https://bioconductor.org/packages/methylGSA and Shiny app is available at: http://www.ams.sunysb.edu/%7epfkuan/softwares.html#methylGSA. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Polygenic prediction of PTSD trajectories in 9/11 responders.

BACKGROUND: Genetics hold promise of predicting long-term post-traumatic stress disorder (PTSD) outcomes following trauma. The aim of the current study was to test whether six hypothesized polygenic risk scores (PRSs) developed to capture genetic vulnerability to psychiatric conditions prospectively predict PTSD onset, severity, and 18-year course after trauma exposure. METHODS: Participants were 1490 responders to the World Trade Center (WTC) disaster (mean age at 9/11 = 38.81 years, s.d. = 8.20; 93.5% male; 23.8% lifetime WTC-related PTSD diagnosis). Prospective longitudinal data on WTC-related PTSD symptoms were obtained from electronic medical records and modelled as PTSD trajectories using growth mixture model analysis. Independent regression models tested whether six hypothesized psychiatric PRSs (PTSD-PRS, Re-experiencing-PRS, Generalized Anxiety-PRS, Schizophrenia-PRS, Depression-PRS, and Neuroticism-PRS) are predictive of WTC-PTSD outcomes: lifetime diagnoses, average symptom severity, and 18-year symptom trajectory. All analyses were adjusted for population stratification, 9/11 exposure severity, and multiple testing. RESULTS: Depression-PRS predicted PTSD diagnostic status (OR 1.37, CI 1.17-1.61, adjusted p = 0.001). All PRSs, except PTSD-PRS, significantly predicted average PTSD symptoms (ß = 0.06-0.10, adjusted p < 0.05). Re-experiencing-PRS, Generalized Anxiety-PRS and Schizophrenia-PRS predicted the high severity PTSD trajectory class (ORs 1.21-1.28, adjusted p < 0.05). Finally, PRSs prediction was independent of 9/11 exposure severity and jointly accounted for 3.7 times more variance in PTSD symptoms than the exposure severity. CONCLUSIONS: Psychiatric PRSs prospectively predicted WTC-related PTSD lifetime diagnosis, average symptom severity, and 18-year trajectory in responders to 9/11 disaster. Jointly, PRSs were more predictive of subsequent PTSD than the exposure severity. In the future, PRSs may help identify at-risk responders who might benefit from targeted prevention approaches.

Multiple-kernel learning for genomic data mining and prediction.

BACKGROUND: Advances in medical technology have allowed for customized prognosis, diagnosis, and treatment regimens that utilize multiple heterogeneous data sources. Multiple kernel learning (MKL) is well suited for the integration of multiple high throughput data sources. MKL remains to be under-utilized by genomic researchers partly due to the lack of unified guidelines for its use, and benchmark genomic datasets. RESULTS: We provide three implementations of MKL in R. These methods are applied to simulated data to illustrate that MKL can select appropriate models. We also apply MKL to combine clinical information with miRNA gene expression data of ovarian cancer study into a single analysis. Lastly, we show that MKL can identify gene sets that are known to play a role in the prognostic prediction of 15 cancer types using gene expression data from The Cancer Genome Atlas, as well as, identify new gene sets for the future research. CONCLUSION: Multiple kernel learning coupled with modern optimization techniques provides a promising learning tool for building predictive models based on multi-source genomic data. MKL also provides an automated scheme for kernel prioritization and parameter tuning. The methods used in the paper are implemented as an R package called RMKL package, which is freely available for download through CRAN at https://CRAN.R-project.org/package=RMKL .

An aberrant NOTCH2-BCR signaling axis in B cells from patients with chronic GVHD.

B-cell receptor (BCR)-activated B cells contribute to pathogenesis in chronic graft-versus-host disease (cGVHD), a condition manifested by both B-cell autoreactivity and immune deficiency. We hypothesized that constitutive BCR activation precluded functional B-cell maturation in cGVHD. To address this, we examined BCR-NOTCH2 synergy because NOTCH has been shown to increase BCR responsiveness in normal mouse B cells. We conducted ex vivo activation and signaling assays of 30 primary samples from hematopoietic stem cell transplantation patients with and without cGVHD. Consistent with a molecular link between pathways, we found that BCR-NOTCH activation significantly increased the proximal BCR adapter protein BLNK. BCR-NOTCH activation also enabled persistent NOTCH2 surface expression, suggesting a positive feedback loop. Specific NOTCH2 blockade eliminated NOTCH-BCR activation and significantly altered NOTCH downstream targets and B-cell maturation/effector molecules. Examination of the molecular underpinnings of this "NOTCH2-BCR axis" in cGVHD revealed imbalanced expression of the transcription factors IRF4 and IRF8, each critical to B-cell differentiation and fate. All-trans retinoic acid (ATRA) increased IRF4 expression, restored the IRF4-to-IRF8 ratio, abrogated BCR-NOTCH hyperactivation, and reduced NOTCH2 expression in cGVHD B cells without compromising viability. ATRA-treated cGVHD B cells had elevated TLR9 and PAX5, but not BLIMP1 (a gene-expression pattern associated with mature follicular B cells) and also attained increased cytosine guanine dinucleotide responsiveness. Together, we reveal a mechanistic link between NOTCH2 activation and robust BCR responses to otherwise suboptimal amounts of surrogate antigen. Our findings suggest that peripheral B cells in cGVHD patients can be pharmacologically directed from hyperactivation toward maturity.

Quantitative decision-making in randomized Phase II studies with a time-to-event endpoint.

One of the most critical decision points in clinical development is Go/No-Go decision-making after a proof-of-concept study. Traditional decision-making relies on a formal hypothesis testing with control of type I and type II error rates, which is limited by assessing the strength of efficacy evidence in a small isolated trial. In this article, we propose a quantitative Bayesian/frequentist decision framework for Go/No-Go criteria and sample size evaluation in Phase II randomized studies with a time-to-event endpoint. By taking the uncertainty of treatment effect into consideration, we propose an integrated quantitative approach for a program when both the Phase II and Phase III trials share a common endpoint while allowing a discount of the observed Phase II data. Our results confirm the argument that an increase in the sample size of a Phase II trial will result in greater increase in the probability of success of a Phase III trial than increasing the Phase III trial sample size by equal amount. We illustrate the steps in quantitative decision-making with a real example of a randomized Phase II study in metastatic pancreatic cancer.

Utility of TERT Promoter Mutations for Cutaneous Primary Melanoma Diagnosis.

Telomerase reverse transcriptase (TERT) promoter mutations are commonly found in malignant melanomas but rare in melanocytic nevi. To assess its potential diagnostic utility for the distinction of melanoma from nevus, we determined the TERT promoter mutation status of 86 primary melanomas, 72 melanocytic nevi, and 40 diagnostically problematic melanocytic proliferations. Of the 86 melanomas, 67 (77.9%) were TERT-positive, defined as harboring a hotspot TERT promoter mutation at positions -124C>T, -124_125CC>TT, -138_139CC>TT, or -146C>T. Of the 72 nevi, only 1 (1.4%) was TERT-positive. Of the 40 diagnostically uncertain melanocytic proliferations, 2 (5.0%) were TERT-positive. TERT positivity as a test for melanoma versus nevus had an accuracy of 87.3% [95% confidence interval (CI), 81.1-92.1], a sensitivity of 77.9% (95% CI, 68.9-85.4), a specificity of 98.6% (95% CI, 95.8-100), a positive predictive value of 98.5% (95% CI, 95.6-100), and a negative predictive value of 78.9% (95% CI, 72.6-85.4). Our results indicate that hotspot TERT promoter mutation status may be a useful ancillary parameter for the diagnosis of melanoma. In particular, the high specificity of these mutations for melanoma indicates the presence of a TERT promoter mutation in a melanocytic neoplasm associated with diagnostic controversy, or uncertainty should increase concern for a melanoma.

Comparison of the restricted mean survival time with the hazard ratio in superiority trials with a time-to-event end point.

With the emergence of novel therapies exhibiting distinct mechanisms of action compared to traditional treatments, departure from the proportional hazard (PH) assumption in clinical trials with a time-to-event end point is increasingly common. In these situations, the hazard ratio may not be a valid statistical measurement of treatment effect, and the log-rank test may no longer be the most powerful statistical test. The restricted mean survival time (RMST) is an alternative robust and clinically interpretable summary measure that does not rely on the PH assumption. We conduct extensive simulations to evaluate the performance and operating characteristics of the RMST-based inference and against the hazard ratio-based inference, under various scenarios and design parameter setups. The log-rank test is generally a powerful test when there is evident separation favoring 1 treatment arm at most of the time points across the Kaplan-Meier survival curves, but the performance of the RMST test is similar. Under non-PH scenarios where late separation of survival curves is observed, the RMST-based test has better performance than the log-rank test when the truncation time is reasonably close to the tail of the observed curves. Furthermore, when flat survival tail (or low event rate) in the experimental arm is expected, selecting the minimum of the maximum observed event time as the truncation timepoint for the RMST is not recommended. In addition, we recommend the inclusion of analysis based on the RMST curve over the truncation time in clinical settings where there is suspicion of substantial departure from the PH assumption.

A systematic evaluation of nucleotide properties for CRISPR sgRNA design.

BACKGROUND: CRISPR is a versatile gene editing tool which has revolutionized genetic research in the past few years. Optimizing sgRNA design to improve the efficiency of target/DNA cleavage is critical to ensure the success of CRISPR screens. RESULTS: By borrowing knowledge from oligonucleotide design and nucleosome occupancy models, we systematically evaluated candidate features computed from a number of nucleic acid, thermodynamic and secondary structure models on real CRISPR datasets. Our results showed that taking into account position-dependent dinucleotide features improved the design of effective sgRNAs with area under the receiver operating characteristic curve (AUC) >0.8, and the inclusion of additional features offered marginal improvement (∼2% increase in AUC). CONCLUSION: Using a machine-learning approach, we proposed an accurate prediction model for sgRNA design efficiency. An R package predictSGRNA implementing the predictive model is available at http://www.ams.sunysb.edu/~pfkuan/softwares.html#predictsgrna .

Variation in chromatin accessibility in human kidney cancer links H3K36 methyltransferase loss with widespread RNA processing defects.

Comprehensive sequencing of human cancers has identified recurrent mutations in genes encoding chromatin regulatory proteins. For clear cell renal cell carcinoma (ccRCC), three of the five commonly mutated genes encode the chromatin regulators PBRM1, SETD2, and BAP1. How these mutations alter the chromatin landscape and transcriptional program in ccRCC or other cancers is not understood. Here, we identified alterations in chromatin organization and transcript profiles associated with mutations in chromatin regulators in a large cohort of primary human kidney tumors. By associating variation in chromatin organization with mutations in SETD2, which encodes the enzyme responsible for H3K36 trimethylation, we found that changes in chromatin accessibility occurred primarily within actively transcribed genes. This increase in chromatin accessibility was linked with widespread alterations in RNA processing, including intron retention and aberrant splicing, affecting ∼25% of all expressed genes. Furthermore, decreased nucleosome occupancy proximal to misspliced exons was observed in tumors lacking H3K36me3. These results directly link mutations in SETD2 to chromatin accessibility changes and RNA processing defects in cancer. Detecting the functional consequences of specific mutations in chromatin regulatory proteins in primary human samples could ultimately inform the therapeutic application of an emerging class of chromatin-targeted compounds.

Breast tumor DNA methylation patterns associated with smoking in the Carolina Breast Cancer Study.

PURPOSE: Tobacco smoking is a risk factor in several cancers, yet its roles as a putative etiologic exposure or poor prognostic factor in breast cancer are less clear. Altered DNA methylation contributes to breast cancer development and may provide a mechanistic link between smoking and gene expression changes leading to cancer development or progression. METHODS: Using a cancer-focused array, we examined methylation at 933 CpGs in 517 invasive breast tumors in the Carolina Breast Cancer Study to determine whether methylation patterns differ by exposure to tobacco smoke. Multivariable generalized linear regression models were used to compare tumor methylation profiles between smokers and never smokers, overall, or stratified on hormone receptor (HR) status. RESULTS: Modest differences in CpG methylation were detected at p < 0.05 in breast tumors from current or ever smokers compared with never smokers. In stratified analyses, HR- tumors from smokers exhibited primarily hypomethylation compared with tumors from never smokers; hypomethylation was similarly detected within the more homogeneous basal-like subtype. Most current smoking-associated CpG loci exhibited methylation levels in former smokers that were intermediate between those in current and never smokers and exhibited progressive changes in methylation with increasing duration of smoking. Among former smokers, restoration of methylation toward baseline (never smoking) levels was observed with increasing time since quitting. Moreover, smoking-related hypermethylation was stronger in HR+ breast tumors from blacks than in whites. CONCLUSIONS: Our results suggest that breast tumor methylation patterns differ with tobacco smoke exposure; however, additional studies are needed to confirm these findings.

Epigenome-wide association of PTSD from heterogeneous cohorts with a common multi-site analysis pipeline.

Compelling evidence suggests that epigenetic mechanisms such as DNA methylation play a role in stress regulation and in the etiologic basis of stress related disorders such as Post traumatic Stress Disorder (PTSD). Here we describe the purpose and methods of an international consortium that was developed to study the role of epigenetics in PTSD. Inspired by the approach used in the Psychiatric Genomics Consortium, we brought together investigators representing seven cohorts with a collective sample size of N = 1147 that included detailed information on trauma exposure, PTSD symptoms, and genome-wide DNA methylation data. The objective of this consortium is to increase the analytical sample size by pooling data and combining expertise so that DNA methylation patterns associated with PTSD can be identified. Several quality control and analytical pipelines were evaluated for their control of genomic inflation and technical artifacts with a joint analysis procedure established to derive comparable data over the cohorts for meta-analysis. We propose methods to deal with ancestry population stratification and type I error inflation and discuss the advantages and disadvantages of applying robust error estimates. To evaluate our pipeline, we report results from an epigenome-wide association study (EWAS) of age, which is a well-characterized phenotype with known epigenetic associations. Overall, while EWAS are highly complex and subject to similar challenges as genome-wide association studies (GWAS), we demonstrate that an epigenetic meta-analysis with a relatively modest sample size can be well-powered to identify epigenetic associations. Our pipeline can be used as a framework for consortium efforts for EWAS.

Covariate adjusted differential variability analysis of DNA methylation with propensity score method.

It has been proposed recently that differentially variable CpG methylation (DVC) may contribute to transcriptional aberrations in human diseases. In large scale epigenetic studies, potential confounders could affect the observed methylation variabilities and need to be accounted for. In this paper, we develop a robust statistical model for differential variability DVC analysis that accounts for potential confounding covariates by utilizing the propensity score method. Our method is based on a weighted score test on strata generated propensity score stratification. To the best of our knowledge, this is the first proposed statistical method for detecting DVCs that adjusts for confounding covariates. We show that this method is robust against model misspecification and achieves good operating characteristics based on extensive simulations and a case study.