Exploring the role of neuronal-enriched extracellular vesicle miR-93 and interoception in major depressive disorder.

Major depressive disorder (MDD) is associated with interoceptive processing dysfunctions, but the molecular mechanisms underlying this dysfunction are poorly understood. This study combined brain neuronal-enriched extracellular vesicle (NEEV) technology and serum markers of inflammation and metabolism with Functional Magnetic Resonance Imaging (fMRI) to identify the contribution of gene regulatory pathways, in particular micro-RNA (miR) 93, to interoceptive dysfunction in MDD. Individuals with MDD (n = 41) and healthy comparisons (HC; n = 35) provided blood samples and completed an interoceptive attention task during fMRI. EVs were separated from plasma using a precipitation method. NEEVs were enriched by magnetic streptavidin bead immunocapture utilizing a neural adhesion marker (L1CAM/CD171) biotinylated antibody. The origin of NEEVs was validated with two other neuronal markers - neuronal cell adhesion molecule (NCAM) and ATPase Na+/K+ transporting subunit alpha 3 (ATP1A3). NEEV specificities were confirmed by flow cytometry, western blot, particle size analyzer, and transmission electron microscopy. NEEV small RNAs were purified and sequenced. Results showed that: (1) MDD exhibited lower NEEV miR-93 expression than HC; (2) within MDD but not HC, those individuals with the lowest NEEV miR-93 expression had the highest serum concentrations of interleukin (IL)-1 receptor antagonist, IL-6, tumor necrosis factor, and leptin; and (3) within HC but not MDD, those participants with the highest miR-93 expression showed the strongest bilateral dorsal mid-insula activation during interoceptive versus exteroceptive attention. Since miR-93 is regulated by stress and affects epigenetic modulation by chromatin re-organization, these results suggest that healthy individuals but not MDD participants show an adaptive epigenetic regulation of insular function during interoceptive processing. Future investigations will need to delineate how specific internal and external environmental conditions contribute to miR-93 expression in MDD and what molecular mechanisms alter brain responsivity to body-relevant signals.

Differential expression of MDGA1 in major depressive disorder.

The identification of gene expression-based biomarkers for major depressive disorder (MDD) continues to be an important challenge. In order to identify candidate biomarkers and mechanisms, we apply statistical and machine learning feature selection to an RNA-Seq gene expression dataset of 78 unmedicated individuals with MDD and 79 healthy controls. We identify 49 genes by LASSO penalized logistic regression and 45 genes at the false discovery rate threshold 0.188. The MDGA1 gene has the lowest P-value (4.9e-5) and is expressed in the developing brain, involved in axon guidance, and associated with related mood disorders in previous studies of bipolar disorder (BD) and schizophrenia (SCZ). The expression of MDGA1 is associated with age and sex, but its association with MDD remains significant when adjusted for covariates. MDGA1 is in a co-expression cluster with another top gene, ATXN7L2 (ataxin 7 like 2), which was associated with MDD in a recent GWAS. The LASSO classification model of MDD includes MDGA1, and the model has a cross-validation accuracy of 79%. Another noteworthy top gene, IRF2BPL, is in a close co-expression cluster with MDGA1 and may be related to microglial inflammatory states in MDD. Future exploration of MDGA1 and its gene interactions may provide insights into mechanisms and heterogeneity of MDD.

Common Data Elements, Scalable Data Management Infrastructure, and Analytics Workflows for Large-Scale Neuroimaging Studies.

Neuroscience studies require considerable bioinformatic support and expertise. Numerous high-dimensional and multimodal datasets must be preprocessed and integrated to create robust and reproducible analysis pipelines. We describe a common data elements and scalable data management infrastructure that allows multiple analytics workflows to facilitate preprocessing, analysis and sharing of large-scale multi-level data. The process uses the Brain Imaging Data Structure (BIDS) format and supports MRI, fMRI, EEG, clinical, and laboratory data. The infrastructure provides support for other datasets such as Fitbit and flexibility for developers to customize the integration of new types of data. Exemplar results from 200+ participants and 11 different pipelines demonstrate the utility of the infrastructure.

Theoretical properties of distance distributions and novel metrics for nearest-neighbor feature selection.

The performance of nearest-neighbor feature selection and prediction methods depends on the metric for computing neighborhoods and the distribution properties of the underlying data. Recent work to improve nearest-neighbor feature selection algorithms has focused on new neighborhood estimation methods and distance metrics. However, little attention has been given to the distributional properties of pairwise distances as a function of the metric or data type. Thus, we derive general analytical expressions for the mean and variance of pairwise distances for Lq metrics for normal and uniform random data with p attributes and m instances. The distribution moment formulas and detailed derivations provide a resource for understanding the distance properties for metrics and data types commonly used with nearest-neighbor methods, and the derivations provide the starting point for the following novel results. We use extreme value theory to derive the mean and variance for metrics that are normalized by the range of each attribute (difference of max and min). We derive analytical formulas for a new metric for genetic variants, which are categorical variables that occur in genome-wide association studies (GWAS). The genetic distance distributions account for minor allele frequency and the transition/transversion ratio. We introduce a new metric for resting-state functional MRI data (rs-fMRI) and derive its distance distribution properties. This metric is applicable to correlation-based predictors derived from time-series data. The analytical means and variances are in strong agreement with simulation results. We also use simulations to explore the sensitivity of the expected means and variances in the presence of correlation and interactions in the data. These analytical results and new metrics can be used to inform the optimization of nearest neighbor methods for a broad range of studies, including gene expression, GWAS, and fMRI data.

Nearest-Neighbor Projected Distance Regression for Epistasis Detection in GWAS With Population Structure Correction.

Nearest-neighbor Projected-Distance Regression (NPDR) is a feature selection technique that uses nearest-neighbors in high dimensional data to detect complex multivariate effects including epistasis. NPDR uses a regression formalism that allows statistical significance testing and efficient control for multiple testing. In addition, the regression formalism provides a mechanism for NPDR to adjust for population structure, which we apply to a GWAS of systemic lupus erythematosus (SLE). We also test NPDR on benchmark simulated genetic variant data with epistatic effects, main effects, imbalanced data for case-control design and continuous outcomes. NPDR identifies potential interactions in an epistasis network that influences the SLE disorder.

Correction: Identification and replication of RNA-Seq gene network modules associated with depression severity.

An amendment to this paper has been published and can be accessed via a link at the top of the paper.

Consensus features nested cross-validation.

SUMMARY: Feature selection can improve the accuracy of machine-learning models, but appropriate steps must be taken to avoid overfitting. Nested cross-validation (nCV) is a common approach that chooses the classification model and features to represent a given outer fold based on features that give the maximum inner-fold accuracy. Differential privacy is a related technique to avoid overfitting that uses a privacy-preserving noise mechanism to identify features that are stable between training and holdout sets.We develop consensus nested cross-validation (cnCV) that combines the idea of feature stability from differential privacy with nCV. Feature selection is applied in each inner fold and the consensus of top features across folds is used as a measure of feature stability or reliability instead of classification accuracy, which is used in standard nCV. We use simulated data with main effects, correlation and interactions to compare the classification accuracy and feature selection performance of the new cnCV with standard nCV, Elastic Net optimized by cross-validation, differential privacy and private evaporative cooling (pEC). We also compare these methods using real RNA-seq data from a study of major depressive disorder.The cnCV method has similar training and validation accuracy to nCV, but cnCV has much shorter run times because it does not construct classifiers in the inner folds. The cnCV method chooses a more parsimonious set of features with fewer false positives than nCV. The cnCV method has similar accuracy to pEC and cnCV selects stable features between folds without the need to specify a privacy threshold. We show that cnCV is an effective and efficient approach for combining feature selection with classification. AVAILABILITY AND IMPLEMENTATION: Code available at https://github.com/insilico/cncv. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Correction: Tissue-Specific Expressed Antibody Variable Gene Repertoires.

[This corrects the article DOI: 10.1371/journal.pone.0100839.].

Nearest-neighbor Projected-Distance Regression (NPDR) for detecting network interactions with adjustments for multiple tests and confounding.

SUMMARY: Machine learning feature selection methods are needed to detect complex interaction-network effects in complicated modeling scenarios in high-dimensional data, such as GWAS, gene expression, eQTL and structural/functional neuroimage studies for case-control or continuous outcomes. In addition, many machine learning methods have limited ability to address the issues of controlling false discoveries and adjusting for covariates. To address these challenges, we develop a new feature selection technique called Nearest-neighbor Projected-Distance Regression (NPDR) that calculates the importance of each predictor using generalized linear model regression of distances between nearest-neighbor pairs projected onto the predictor dimension. NPDR captures the underlying interaction structure of data using nearest-neighbors in high dimensions, handles both dichotomous and continuous outcomes and predictor data types, statistically corrects for covariates, and permits statistical inference and penalized regression. We use realistic simulations with interactions and other effects to show that NPDR has better precision-recall than standard Relief-based feature selection and random forest importance, with the additional benefit of covariate adjustment and multiple testing correction. Using RNA-Seq data from a study of major depressive disorder (MDD), we show that NPDR with covariate adjustment removes spurious associations due to confounding. We apply NPDR to eQTL data to identify potentially interacting variants that regulate transcripts associated with MDD and demonstrate NPDR's utility for GWAS and continuous outcomes. AVAILABILITY AND IMPLEMENTATION: Available at: https://insilico.github.io/npdr/. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

A Humanized Yeast Phenomic Model of Deoxycytidine Kinase to Predict Genetic Buffering of Nucleoside Analog Cytotoxicity.

Knowledge about synthetic lethality can be applied to enhance the efficacy of anticancer therapies in individual patients harboring genetic alterations in their cancer that specifically render it vulnerable. We investigated the potential for high-resolution phenomic analysis in yeast to predict such genetic vulnerabilities by systematic, comprehensive, and quantitative assessment of drug-gene interaction for gemcitabine and cytarabine, substrates of deoxycytidine kinase that have similar molecular structures yet distinct antitumor efficacy. Human deoxycytidine kinase (dCK) was conditionally expressed in the Saccharomycescerevisiae genomic library of knockout and knockdown (YKO/KD) strains, to globally and quantitatively characterize differential drug-gene interaction for gemcitabine and cytarabine. Pathway enrichment analysis revealed that autophagy, histone modification, chromatin remodeling, and apoptosis-related processes influence gemcitabine specifically, while drug-gene interaction specific to cytarabine was less enriched in gene ontology. Processes having influence over both drugs were DNA repair and integrity checkpoints and vesicle transport and fusion. Non-gene ontology (GO)-enriched genes were also informative. Yeast phenomic and cancer cell line pharmacogenomics data were integrated to identify yeast-human homologs with correlated differential gene expression and drug efficacy, thus providing a unique resource to predict whether differential gene expression observed in cancer genetic profiles are causal in tumor-specific responses to cytotoxic agents.

Multi-Level Model to Predict Antibody Response to Influenza Vaccine Using Gene Expression Interaction Network Feature Selection.

Vaccination is an effective prevention of influenza infection. However, certain individuals develop a lower antibody response after vaccination, which may lead to susceptibility to subsequent infection. An important challenge in human health is to find baseline gene signatures to help identify individuals who are at higher risk for infection despite influenza vaccination. We developed a multi-level machine learning strategy to build a predictive model of vaccine response using pre-vaccination antibody titers and network interactions between pre-vaccination gene expression levels. The first-level baseline-antibody model explains a significant amount of variation in post-vaccination response, especially for subjects with large pre-existing antibody titers. In the second level, we clustered individuals based on pre-vaccination antibody titers to focus gene-based modeling on individuals with lower baseline HAI where additional response variation may be predicted by baseline gene expression levels. In the third level, we used a gene-association interaction network (GAIN) feature selection algorithm to find the best pairs of genes that interact to influence antibody response within each baseline titer cluster. We used ratios of the top interacting genes as predictors to stabilize machine learning model generalizability. We trained and tested the multi-level approach on data with young and older individuals immunized against influenza vaccine in multiple cohorts. Our results indicate that the GAIN feature selection approach improves model generalizability and identifies genes enriched for immunologically relevant pathways, including B Cell Receptor signaling and antigen processing. Using a multi-level approach, starting with a baseline HAI model and stratifying on baseline HAI, allows for more targeted gene-based modeling. We provide an interactive tool that may be extended to other vaccine studies.

The role of systems biology approaches in determining molecular signatures for the development of more effective vaccines.

INTRODUCTION: Emerging infectious diseases are a major threat to public health, and while vaccines have proven to be one of the most effective preventive measures for infectious diseases, we still do not have safe and effective vaccines against many human pathogens, and emerging diseases continually pose new threats. The purpose of this review is to discuss how the creation of vaccines for these new threats has been hindered by limitations in the current approach to vaccine development. Recent advances in high-throughput technologies have enabled scientists to apply systems biology approaches to collect and integrate increasingly large datasets that capture comprehensive biological changes induced by vaccines, and then decipher the complex immune response to those vaccines. AREAS COVERED: This review covers advances in these technologies and recent publications that describe systems biology approaches to understanding vaccine immune responses and to understanding the rational design of new vaccine candidates. EXPERT OPINION: Systems biology approaches to vaccine development provide novel information regarding both the immune response and the underlying mechanisms and can inform vaccine development.

EpistasisRank and EpistasisKatz: interaction network centrality methods that integrate prior knowledge networks.

MOTIVATION: An important challenge in gene expression analysis is to improve hub gene selection to enrich for biological relevance or improve classification accuracy for a given phenotype. In order to incorporate phenotypic context into co-expression, we recently developed an epistasis-expression network centrality method that blends the importance of gene-gene interactions (epistasis) and main effects of genes. Further blending of prior knowledge from functional interactions has the potential to enrich for relevant genes and stabilize classification. RESULTS: We develop two new expression-epistasis centrality methods that incorporate interaction prior knowledge. The first extends our SNPrank (EpistasisRank) method by incorporating a gene-wise prior knowledge vector. This prior knowledge vector informs the centrality algorithm of the inclination of a gene to be involved in interactions by incorporating functional interaction information from the Integrative Multi-species Prediction database. The second method extends Katz centrality to expression-epistasis networks (EpistasisKatz), extends the Katz bias to be a gene-wise vector of main effects and extends the Katz attenuation constant prefactor to be a prior-knowledge vector for interactions. Using independent microarray studies of major depressive disorder, we find that including prior knowledge in network centrality feature selection stabilizes the training classification and reduces over-fitting. AVAILABILITY AND IMPLEMENTATION: Methods and examples provided at https://github.com/insilico/Rinbix and https://github.com/insilico/PriorKnowledgeEpistasisRank. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

STatistical Inference Relief (STIR) feature selection.

MOTIVATION: Relief is a family of machine learning algorithms that uses nearest-neighbors to select features whose association with an outcome may be due to epistasis or statistical interactions with other features in high-dimensional data. Relief-based estimators are non-parametric in the statistical sense that they do not have a parameterized model with an underlying probability distribution for the estimator, making it difficult to determine the statistical significance of Relief-based attribute estimates. Thus, a statistical inferential formalism is needed to avoid imposing arbitrary thresholds to select the most important features. We reconceptualize the Relief-based feature selection algorithm to create a new family of STatistical Inference Relief (STIR) estimators that retains the ability to identify interactions while incorporating sample variance of the nearest neighbor distances into the attribute importance estimation. This variance permits the calculation of statistical significance of features and adjustment for multiple testing of Relief-based scores. Specifically, we develop a pseudo t-test version of Relief-based algorithms for case-control data. RESULTS: We demonstrate the statistical power and control of type I error of the STIR family of feature selection methods on a panel of simulated data that exhibits properties reflected in real gene expression data, including main effects and network interaction effects. We compare the performance of STIR when the adaptive radius method is used as the nearest neighbor constructor with STIR when the fixed-k nearest neighbor constructor is used. We apply STIR to real RNA-Seq data from a study of major depressive disorder and discuss STIR's straightforward extension to genome-wide association studies. AVAILABILITY AND IMPLEMENTATION: Code and data available at http://insilico.utulsa.edu/software/STIR. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

A Nonlinear Simulation Framework Supports Adjusting for Age When Analyzing BrainAGE.

Several imaging modalities, including T1-weighted structural imaging, diffusion tensor imaging, and functional MRI can show chronological age related changes. Employing machine learning algorithms, an individual's imaging data can predict their age with reasonable accuracy. While details vary according to modality, the general strategy is to: (1) extract image-related features, (2) build a model on a training set that uses those features to predict an individual's age, (3) validate the model on a test dataset, producing a predicted age for each individual, (4) define the "Brain Age Gap Estimate" (BrainAGE) as the difference between an individual's predicted age and his/her chronological age, (5) estimate the relationship between BrainAGE and other variables of interest, and (6) make inferences about those variables and accelerated or delayed brain aging. For example, a group of individuals with overall positive BrainAGE may show signs of accelerated aging in other variables as well. There is inevitably an overestimation of the age of younger individuals and an underestimation of the age of older individuals due to "regression to the mean." The correlation between chronological age and BrainAGE may significantly impact the relationship between BrainAGE and other variables of interest when they are also related to age. In this study, we examine the detectability of variable effects under different assumptions. We use empirical results from two separate datasets [training = 475 healthy volunteers, aged 18-60 years (259 female); testing = 489 participants including people with mood/anxiety, substance use, eating disorders and healthy controls, aged 18-56 years (312 female)] to inform simulation parameter selection. Outcomes in simulated and empirical data strongly support the proposal that models incorporating BrainAGE should include chronological age as a covariate. We propose either including age as a covariate in step 5 of the above framework, or employing a multistep procedure where age is regressed on BrainAGE prior to step 5, producing BrainAGE Residualized (BrainAGER) scores.

Identification and replication of RNA-Seq gene network modules associated with depression severity.

Genomic variation underlying major depressive disorder (MDD) likely involves the interaction and regulation of multiple genes in a network. Data-driven co-expression network module inference has the potential to account for variation within regulatory networks, reduce the dimensionality of RNA-Seq data, and detect significant gene-expression modules associated with depression severity. We performed an RNA-Seq gene co-expression network analysis of mRNA data obtained from the peripheral blood mononuclear cells of unmedicated MDD (n = 78) and healthy control (n = 79) subjects. Across the combined MDD and HC groups, we assigned genes into modules using hierarchical clustering with a dynamic tree cut method and projected the expression data onto a lower-dimensional module space by computing the single-sample gene set enrichment score of each module. We tested the single-sample scores of each module for association with levels of depression severity measured by the Montgomery-Åsberg Depression Scale (MADRS). Independent of MDD status, we identified 23 gene modules from the co-expression network. Two modules were significantly associated with the MADRS score after multiple comparison adjustment (adjusted p = 0.009, 0.028 at 0.05 FDR threshold), and one of these modules replicated in a previous RNA-Seq study of MDD (p = 0.03). The two MADRS-associated modules contain genes previously implicated in mood disorders and show enrichment of apoptosis and B cell receptor signaling. The genes in these modules show a correlation between network centrality and univariate association with depression, suggesting that intramodular hub genes are more likely to be related to MDD compared to other genes in a module.

Altered populations of natural killer cells, cytotoxic T lymphocytes, and regulatory T cells in major depressive disorder: Association with sleep disturbance.

A subset of individuals with major depressive disorder (MDD) have impaired adaptive immunity characterized by a greater vulnerability to viral infection and a deficient response to vaccination along with a decrease in the number and/or activity of T cells and natural killer cells (NKC). Nevertheless, it remains unclear which specific subsets of lymphocytes are altered in MDD, a shortcoming we address here by utilizing an advanced fluorescence-activated cell sorting (FACS) method that allows for the differentiation of important functionally-distinct lymphocyte sub-populations. Furthermore, despite evidence that sleep disturbance, which is a core symptom of MDD, is itself associated with alterations in lymphocyte distributions, there is a paucity of studies examining the contribution of sleep disturbance on lymphocyte populations in MDD populations. Here, we measured differences in the percentages of 13 different lymphocytes and 6 different leukocytes in 54 unmedicated MDD patients (partially remitted to moderate) and 56 age and sex-matched healthy controls (HC). The relationship between self-reported sleep disturbance and cell counts was evaluated in the MDD group using the Pittsburgh Sleep Quality Index (PSQI). The MDD group showed a significantly increased percentage of CD127low/CCR4+ Treg cells, and memory Treg cells, as well as a reduction in CD56+CD16- (putative immunoregulatory) NKC counts, the latter, prior to correction for body mass index. There also was a trend for higher effector memory CD8+ cell counts in the MDD group versus the HC group. Further, within the MDD group, self-reported sleep disturbance was associated with an increased percentage of effector memory CD8+ cells but with a lower percentage of CD56+CD16- NKC. These results provide important new insights into the immune pathways involved in MDD, and provide novel evidence that MDD and associated sleep disturbance increase effector memory CD8+ and Treg pathways. Targeting sleep disturbance may have implications as a therapeutic strategy to normalize NKC and memory CD8+ cells in MDD.

Differential privacy-based evaporative cooling feature selection and classification with relief-F and random forests.

MOTIVATION: Classification of individuals into disease or clinical categories from high-dimensional biological data with low prediction error is an important challenge of statistical learning in bioinformatics. Feature selection can improve classification accuracy but must be incorporated carefully into cross-validation to avoid overfitting. Recently, feature selection methods based on differential privacy, such as differentially private random forests and reusable holdout sets, have been proposed. However, for domains such as bioinformatics, where the number of features is much larger than the number of observations p≫n , these differential privacy methods are susceptible to overfitting. METHODS: We introduce private Evaporative Cooling, a stochastic privacy-preserving machine learning algorithm that uses Relief-F for feature selection and random forest for privacy preserving classification that also prevents overfitting. We relate the privacy-preserving threshold mechanism to a thermodynamic Maxwell-Boltzmann distribution, where the temperature represents the privacy threshold. We use the thermal statistical physics concept of Evaporative Cooling of atomic gases to perform backward stepwise privacy-preserving feature selection. RESULTS: On simulated data with main effects and statistical interactions, we compare accuracies on holdout and validation sets for three privacy-preserving methods: the reusable holdout, reusable holdout with random forest, and private Evaporative Cooling, which uses Relief-F feature selection and random forest classification. In simulations where interactions exist between attributes, private Evaporative Cooling provides higher classification accuracy without overfitting based on an independent validation set. In simulations without interactions, thresholdout with random forest and private Evaporative Cooling give comparable accuracies. We also apply these privacy methods to human brain resting-state fMRI data from a study of major depressive disorder. AVAILABILITY AND IMPLEMENTATION: Code available at http://insilico.utulsa.edu/software/privateEC . CONTACT: brett-mckinney@utulsa.edu. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

The Integration of Epistasis Network and Functional Interactions in a GWAS Implicates RXR Pathway Genes in the Immune Response to Smallpox Vaccine.

Although many diseases and traits show large heritability, few genetic variants have been found to strongly separate phenotype groups by genotype. Complex regulatory networks of variants and expression of multiple genes lead to small individual-variant effects and difficulty replicating the effect of any single variant in an affected pathway. Interaction network modeling of GWAS identifies effects ignored by univariate models, but population differences may still cause specific genes to not replicate. Integrative network models may help detect indirect effects of variants in the underlying biological pathway. In this study, we used gene-level functional interaction information from the Integrative Multi-species Prediction (IMP) tool to reveal important genes associated with a complex phenotype through evidence from epistasis networks and pathway enrichment. We test this method for augmenting variant-based network analyses with functional interactions by applying it to a smallpox vaccine immune response GWAS. The integrative analysis spotlights the role of genes related to retinoid X receptor alpha (RXRA), which has been implicated in a previous epistasis network analysis of smallpox vaccine.

Recursive Indirect-Paths Modularity (RIP-M) for Detecting Community Structure in RNA-Seq Co-expression Networks.

Clusters of genes in co-expression networks are commonly used as functional units for gene set enrichment detection and increasingly as features (attribute construction) for statistical inference and sample classification. One of the practical challenges of clustering for these purposes is to identify an optimal partition of the network where the individual clusters are neither too large, prohibiting interpretation, nor too small, precluding general inference. Newman Modularity is a spectral clustering algorithm that automatically finds the number of clusters, but for many biological networks the cluster sizes are suboptimal. In this work, we generalize Newman Modularity to incorporate information from indirect paths in RNA-Seq co-expression networks. We implement a merge-and-split algorithm that allows the user to constrain the range of cluster sizes: large enough to capture genes in relevant pathways, yet small enough to resolve distinct functions. We investigate the properties of our recursive indirect-pathways modularity (RIP-M) and compare it with other clustering methods using simulated co-expression networks and RNA-seq data from an influenza vaccine response study. RIP-M had higher cluster assignment accuracy than Newman Modularity for finding clusters in simulated co-expression networks for all scenarios, and RIP-M had comparable accuracy to Weighted Gene Correlation Network Analysis (WGCNA). RIP-M was more accurate than WGCNA for modest hard thresholds and comparable for high, while WGCNA was slightly more accurate for soft thresholds. In the vaccine study data, RIP-M and WGCNA enriched for a comparable number of immunologically relevant pathways.