FL benchmark

FL is one of the most common B-cell non-Hodgkin's lymphomas (NHLs), the key genetic lesion (found in ~90% of FL samples) is the t(14;18) rearrangement. This translocation causes the constitutive expression of the antiapoptotic BCL2 oncogene (Bende et al, 2007). FL shows a relatively small network dysregulation signature, with only 192 LoC/GoC interactions. BCL2, which supports eight of those interactions, is ranked first by our enrichment analysis method. By comparison, differential expression analysis between FL samples and GC samples (the normal FL counterpart) ranks BCL2 in the fifty-ninth position. Furthermore, the analysis identified the SMAD1 gene, ranked sixth. This gene, although not detectable by differential expression analysis in our data set, has been shown to have an aberrant pathway activation in FL and other NHL phenotypes, mediated by tumor-transforming growth factor-β (Munoz et al, 2004).

BL benchmark

BL is endemic among children in equatorial Africa and occurs sporadically in other geographic areas, where it also affects adults (Bellan et al, 2003). In these malignancies, a key oncogenic lesion is the translocation of the proto-oncogene MYC from chromosome 8 to either the immunoglobulin heavy-chain region on chromosome 14, or one of the light-chain regions on chromosome 2 or chromosome 22. MYC has been shown to have a global regulatory role in BL (Li et al, 2003). MYC is also one of the most connected hubs in the BCI, having 4079 probe-based interactions. Sixty of these interactions were dysregulated, giving this gene the fifteenth most significant enrichment score. By differential expression analysis between BL and GC cells (BL's normal counterpart), MYC has a rank of thirty-two (see Table II). While this result is encouraging per se, our method was also successful in identifying other key effectors of MYC in BL. In particular, MTA1, an established target of MYC, was ranked third, even though it is not even ranked in the top 1000 genes by differential expression. MTA1 was recently identified as a primary downstream effector of MYC function. Specifically, its silencing blocks the ability of MYC to produce a pathologic transformation (Zhang et al, 2005).

MCL benchmark

MCL is an aggressive type of NHL that generally occurs in middle-aged and elderly people. Cyclin D1/BCL1 (CCND1) is a cell-cycle protein that is overexpressed in MCL as a result of the translocation t(11;14) involving the immunoglobulin heavy-chain gene on chromosome 14 and a region on chromosome 11 harboring CCND1. (Miranda et al, 2000). In the BCI, cyclin D1 is connected to six dysregulated interactions, ranking it eighteenth in our list. By differential expression analysis with non-GC samples (MCL's normal counterpart) CCND1 has a rank of six (see Table II). In addition, our analysis ranked HDAC1 third among all candidates. Histone deacetylases inhibitors have recently been suggested as potentially useful in the therapy of MCL (Heider et al, 2006), so this finding is another piece of supporting evidence that our method identifies the correct patterns. HDAC1 is also highly differentially expressed, and ranked fourteenth. These results indicate that in some cases conventional analysis do indeed capture the correct gene(s). However, as shown, our method seems to consistently identify these key genes as well as effectors, which may be undetectable by differential expression.

In these three cases, it is important to note that we expect the translocated gene to be differentially expressed. It is significant therefore, that against a benchmark where differential expression should be very useful, our method still outperforms it in two out of three cases, and consistently ranks these genes at the very top throughout.

Interestingly, when the scores for these phenotypes are shown distinctly for LoC and GoC interactions (see Table II), MYC appears heavily weighted toward GoC, BCL2 toward LoC, and CCND1 shows a mixed mode of both. These results may indicate that the progression of these lymphomas is marked by distinct types of changes in the network.

Biochemical validation

Although the above examples provide some evidence that our method can correctly identify key regulators and effectors in three separate tumors, a more robust form of validation can be provided by a biochemical perturbation of a specific pathway. We proceeded to analyze a set of samples from Ramos (BL) cell lines stimulated with CD40 ligand or antibody against a non-stimulated set. To quantitatively measure the performance of the method, we considered an established signature of 41 genes in the CD40 pathway and used the gene set enrichment analysis (GSEA) (Subramanian et al, 2005) to compare our method to differential expression analysis.

Our method ranked 379 probes as having a non-zero score. Using GSEA, this ranked list produced a nominal enrichment P-value of 0 (P<1e−3 given 1000 permutations), with 13 of the CD40 pathway genes appearing in the list, many of them clustered at the very top. Remarkably, of the top 10 genes five are in the CD40 pathway set, including CD40 itself, which is ranked ninth. The other four CD40 pathway genes include NFKB1 (second), NFKBIA (third), NFKBIE (fifth), and NFKB2 (tenth), all known to be key effectors of CD40 signaling. Since our method produces a score of zero for all genes that do not participate in any dysregulated interactions, it is not possible to analyze enrichment beyond these 379 probes. When compared with differential expression using the same cutoff of 379 probes, GSEA produces a nominal P-value of 0.12, showing no statistically significant enrichment of the CD40 pathway gene list. CD40 itself is ranked twenty-fourth. Furthermore, in our analysis, we find eight CD40 pathway genes in the top 25 (P-value=0 by Fisher's exact test, below machine precision), compared with only 4 of 25 by differential expression analysis (P-value <2e−5). Although both approaches show significant enrichment, the new method captures twice as many relevant genes within the top 25, while finding the actual perturbation target within the top 10. This further supports the use of our method for the identification of targets of compounds of unknown activity. When looking at these results, the extreme enrichment of the CD40 pathway members, both in the top 10 and 25 genes is likely to make the difference between identifying and missing the perturbation MOA. Note that, similar to the other benchmarks, CD40 itself is upregulated upon binding the CD40 ligand. Thus, as expected, differential expression analysis appears partially effective. However, as shown for MTA1, SMAD1, and other effectors (see Figure 2), IDEA does not require the gene to be differentially regulated in order to be identified as a likely candidate.

Open in a separate window

Figure 2

BL module: A network visualization of the top 25 scoring genes in BL. Transcription factors are shown as circles, whereas other proteins are shown as squares. Protein–protein interactions are also shown in beige, protein–DNA interactions are black with an arrowhead, and transcription factor-modulated interactions are shown in blue with a circular endpoint. Red/green indicates overexpression or underexpression (P<1e−8), respectively in BL versus GC cells. There are some notable characteristics of this figure. First, all 25 genes form a connected module, which would not occur by chance. Second, MYC appears to be a central regulator of this module, as a full 21 out of the 25 members are MYC targets. MYC also appears regulated by MYC-associated zinc-finger protein (MAZ), which is also not differentially expressed. Third, there are interesting sets of genes that emerge, such as SMAD1, which is known to be associated with some NHL, and members of the NFAT family, including NFATC3, NFATC4, and NFAT5 (these proteins are members of the Wnt-signaling pathway). There also appears to be a protein complex of COL1A2, COL6A1 COL6A2, and FN1, which are all upregulated (and members of the cell signaling and ECM–receptor interaction pathway). These module diagrams can serve as a useful platform for further hypothesis generation and biochemical investigation.

Visualization and interpretation

One benefit of a network-based approach candidate is that gene lists can be viewed in a network context. When we map the top scoring genes from the phenotypes listed above across the network, they tend to tightly cluster in specific areas. Figure 2 shows a visualization of the top 25 genes predicted in BL, which form a connected module. Of interest is the fact that MYC is a key regulator of this module (with 21 of 25 genes being its target, including MTA1). These ‘cancer module' diagrams provide more context than a ranked list of genes, and as shown, can effectively complement existing methods such as differential expression.

IDEA is useful for generating testable hypotheses in a number of different contexts. In the first case, ranked genes can be viewed in a network module to identify key regulators. As discussed in Figure 2, this approach would identify MYC, which upon visualization clearly controls the vast majority of top ranked genes. These candidate driver genes could be experimentally validated using siRNA knockdowns or other perturbation assays. Second, these lists can be analyzed for enrichment in specific pathways. We compared the ranked output to a set of Kyoto Encyclopedia of Genes and Genomes, or KEGG (Kanehisa et al, 2006), pathway annotations. For BL, this method identified focal adhesion (P=0) and the ECM–receptor interaction pathway (P=0), which contain similar sets of genes, which are more commonly associated with solid tumors. Also identified were the B-cell receptor-signaling pathway (P=0.006) and the Jak-Stat-signaling pathway (P=0.057), which has been associated with several different cancer phenotypes. Lastly, genes that score high across multiple phenotypes could be identified pertaining to common mechanisms. When the scores across all phenotypes are averaged, the top scoring genes contain several key oncogenic regulators. Included in the top of this list are MYC, the tumor repressor PRDM2, JAK3, the transcriptional repressor DRAP1, and the estrogen receptor ESR1. Ranked second was the transcription factor POU6F1, which is known to have a role in several eukaryotic development processes, but has not been previously associated with lymphoma, and may warrant further investigation.

We applied this approach to the analysis of chronic lymphocytic leukemia (CLL), a complex tumor phenotype, for which oncogenic lesions have not been identified. The top-ranked genes include PRDM2, MYC, and MLL, which are known to be translocated in different subtypes of leukemia, and SMAD3, which is active in several NHL phenotypes. The top 25 genes also form a tightly connected cluster, with almost half the connections being modulated interactions. Pathway enrichment identified the cell-cycle (P=0), B-cell receptor (P=0.0007), TGFβ (P=0.038), and P53-signaling pathways (P=0.05). These pathways are commonly associated with B-cell lymphomas and this is not surprising, but the presence of MYC, MLL, and PRDM2, all strong oncogenic effectors, may be worthy of inquiry in CLL, as they have not previously been associated with this malignant phenotype. MYC shows a high level of connectivity in the module diagram, connecting to 18 out of the 24 other genes. It is also predicted to be a regulator and modulator of PRDM2. As translocations of MYC and MLL are exceedingly rare in CLL (Reddy et al, 2006), it is unclear what role they have in this specific cancer.

Network assembly

The BCI is a mixed-interaction network composed of protein–protein (PP) and protein–DNA (PD) interactions in a human B-cell context (Lefebvre et al, 2007). The former include both same-complex protein interactions and transient ones, such as those supporting signaling pathways. This network has since been enhanced (C Lefebvre et al, in preparation) to include additional post-translational interactions predicted by the MINDy algorithm (Wang et al, 2006). These interactions include those cases where the ability of a transcription factor (TF) to regulate its target(s) (T) is modulated by a third protein (M) (e.g., an activating kinase). The BCI is generated using ‘gold-standard' evidences from curated databases, by applying a Naïve Bayes classifier to integrate a large number of experimental and computational evidence. Evidence is drawn from several sources, including literature mining from GeneWays (Rzhetsky et al, 2004), transcription factor-binding motif enrichment, orthologous interactions from model organisms, and reverse engineering algorithms, including ARACNe (Basso et al., 2005; Margolin et al, 2006) and MINDy for regulatory and post-translational interactions, respectively. A likelihood ratio (LR) for each evidence source was generated using the positive and negative gold-standard sets. Individual LRs are then combined into a global LR for each interaction. A threshold corresponding to a posterior probability P0.5 was used to qualify interactions as present or absent. See the Supplementary Information for full details of the method.

Dysregulation analysis

Analysis was performed using a large compendium of over 200 microarray expression profiles in B cells (BCGEP), including primary tissue as well as cell line samples, available in the NIH Gene Expression Omnibus (GSE2350). Samples in this set were hybridized to the Affymetrix HG-U95Av2 GeneChip^®. After filtering for uninformative probes (those having less than a mean of 50 and a coefficient of variation less than 0.3 in the BCGEP), 7907 remained for analysis. Hierarchical clustering was performed to identify relatively homogeneous phenotype groups suitable for this analysis. The three benchmarking phenotypes used included BL (26 purified and unpurified samples), FL (six purified samples), and MCL (eight purified samples). Other phenotypes represented in this data set included germinal center (GC), naïve (N), memory (M), CLL from mutated (CLL-mut) and unmutated (CLL-unmut) subsets, diffuse large B-cell lymphoma (DLCL), and primary effusion lymphoma (PEL). A list of the analyzed cancerous phenotypes can be seen in Table II. For the CD40 perturbation analysis, a set of 24 CD40-stimulated Ramos cell line samples was used against a background of 43 Ramos samples.

For each phenotype, each BCI interaction was analyzed in sequence to determine if it could be classified as either a GoC, LoC, or no change (NC). The test was based on the estimate of the MI between the expression profiles of the two genes in the interaction. MI is an information theoretic measure of statistical dependence, which is zero if and only if two variables are statistically independent. It was calculated using Gaussian kernel estimation (Margolin et al, 2006). Specifically, we tested whether the MI increased (LoC) or decreased (GoC) when the samples corresponding to the specific phenotype were removed from the entire compendium (used to compute the background MI). A null distribution was computed to assess the statistical significance of an MI change as a function of the background MI and of the number of removed samples. More detailed interpretations of LoC and GoC events are shown in Box 1. See full details in the Supplementary Information.

Scoring

Genes were scored by the enrichment of their direct network neighborhood in GoC/LoC interactions, using a Fisher' exact test. Specifically, for both LoC and GoC, two partial P-values were separately computed, based on the number of dysregulated interactions a gene was directly involved in or it was modulating within its direct neighborhood. A global P-value was computed as the product of all four partial P-values. All scoring totals can be seen in the Supplementary Information, where the score is the negative log of the global P-value.

This work was supported by the NCI (R01CA109755), the NIAID (R01AI066116), and the National Centers for Biomedical Computing NIH Roadmap initiative (U54CA121852). KMM is supported by the NLM Informatics Research Training Program (5 T15 LM007079-15).