The presence of lymph node invasion is one of the strongest indicators for prognoses of distant metastasis and survival in most cancers [1, 2]. In the multi-step process of cancer metastasis development, invasion into a vascular or a lymphatic system has generally been believed to be a key step of tumor cell dissemination [3–5]. Once tumor cells acquire abilities of intravasation and survival in an unfavorable vascular environment, they circulate around the whole body parts to form new tumors at the secondary site [6]. While the exact mechanisms of cancer metastasis through blood vessels and lymph nodes are still being studied, it is necessary to explain the processes in a genetic level as a key factor of cancer patients’ prognosis.

Many researchers have devoted their efforts to understand lymph node invasion in breast cancers, because regional lymph nodes are frequently observed as the first site of metastasis [7]. Survival analyses with clinical features showed that lymph node status is generally marked as a top significant factor among conventional clinical features [8–10]. Studies of finding molecular markers using genome-wide expression profiles identified various genetic signatures for prediction of lymph node and distant metastasis [11–19]. However, the associations between conventional clinical features including tumor size, lymph node involvement and distant metastasis (TNM staging [20]) and prognosis are not yet identified in a genetic level. Moreover, the existence of a common gene set for lymph node metastasis in a transcriptional level is unclear [21].

So far, t-test based differential expression analysis or clustering methods between lymph node negative and positive samples have been used to detect corresponding gene sets [17–19, 21]. Although these methods are straightforward and intuitive, there are several inherent problems in them. First, direct comparison within two classes (lymph node negatives and positives) may simplify the subtle changes over cancer progression. Usually, four-stage pathological N (N0~N3) is used to indicate the degree of lymph node invasion in breast cancer; N0 denotes no lymph node invasion observed. Regarding the expression values as longitudinal data to find patterns over lymph node progression might benefit from utilizing known biomedical information. For instance, a gene whose expression is significantly high only at a certain stage (e.g. N2) is hardly accepted as a closely related gene from current metastasis model. However, a two-class comparison (e.g. N0 vs. others) would mark it differently expressed. Second, the effect of a factor (e.g. lymph node invasion) should be separated from the effect of the others (e.g. tumor size or histological subtype). These factors are generally not independent and will lead to false findings unless carefully analysed. And, the validation of inferred gene signatures should be performed on sufficient number of independent sets in a strict statistical manner. The data statistics and characteristics including inherent biases should be recognized, appropriately treated and be properly analyzed in a meta-analytic way.

There are several statistical models applicable for multivariate correlation scoring (instead of two-class based scoring). Linear/Non-linear (multiple) regression and analysis of variance models (two-way ANOVA and MANOVA) have been widely used in various fields. Both models (linear and non-linear), however, have a few weaknesses; a gene expression pattern over lymph node progression is not necessarily linear, and the data has too few time points to be assessed in a non-linear way. ANOVA models are usually used to test if there is a significant difference among the mean values, so it is not robust to inconsistent fluctuations of expression values. In time series analyses, autoregressive moving average model and its variants (ARMA, ARMAX and ARIMA) are widely used especially in electronic engineering and system identification fields, and some unit root tests (for stationarity test in time series including Augmented Dickey Fuller test [22]) have been used in statistics and econometrics as well. However, there are a few difficulties in adapting these models to our problem; the number of time points is very few, intervals are not regular and the stage is a pseudo-time. After reviewing the conventional models, we developed a new multivariate correlation measure specially designed for non-linear and small data point analysis. Nevertheless, the conventional models were applied as well and tested to compare with our measure and two-class based analyses.

Our method, monotonically expressed gene analysis (MEGA), scores gene expression patterns with their non-linear monotonicity over a stage progression of interest. It accumulates all the normalized expressional differences between two consecutive stages (see Methods). If the direction of expressional change is consistently positive or negative, the score increases; otherwise, the sum of differences will be cancelled out. Because there are two non-independent factors (stage T and N), one variable should be fixed while the other variable is being used. In MEGA, a two dimensional matrix is constructed, each dimension of which is composed of four points (N0~N3 and T1~T4, T0 is excluded due to the lack of data) generating totally 16 data points per a gene. So, applying the scoring function to each row or column represents calculating the cumulative expressional changes over one factor while the other is fixed. MEGA also has a weight parameter to emphasize a specific stage transition (e.g. N1→ N2) to capture genes activated or repressed in a particular time range. After calculating scores, top k genes are collected and named N-wise monotonically expressed genes (N-MEG) or T-wise monotonically expressed genes (T-MEG) depending on which factor is used for the analysis. Validation of inferred gene sets can be done in a retrospective way to see how accurately the gene sets classify prognostic outcomes in other independent data. P-values from each test data are integrated by meta-analysis to report more confident accuracy of the gene sets. This is basically one of the most unbiased ways for evaluating usefulness of inferred gene signatures.

If the gene sets show consistence and confident accuracy, a series follow-up analysis can be used for reasoning biological meaning (e.g. common pathways or transcription factors). First, gene set analysis can discover some biological pathways involving in metastasis progression. Considering pathways instead of individual genes as an acting unit of biological phenomena explains how different gene sets are sometimes associated with same conditions. And we can find more succinct way to describe the whole processes. Second, the fact that the genes show similar expression patterns as the cancer metastasis progresses leads us to a hypothesis that some common transcription factors play a crucial role in the process. Here, all the genes are not necessarily causative; rather, they are effect from changes of a fewer number of genes in upper hierarchy. In this case, finding frequently represented motifs from the promoter regions of the gene sets might be a good analysis for discovering the transcription factors. This would be more powerful information in practical applications such as pharmaceutical research and patient treatment.

Totally four gene sets of size 20 are constructed from 278 breast tumor gene expression data (expO database) by applying monotonically expressed gene analysis (MEGA). They are N-wise monotonically expressed genes (N-MEG) and T-wise monotonically expressed genes (T-MEG), which are further divided into positive and negative correlation sets. Given these four gene sets (N-MEG+, N-MEG-, T-MEG+, and T-MEG-), we tested on65 independent breast cancer prognosis data sets downloaded from ONCOMINE database (See Methods for details) how much the expression values of the genes are correlated with prognostic outcomes.

Lymph node-wise monotonically expressed genes (N-MEG)

The result of meta-analysis test with N-MEG+ and N-MEG- is shown in Figure 1. Two gene sets are divided by the vertical separator. Three major analysis types (PRG, STG, and GRD) and seven minor analysis types are denoted in the first and second columns. Each row corresponds to an experiment and each column corresponds to a gene. So a value in a cell is a p-value of a gene in the corresponding experiment. Cells are colored blue when the genes are significantly up-regulated at the study, yellow when down-regulated, black when not significantly regulated, and grey when the genes could not be found in the corresponding experiments; here, up-regulation means genes are up-regulated in bad-prognoses, higher stages, and higher grades.

It is easily shown that N-MEG+ genes are positively correlated with worse prognoses, higher tumor stages and higher tumor grades. Similarly, N-MEG- is negatively correlated. From the p-value matrix, we can calculate integrated p-values using three meta-analysis methods over ten test classes (Figure 2). It is easily found that the N-MEG is highly significant in all types of prognosis analyses (p-values less than 10^-14 in any methods). Except the test for stage M (current status of metastasis), all p-values were less than 0.01. The study of stage M is designed for elucidating differences of gene expression profile between primary tumors and metastatic tumors. The conceptual difference from the prognosis study of metastasis is that while the former describes the status of ‘metastasis occurred’, the latter describes ‘metastasis will occur’. The results of stage N and stage T were intermediately significant; five of the seven studies in the stage N are two class comparisons (N0 vs. others).Correlations with the tumor grade studies were extremely significant. It is also shown that the Stouffer’s Z method gives relatively more conservative results. As the Stouffer’s Z method has been proven to be more robust to a few extreme values [23] and correctable here (see Methods), we will use the corrected version of this method for rest of the study.

In a comparison with gene sets from previous work, the N-MEG showed the highest association with cancer prognoses (Figure 3).Gene sets from a multiple regression and a two-way ANOVA model followed it and other two gene sets (Suzuki et al and Ellsworth et al) showed relatively lower significance. This result implies that the pattern based methods (MEGA, two-way ANOVA and multiple regression models) are more effective than two class direct comparison methods (t-test and clustering) in finding prognosis associated genes. On the other hand, ANOVA and Suzuki set showed the best score with the N stage. Like we already mentioned, most of the existing N stage test sets are based on a t-test within two classes, which is the same method as what Suzuki et al used. In other analysis types including M and T stage grouping and tumor grade, we could not find significant differences among five methods. Abba set was also tested even though the gene set was already had a selection step using prognosis data (selecting 46 top ranked genes from 300 genes, see Methods). The test showed that our gene set was comparable to it (better in metastasis, relapse and overall prognosis) in spite of a significant degree of unfairness.

Overall aspects of gene expression progression along the N stage give explicit explanations of differences among the candidate gene sets (Figure 4). In the N-MEG and multiple regression model-based gene sets show consistent increase or decrease along the N stage independent to the T stage (Figure 4A and 4B). However, gene sets from two class direct comparison methods (t-test and Mann-Whitney test) show certain degree of inconsistency and discrepancy between lymph node phenotypes and gene expression patterns (Figure 4D, 4E, and 4F). This result shows that those gene sets (lymph node positive vs. negative) may contain false positives from abstracting detailed pattern information, and also implicates the reason why N-MEG showed relatively high significance in the prognosis test.

Classification and survival analysis

To show the classification power of the N-MEGs, we conducted a test for 5-year metastasis free survival data from Wang et al [15]. Because the meaning of N-MEG+ and N-MEG- is so clear, we scored the sum of row-normalized z-scores of corresponding genes; adding for 20 N-MEG+ genes and subtracting for 20 N-MEG- genes. For the 286 primary breast samples (91 metastasis in 5-years), the mean score was nearly zero (6.5x10^-13) and the standard deviation was 7.8. From the 51 patients whose scores were bigger than the mean plus one standard deviation, 37 had metastasis in 5-years giving 0.35 of sensitivity and 0.92 of specificity. The overall accuracy was 0.71. An ROC curve was drawn to compare the N-MEGs with other gene sets (Figure 5A). The N-MEGs showed the best classification power. Interestingly, while the two statistical approaches using stage progression (multiple regression and two-way ANOVA) managed to prove a certain degree of usefulness, the studies using two class comparisons did not. Although the result may be further improved by other fancy classifiers with optimization procedures, we can tentatively conclude that observing the signatures of stage progression gives better results. A set of area under curve (AUC) were denoted in Figure 5B. The AUC of MEGA was 0.69 (0.626 ~ 0.757 in 95% confidence limits). To conduct a survival analysis we divided the all 286 patients into three groups of equal size (n=95 for good and poor group, 96 for intermediate group). It is shown that the three groups have distinct metastasis free survival and hazard rate in Kaplan Meier estimation (Figure 5 C and D).

Common TF binding site prediction for N-MEG

Because the N-MEG+ and N-MEG- are already selected from their expression patterns along the lymph node invasion progression, we can hypothesize that they are co-regulated by several core transcription factors. To find candidate common transcriptional regulators, we analyzed upstream regions of the N-MEG and selected significantly over-represented motif binding sites using Pscan program [35]. We selected top 20, 30, and 50 N-MEG and obtained matching mRNA RefSeq sequences from DAVID database [36, 37]. By running Pscan on [-450, +50] upstream regions onto the JASPAR database [38], we found that ELK4 and ELK1 binding sites were significantly over-represented (p-values 7.8×10^-7 and 8.1×10^-6 respectively, data shown in Figure 6). ELK4 and ELK1 (E26 Like Transcription Factor) are previously known as members of ternary complex factor (TCF) subfamily, which forms a ternary complex by binding to the serum response factor and the serum response element in a promoter region of the c-fos proto-oncogene [39] (SAP1 is a previous name for ELK4). This finding supports the results of the pathway analyses in the previous section implicating that ELK4 based serum response mechanism might be a driving force for breast cancer lymph node invasion and metastasis.

We also conducted the same analysis with N-MEG^0→1 and N-MEG^2→3 (β=10). While there was no significant common binding sites in N2→N3 progression, we found STAT1, IRF2, and IRF1 can be common binding transcription factors in 50 N-MEG^0→1 (p-values of 4.4×10^-13,1.1×10^-12, and7.4×10^-12 respectively, data shown in Figure 6). We found that IRF1, who plays a tumor suppressing role, is negatively regulated by competitive transcriptional binding of IRF2, both of which were significantly correlated with tumor stage (p-value 0.001), depth of tumor infiltration (p-value 0.006) and lymph node metastasis (p-value 0.015) in human esophageal cancers [40]. Here, we also suggest that IRF2-IRF1 pathway is likely to be involved in lymph node invasion and metastasis progression in human breast cancers with well-known activities of STAT1 [41–44].

In this study, we proposed a monotonically expressed gene analysis (MEGA) for extracting genes that are related to lymph node invasion and tumor size in breast cancer. We analysed expression patterns over a two dimensional N×T space and provided results of meta-analysis to evaluate the gene sets. The test has been conducted on completely independent data sets. We showed that gene sets selected from the suggested LE’ and TE’ functions are strongly correlated with cancer prognoses including metastasis, relapse and survival, and showed significantly better results than conventional approaches. These functions are specially designed to capture expressional differences between two consecutive stages and consistency of expression patterns as well. The MEGA model also enabled us to analyze the impact of each clinical factor independently, and to inspect a specific stage transition in a cancer progression.

Before concluding our report, it is necessary to reconsider the meaning of linking clinical factors and cancer prognosis in a molecular level. A general relationship among clinical and genetic factors is described in Figure 7. Since a primary tumor first occurred, accumulated genetic and epigenetic aberrations drive the tumor’s progression. The relationship between genetic aberrations and gene expression changes or protein function aberration is strongly established. So it is obvious that progressed tumors have different gene expressions or abnormal proteins. In our MEGA analysis, we focused on the gene expression part. If we found a candidate gene set which connects gene expression changes and lymph node invasion (A in Figure 7), it should be able to explain the relationship to cancer prognosis (C in Figure 7), because the correlation between lymph node invasion and cancer prognosis has been firmly proven (B in Figure 7). We found that previous candidate genes in A rarely found proper explanations of C. Here, the meaning of our study is summarized in two points. First, we tried to improve accuracy in finding candidate genes in A by interpreting gene expression patterns over lymph node progression (MEGA). Second, we provided a credible meta-analysis test procedure to validate the relationship in C. As there still unexplained important factors remain, we have to integrate additional data of other levels to finalize the lymph node invasion related (or causing) genes. But we suggest the future work also should be cross-validated in the different types of relations.

Although the MEGA analysis provided a feasible link to clinical factors and cancer prognoses in a genetic level, some parts remain to be improved. First, the monotonicity can be defined in various ways. Currently we can rarely determine the activities of genes in an absolute expression level. Some genes have higher saturation level so that the gene expression pattern might show a monotonic increase through all the clinical stages. On the other hand, if an activity of a gene is easily saturated, the gene expression pattern above a certain degree would not be informative anymore. Handling and determining the optimized pattern on every gene is almost impossible, but other heuristics will be available using kinetics and text-mining. Second, integration with other omics data always count. It is still a big question; how to connect two different types and levels of information. As shown in Figure 7, protein level information explains another big part of clinical outcomes. In our opinion, those information sets should be integrated to an augmented ‘gene’ entity with other available information like point mutations, SNPs, and CNVs. In this case, the MEGA model has to be revised in its scoring functions. Lastly, finding for driving mechanisms of progression is one of the ultimate goals in this field. We tried to elucidate these mechanisms through pathway analyses and commonly binding transcription factor analyses here, but it is yet to be a striking discovery. After we solve the prior questions, our approach might be more helpful in clarifying the core genes or genetic events that can be essential in therapeutic applications.