N-of-1-pathways MixEnrich: advancing precision medicine via single-subject analysis in discovering dynamic changes of transcriptomes

Background Transcriptome analytic tools are commonly used across patient cohorts to develop drugs and predict clinical outcomes. However, as precision medicine pursues more accurate and individualized treatment decisions, these methods are not designed to address single-patient transcriptome analyses. We previously developed and validated the N-of-1-pathways framework using two methods, Wilcoxon and Mahalanobis Distance (MD), for personal transcriptome analysis derived from a pair of samples of a single patient. Although, both methods uncover concordantly dysregulated pathways, they are not designed to detect dysregulated pathways with up- and down-regulated genes (bidirectional dysregulation) that are ubiquitous in biological systems. Results We developed N-of-1-pathways MixEnrich, a mixture model followed by a gene set enrichment test, to uncover bidirectional and concordantly dysregulated pathways one patient at a time. We assess its accuracy in a comprehensive simulation study and in a RNA-Seq data analysis of head and neck squamous cell carcinomas (HNSCCs). In presence of bidirectionally dysregulated genes in the pathway or in presence of high background noise, MixEnrich substantially outperforms previous single-subject transcriptome analysis methods, both in the simulation study and the HNSCCs data analysis (ROC Curves; higher true positive rates; lower false positive rates). Bidirectional and concordant dysregulated pathways uncovered by MixEnrich in each patient largely overlapped with the quasi-gold standard compared to other single-subject and cohort-based transcriptome analyses. Conclusion The greater performance of MixEnrich presents an advantage over previous methods to meet the promise of providing accurate personal transcriptome analysis to support precision medicine at point of care. Electronic supplementary material The online version of this article (doi:10.1186/s12920-017-0263-4) contains supplementary material, which is available to authorized users.


Background
Technologies, such as RNA-Seq, provide precise, timely, and cost-effective quantification of whole genome expression [1]. However, analytic tools remain underdeveloped for providing personal transcriptome profiling and individualized biological interpretation. Conventional transcriptome methods have been designed to uncover common mRNA and pathway signatures across a large cohort of patients, overlooking signals that differentiate one patient from another [2,3]. The analysis of dynamic transcriptomes of a single subject has the potential to capture and inform gene expression changes reflective of personal physiological modifications, disease progression, and response to therapies in ways that genetic information cannot. Indeed, the majority of disorders with complex inheritance results from a combination of genetic risks and environmental factors unique to each patient that dynamically influence the course of disease. These dynamic biological changes that are genome X environment interactions between two conditions can be measured at the transcriptome level; however, current cohort-based statistics, which average signals across patients, are not applicable for the analysis of personal transcriptome dynamics [4]. Although in vitro assays were used to assess dynamic gene expression changes to predict experimental outcomes and disease progression at the patient level, these analyses remain limited and biased as they only assess a handful of gene candidates pertaining to known pathways [5]. However, scaling-up these assays and analyses to measure whole genome expression changes of a single subject (e.g., before and after treatment) has the advantage to unbiasedly discover dysregulated pathways unique to each individual.
Recognizing the limitations of conventional methods, we recently designed and validated in different disease contexts the N-of-1-pathways, which is a novel framework for single-subject transcriptome analysis based on a pair of samples (e.g., healthy and tumor, before and after therapy) from the same individual [6][7][8][9][10]. N-of-1pathways relies on three principles: (1) the sole unit of observation is a single patient (case and control); (2) gene-level information are aggregated into gene sets (pathways); and (3) pathway results are summarized into personal biological profiling for clinical interpretation. Two methods under N-of-1-pathways framework were developed, N-of-1-pathways Wilcoxon (Wilcoxon) [6][7][8] using a Wilcoxon signed-rank test [11] and the N-of-1pathways Mahalanobis distance (MD) [10,12] using a statistical distance from a model of equal expression. The N-of-1-pathways Wilcoxon and MD analyze the dynamic change of mRNA expression and uncover dysregulated pathways (gene sets) from single-subject paired samples. The use of gene sets derived from gene ontology [13] provides computational advantage by reducing data dimension while providing mechanistic interpretation [14,15]. While both methods have shown promise in single-subject transcriptome analysis, they were not designed to identify pathways (gene sets) with both upregulated and down-regulated mRNA expressions and, therefore, take into account only concordantly dysregulated mRNAs within a pathway. In addition, Wilcoxon and MD are both self-contained methods [16] analyzing only mRNAs within a gene set and do not account for background noise due to technical and experimental artifacts [17][18][19].
To address the shortcomings of the current singlesubject transcriptome analysis methods, we developed a novel approach within the N-of-1-pathways framework: N-of-1-pathways MixEnrich (MixEnrich) using a mixture model (mixture of two distributions: dysregulated vs. unaltered mRNAs) followed by a competitive-based [16] enrichment test. Self-contained (non-competitive) methods use exclusively the gene expression values of a gene set, while competitive methods utilize the entire transcriptome as a background [16]. MixEnrich is designed to cluster all mRNAs expression into two groups, unaltered and dysregulated (including up-and down-regulated), using mixture modeling [20]. Then pathways enriched with bidirectionally dysregulated mRNAs are identified using Fisher's exact test [21]. Notably, this method builds on the work of Piccolo and his colleagues who have successfully applied mixture modeling in single samples for a different problem: to identify expressed vs. non-expressed mRNAs [22]. To test the performance of N-of-1-pathways MixEnrich in comparison to the only other single-subject pairedsample gene set tests (Wilcoxon and MD), we performed a simulation study and validation case study. We show that MixEnrich outperforms Wilcoxon and MD under various scenarios of simulated dysregulated pathways. This synthetic result was validated in a case study using head and neck squamous cell carcinomas (HNSCCs) RNA-Seq dataset, where MixEnrich uncovered biological relevant dysregulated pathways.

Datasets
Transcriptome datasets (Table 1) An RNA-Seq dataset of 55 normal lung tissue samples from The Cancer Genome Atlas (TCGA) [23] was used to estimate expression means for each mRNA in the simulation study. To validate N-of-1-pathways MixEnrich, we used another RNA-Seq dataset derived from paired samples of head and neck squamous cell carcinomas (HNSCCs) patients [24].

Knowledge-base dataset
In the HNSCCs case study, gene sets were defined using Gene Ontology Biological Process, GO-BP [13,25]. The GO-BP dataset was retrieved in June 2015 using the org.Hs.eg.db package from Bioconductor [26]. Note, the two terms 'GO-BP' and 'pathway' are interchangeably used in this present study.

An Overview of the methodology of N-of-1-pathways MixEnrich
We propose a novel method, MixEnrich, under the framework of N-of-1-pathways. MixEnrich identifies dysregulated pathways by: (1) clustering mRNAs as unaltered and dysregulated mRNAs and (2) detecting gene sets enriched with dysregulated genes. We named this two-stage procedure MixEnrich for Mixture model clustering followed by an Enrichment analysis. As illustrated in Fig. 1, from a pair of transcriptome derived from a single subject, we constructed a mixture model by modeling the absolute value of log2 transformed fold changes for all mRNAs as a mixture of two distributions: a distribution of dysregulated mRNAs and a distribution of non-dysregulated (unaltered) mRNA expression. We then performed a Fisher's Exact Test (FET) to determine the over-representation of dysregulated mRNAs (both directions) in each pathway [21].

Clustering using the mixture model
For each mRNA, we calculated its absolute value of log 2 transformed fold change, |log 2 FC|, as |log 2 (E 2 /E 1 )|, where E 1 is the expression level of this mRNA in condition 1 (e.g., normal tissue) and E 2 is the expression level in condition 2 (e.g., tumor tissue). Under the mixture model, each mRNA is assumed to belong to a cluster k (unaltered mRNA or dysregulated mRNA) with a prior probability π k . The cluster membership of each mRNA is a Bernoulli trial (Eq. 1).
where Z i is a latent variable and G is the total number of mRNAs in the transcriptome. An mRNA for a gene index i is a member of cluster k when Z i is equal to k. We use x i to represent |log 2 FC i |, and in cluster k, the absolute value of log fold change, x i follows a certain distribution whose parameters need to be estimated. For simplicity, we assumed that the distribution of x i in each cluster followed a normal distribution, whose probability density function is denoted as φ (Eq. 2).
Here μ k andσ k are the mean and standard deviation of the normal distribution for the cluster k. The marginal distribution of X can be obtained by the sum of two weighted normal distributions, hence providing the (discrete) mixture model (Eq. 3).
The estimation of the parameters of the mixture model is implemented by maximum likelihood using an Expectation-Maximization (EM) algorithm [27]. The likelihood that each mRNA belongs to one cluster or the other is assessed by the posterior probability using Bayes rule (Eq. 4).
We defined an mRNA as dysregulated when its posterior probability of belonging to the dysregulated cluster is above 0.5, where the dysregulated cluster is defined as the cluster with the larger mean.

Enrichment of the dysregulated mRNAs
After assigning mRNAs to clusters, a Fisher's Exact Test (FET) was applied to detect the gene sets (pathways) enriched with dysregulated mRNAs [21]. Assume one pathway consists of M mRNAs among which d mRNAs are dysregulated; while the entire genome consists of N mRNAs among which D mRNAs are dysregulated (summarized by a contingency table, Table 2). By this construction, pathway dysregulation is determined relative to the dysregulation of the entire transcriptome as the background. Since different pathways may not be independent due to overlapping mRNAs between them, the p-values resulting from FETs were adjusted for multiple hypothesis testing using the approach developed by Benjamini and Yekutieli [28] that accounts for correlated p-values.
Performance evaluation of the three single-subject methods by simulation Generation of the simulated dataset Single-subject paired RNA-Seq data were simulated to evaluate the performance of three N-of-1-pathways methods: MixEnrich, Wilcoxon, and MD. It has been shown, for biological replicates collected from different subjects, that a negative binomial distribution [29,30] models the distribution of RNA-Seq read counts more adequately than a Poisson distribution [31,32], as the negative binomial distribution accounts for the overdispersion (biological variation) of mRNA expression. However, the overdispersion is assumed to be negligible under N-of-1-pathways framework since it analyzes the paired samples from the same tissue of the same individual [33]. Therefore, the Poisson distribution is employed to simulate 'virtual patients' with various scenarios of dysregulation.
By varying six simulation parameters listed in Table 3, we investigated 107,640 different scenarios of pathway dysregulation. Specifically, for each scenario, we simulated 100 'virtual patients'. Each virtual patient has one dysregulated pathway and one unaltered pathway with the same size. The simulation process is as follows: 1) Estimate the expression mean for every mRNA, g, from 55 RNA-Seq normal lung samples downloaded from TCGA (Table 1). 2) Generate a pair of expression values, Y g1 and Y g2 , for each mRNA g in two conditions (normal vs. tumor) using Poisson distribution, Poisson (λ g ).
3) Generate a dysregulated pathway: a) Randomly sample a proportion (bg.dPct) of mRNAs, in the second transcriptome (tumor), without replacement, and then replace their values by their corresponding values in the first transcriptome (normal) multiplied by a fold change (bg.FC). b) Designate the target pathway by randomly sampling mRNAs (the number of sampled mRNAs = p.S) from the transcriptome without replacement. c) Randomly sample mRNAs (the number of sampled mRNAs = p.S × p.dPct) from the target pathway without replacement, and designate the sampled mRNAs as dysregulated.

Comparing the performance of MixEnrich with Wilcoxon and MD
Using the simulated datasets, we compared the proposed method N-of-1-pathways MixEnrich with two other single-subject methods: N-of-1-pathways Wilcoxon [6] and MD [9]. We evaluated the performance of the three methods by the following measurements:

Area under the ROC Curve (AUC)
For each scenario of pathway dysregulation, we calculated an Area Under the receiver operating characteristic Curve (AUC) value as follows: Each scenario corresponds to 100 'virtual patients' , and each 'virtual patient' possesses one dysregulated pathway and one unaltered pathway. At a given p-value threshold, among the 100 dysregulated pathways (p.dPct > 0), those identified as dysregulated are true positives (TP) and those identified as unaltered are false negatives (FN). Similarly, among the 100 unaltered pathways (p.dPct = 0), those identified as unaltered are true negatives (TN) and those identified as dysregulated are false positives (FP). We calculated the true positive rate (TPR, or sensitivity) and false positive rate (FPR, or Type I error rate) by equation 5 and equation 6: Receive Operating Characteristic (ROC) curves were generated by plotting FPR against TPR at various pvalue thresholds. Areas under the ROC curves (AUCs) were computed approximately using Riemann sum in R.
Area above the 95% contour curve (AAC 95% ) We investigated the interaction effect of two simulation parameters, p.S and p.dPct, on method performance. A contour plot was used to present the joint impact on AUCs induced by the two parameters, p.S and p.dPct, while fixing the other four simulation parameters listed in Table 3. Each point on the contour plot corresponds to an AUC value of a particular scenario of pathway dysregulation. Then the Area Above the 95% contour Curve (AAC 95% ) was calculated as an overall measure of method accuracy when the two simulation parameters vary simultaneously. Specifically, using color-coded values, we plotted AUCs corresponding to any combination of the two parameters p.S and p.dPct while fixing the four other parameters, p.Fc, p.upPct, bf.FC, and bg.dPct. The horizontal and vertical axes in the contour plot represent the values of p.S and p.dPct, respectively. AUC values on the contour plot are indicated by color gradient. All points with an AUC value of 95% on the contour plot were connected to construct the 95% curve, demarcating the ACC 95% boundary.

Validation case study of head and neck cell carcinoma patients
We further evaluated the performance of N-of-1-pathways MixEnrich, in the context of head and neck squamous cell carcinomas (HNSCCs) (Datasets), using paired RNA-Seq data (tumor vs. healthy) from 45 HNSCC patients. Since a vetted gold standard for HNSCCs does not exist and would require experimentally testing pathways, we established 'quasi-gold  [28] for adjusting multiple hypothesis testing. The quasi-gold standard was constructed as the set of all pathways with FDR BY < 0.05. Employing the quasi-gold standard, we compared the accuracy of MixEnrich with that of MD, Wilcoxon, GSEA and DESeq + Enrichment. N-of-1-pathways methods, MixEnrich, MD, and Wilcoxon, are singlesubject methods and were conducted on every single patient of the 15 testing patients. 15 area under the ROC curves (AUCs) were calculated for each N-of-1-pathways methods. Since GSEA and DESeq + Enrichment can only perform on a group of patients, they were evaluated on 50 distinct subsets, which contain 3, 6, or 12 patients, of the 15 testing patients. Taking the subset of 3 patients as an example, 15 testing patients can yield 455 distinct combinations of three patients. To mitigate computational burden, we randomly chose 50 distinct combinations from the 455 combinations as a test set. GSEA and DESeq + Enrichment were conducted on every combination of the 50 distinct patient combinations, which yielded 50 AUCs for each method when compared to the quasi-gold standard. The AUCs resulted from each N-of-1-pathways methods and the AUCs resulted from cohort-based methods performed on 3, 6, or 12 patients were plotted by boxplots. With the same strategy, we also evaluated MixEnrich in the context of breast invasive carcinoma (BRCA) (Datasets) using paired RNA-Seq data (tumor vs. healthy) from 112 BRCA patients (Additional file 1: Figure S1).

Simulation study
To evaluate the performance of N-of-1-pathways Mix-Enrich, we produced synthetic datasets corresponding to 107,640 scenarios of pathway dysregulation by varying six simulation parameters (Table 3). We compared N-of-1-pathways MixEnrich with two other single-subject methods, Wilcoxon, and MD, based on (i) the overall performance across all types of dysregulated pathways (Global comparison of the three N-of-1-pathways methods); (2) change in performance as the value of a single simulation parameter varies (MixEnrich is robust against background noise and bidirectional dysregulation), and (3) the change in accuracy as two critical parameters, pathway size (p.S) and percentage of the dysregulated mRNAs in the target pathway (p.dPct), vary simultaneously (MixEnrich outperforms MD and Wilcoxon when studying the joint effect of pathway size and proportion of dysregulated mRNAs).

Global comparison of the three N-of-1-pathways methods
We compared N-of-1-pathways MixEnrich with MD and Wilcoxon for their overall performance across all types of pathway dysregulation by combing all 107,640 AUCs (Comparing the performance of MixEnrich with Wilcoxon and MD, Fig. 2). Using Wilcoxon signed-rank test [34], the AUCs of MixEnrich are significantly higher than the ones of N-of-1-pathways Wilcoxon (p-value < 1 × 10 −10 ) and MD (p-value < 1 × 10 −10 ). This result is further supported by the boxplots (Fig. 2b) comparing the overall performance across all simulated pathway dysregulation scenarios (107,640 AUCs for each method) suggesting that MixEnrich is preferable to Wilcoxon and MD for single-subject transcriptome analysis to evaluate the dynamic change in gene expression in the presence of background noise or to uncover bidirectionally dysregulated pathways, as detailed in the subsequent sections.

MixEnrich is robust against background noise and bidirectional dysregulation
We further explored the relative effect of each of the six simulation parameters (Table 3)   The increase in percentage of background dysregulated mRNAs (indicative of greater transcriptome noise) does not affect the performance of MixEnrich (Fig. 3, sixth column). Compared to MixEnrich, Wilcoxon and MD are less accurate at various percentages of dysregulated mRNAs in the background (Fig. 3, sixth  column). MixEnrich is a competitive model [16] that compares mRNAs in a pathway against the background transcriptome and distinguishes pathways that are significantly more dysregulated than the background transcriptome. In contrast, MD and Wilcoxon are selfcontained methods [16], implying that they only use mRNA (gene) expression within a given pathway. Therefore, MixEnrich is expected to have a lower false positive rate when there are a lot of dysregulated mRNAs in the background noise such as biological variation and technical artefacts [16,35]. However, as the percentage of background noise increases, the performance declines of Wilcoxon and MD are moderate. The data suggest that bidirectional dysregulation decreases the performance of Wilcoxon and MD more severely than the background noise, and therefore the degenerate effect of background noise is hidden by the effect of bidirectional dysregulation (data not shown). Notably, all three methods perform better as the percentage of dysregulated mRNAs in a pathway (p.dPct), pathway size (p.S), or fold change of the dysregulated mRNAs in a pathway (p.FC) increase (Fig. 3, first, second, and third columns).

MixEnrich outperforms MD and Wilcoxon when studying the joint effect of pathway size and proportion of dysregulated mRNAs
The number of mRNAs in the pathway (p.S) and the proportion of these mRNAs that are dysregulated (p.dPct) are two factors most relevant to biology. A comparison ( Fig. 4 Panel b) of the AAC 95% (Comparing the performance of MixEnrich with Wilcoxon and MD) distributions for the three single-subject methods demonstrates that MixEnrich produced an overall better performance when two parameters, p.S and p.dPct, chang simultaneously. Using Wilcoxon signed-rank test to compare AAC 95% , MixEnrich outperformed both Nof-1-pathways MD and N-of-1-pathways Wilcoxon (p <1 × 10 −10 and p <1 × 10 −10 , respectively). MixEnrich obtained an AAC 95% > 0.8 for 228 of the 234 tested scenarios while N-of-1-pathways MD and N-of-1-pathways Wilcoxon yielded AAC 95% > 0.8 for 15 and 22 of scenarios, respectively. In the scenarios in which MixEnrich yielded AAC 95% < 0.8, the fold change of dysregulated mRNAs was small (1.3 FC) in both the target pathway (p.FC) and the background noise (bg.FC).

Validation case study: pathways uncovered by MixEnrich agree with the quasi-gold standard
We investigated the biological relevance of the dysregulated pathways uncovered by N-of-1-pathways MixEnrich using a biological dataset of RNA-seq paired samples (healthy and cancer tissues) derived from head and neck squamous cell carcinoma patients, HNSCCs [24], presented in Table 1 cohort-based methods GSEA and DESeq + Enrichment in uncovering dysregulated GO-BP terms for HNSCCs.
Since it is not feasible to biologically test each pathway to determine the truly dysregulated pathways and unaltered pathways, we conducted DESeq + Enrichment on the 30 patients of HNSCCs cohort to produce a quasigold standard (Validation case study of head and neck cell carcinoma patients). DESeq identified 4061 differentially expressed genes from the 30 patients, and a big proportion of these genes were recapitulated by the intermediate step of MixEnrich (Additional file 1: Table  S1). The quasi-gold standard consisted of 251 dysregulated GO-BP terms out of the total 3,485 GO-BP terms. MixEnrich achieved higher AUCs (Validation case study of head and neck cell carcinoma patients) in general on predicting the quasi-gold standard in comparison to MD and Wilcoxon (Fig. 5) as well as when compared to AUCs yielded by cohort-based methods conducted across 3, 6 and 12 patients. The superior performance of MixEnrich over cohort-based methods is likely attributed to two reasons: (i) cohort-based methods are underpowered when the sample size is small, and (ii) MixEnrich detects patient-specific signals in addition to the common signals shared among the three patients.
We then tested the hypothesis that single-subject method MixEnrich can capture the individual signals in addition to the common signals shared by all patients. Interestingly, an outlier (patient ID: A6H7) presents in the MixEnrich results, which carries a lower AUC of 0.707. We investigated the dysregulated pathways identified by MixEnrich from patient A6H7 but are not present in the quasi-gold standard (Additional file 1: Table S2). Most of those pathways are related to cell cycle, DNA damage repair, and inflammatory response. Further, all of the GO-BP terms that are identified as dysregulated by MixEnrich from all 15 testing patients exist in the quasi-gold standard (Additional file 1: Table S3).
We also performed another case study using a dataset of matched healthy and cancer RNA-seq samples derived from 112 breast invasive carcinoma patients (Table 1) and again observed the superior performance of MixEnrich (Additional file 1: Figure S1).

Limitations and future work
As noted in 3.1.3, MixEnrich does not perform well when the FC of dysregulated mRNAs is small in both the background and the target pathway. In addition, the Cohort-based methods were performed across 3, 6 and 12 patients (Pt). The number of distinct subjects is shown below the horizontal axis as human icons to further illustrate how many distinct subjects are required in cohort-based analyses to obtain improvements of the AUC (vertical axis). In addition, the three single-subject analyses predict between 200-300 candidate pathways at FDR = 1%, while cohort-based statistics operating on 3 to 12 individuals predict only 50 pathways at FDR = 5% and over 200 at FDR = 20% (data not shown), which explains in part the observed differences in accuracies use of a Poisson distribution in the simulation study can be considered a limitation of our work, as negative binomial distribution could have been employed to introduce more noise in the simulated paired samples. Another possible limitation of our study is the choice of the Expectation Maximization (EM) algorithm [27] for estimating the parameters of the mixture model. This algorithm is not guaranteed to converge towards the global optimum. Since MixEnrich operates on the log 2 transformed mRNA expression fold changes, it may have higher tendency to discover lowly expressed genes as dysregulated, although making inference on gene sets mitigates the bias towards lowly expressed genes (Additional file 1: Table S7). The datasets used in both the simulation study and the validation case studies contain a large amount of lowly expressed genes (Additional file 1: Table S4); the genes annotated to the dysregulated GO-BP terms identified from each of the 15 testing patients have similar distributions compared to the genes annotated all GO-BP terms investigated in the HNSCCs case study (Additional file 1: Table S5-S6). The twostage process of clustering and enrichment can be viewed as a general framework for paired single-subject analysis. We speculate that more elaborate statistical models could improve the performance of the clustering. Future studies could employ a more general gamma distribution kernel and explore techniques that automatically determine the number of clusters.
Importantly, the simulation study results highlight that MixEnrich detects pathways more dysregulated than the background. This is not addressed by the self-contained methods of MD and Wilcoxon. The fact that selfcontained tests do not perform well according to the criteria imposed in the simulation does not invalidate the use of these approaches. Further, self-contained tests are useful to test a small panel of genes such as obtained by real-time PCR.
The current study is based on two paired samples from a single subject. Further improvements and new features of the N-of-1-pathways analytic tools can provide more statistical power as more comprehensive Nof-1 experimental studies and assays may be conceived. For example, future studies may include (i) multiple biological and technical replicates of both tumor and control samples from a single subject or (ii) multiple omics measurements beyond the transcriptome (e.g., proteome, methylome, etc.). Future improvements will need to address N-of-1 studies designed with time-series datasets using multi-gene measurement and genomic information based on data derived from normal, treated, and withdrawn treatment samples from a single patient. As single-cell transcriptome datasets from a single patient are increasingly being studied [36], N-of-1-pathways framework can be applied and further improved as demonstrated by our recent work in profiling circulating tumor cells using N-of-1-pathways MD [10]. As we strive for precision medicine, we must tackle these challenges to accurately provide personal transcriptome analysis at point of care for diagnosis and prognosis.

Conclusion
Compared to our previously developed N-of-1-pathways methods, Wilcoxon and MD, N-of-1-pathways MixEnrich is more effective in detecting non-concordant pathway dysregulation, better reflecting what one would expect in biological pathways. Moreover, this novel twostage competitive gene set testing strategy provides more resistance to background noise, which is ubiquitous in biological systems. Results based on the head and neck squamous cell carcinomas study demonstrate that the dysregulated pathways discovered using MixEnrich overlap highly with the quasi-gold standard compared to the two single-subject methods (Wilcoxon and MD). In addition, we have shown the robust performance of Nof-1-pathways MixEnrich operating on single subjects in identifying dysregulated pathways when compared to small-sample, cohort-based methods (DESeq + Enrichment and GSEA).
In this era of precision medicine, it becomes crucial to develop unbiased and personalized transcriptome analytics for single-subject diagnosis and prognosis, rather than using methods that aggregate signals across heterogeneous patients. N-of-1-pathways MixEnrich is an innovative framework that bridges this gap by analyzing paired samples, one patient at a time, and is ostensibly extensible to other quantitative 'omics measurements (e.g., methylome and proteome). MixEnrich is a valuable tool for studying rare and orphan diseases for which sample sizes remain small whereas cohort-based methods are underpowered in that setting. Lastly, the mRNA-and pathway-level analysis performed patientby-patient by N-of-1-pathways MixEnrich offers more interpretable results for biologists and physicians such as dysregulated mRNAs of interest that can be potentially validated and identified as biomarker candidates for diagnosis.

Additional file
Additional file 1: Figure S1. MixEnrich shows higher performance than other single-subject methods. We repeated the case study using another dataset that contains matched tumor and normal samples for 142 breast invasive carcinoma patients. Each boxplot corresponding to the N-of-1pathways methods (MixEnrich in purple, MD in green, and Wilcoxon in orange) consists of 15 AUCs resulting from 15 testing patients. Table S1. Overlap of dysregulated genes (DEG) between the ones in quasi-gold standard and the ones discovered from single patients. Table S2. GO-BP terms do not exist in the quasi-gold standard but are identified as dysregulated by MixEnrich from patient A6H7. Table S3. GO-BP terms identified as dysregulated by MixEnrich from all 15 head and neck squamous cell carcinoma patients (HNSCCs) patients. Table S4. Summary statistics of expression levels for the three data sets (Datasets). According to the first quartile of the three datasets, all three contain a large amount of lowly expressed genes. Table S5. Summary statistics of the expression levels for the genes annotated to the dysregulated GO-BPs identified from each of the 15 testing patients. Table S6. Summary statistics of the expression levels for the genes annotated to the any GO-BPs investigated in the HNSC validation case study. Table S7. True positive rate (TPR) and false positive rate (FPR) of MixEnrich at the final enrichment step (pathways) are improved as compared to the initial clustering step (mRNAs). We randomly chose 1000 pathway dysregulation scenarios from the simulation study (Generation of the simulated dataset). For each dysregulation scenario, we ran MixEnrich and computed TPR and FPR at the enrichment step as described in Comparing the performance of MixEnrich with Wilcoxon and MD. In addition, we computed the TPR and FPR at the clustering step, at which MixEnrich defines a positive case (dysregulated mRNA) as one whose posterior probability of being dysregulated (Eq. 4) is greater than its posterior probability of being unaltered. True positive mRNAs (TP) at the clustering step are the ones that were dysregulated in simulation and also identified as dysregulated by MixEnrich. False positive mRNAs (FP) at this step are the ones that were not dysregulated in simulation but identified as dysregulated by MixEnrich. (DOCX 36 kb)