Phylogenies Derived from Matched Transcriptome Reveal the Evolution of Cell Populations and Temporal Order of Perturbed Pathways in Breast Cancer Brain Metastases

Abstract

Metastasis is the mechanism by which cancer results in mortality and there are currently no reliable treatment options once it occurs, making the metastatic process a critical target for new diagnostics and therapeutics. Treating metastasis before it appears is challenging, however, in part because metastases may be quite distinct genomically from the primary tumors from which they presumably emerged. Phylogenetic studies of cancer development have suggested that changes in tumor genomics over stages of progression often results from shifts in the abundance of clonal cellular populations, as late stages of progression may derive from or select for clonal populations rare in the primary tumor. The present study develops computational methods to infer clonal heterogeneity and temporal dynamics across progression stages via deconvolution and clonal phylogeny reconstruction of pathway-level expression signatures in order to reconstruct how these processes might influence average changes in genomic signatures over progression. We show, via application to a study of gene expression in a collection of matched breast primary tumor and metastatic samples, that the method can infer coarse-grained substructure and stromal infiltration across the metastatic transition. The results suggest that genomic changes observed in metastasis, such as gain of the ErbB signaling pathway, are likely caused by early events in clonal evolution followed by expansion of minor clonal populations in metastasis (Algorithmic details, parameter settings, and proofs are provided in an Appendix with source code available at https://github.com/CMUSchwartzLab/BrM-Phylo).

Appendices

Appendix

A1 Transcriptome Data Preprocessing

We applied our methods to raw bulk RNA-Sequencing data of 44 matched primary breast and metastatic brain tumors from 22 patients (each patient gives two samples) [31, 38], where six patients are from the Royal College of Surgeons (RCS) and sixteen patients from the University of Pittsburgh (Pitt). These data profiled the expression levels of approximately 60,000 transcripts. We removed the genes that are not expressed in any sample. We also considered only protein-coding genes in the present study. We conducted quantile normalization across samples using the geometric mean to remove possible artifacts from different experiment batches [1]. The top 2.5% and bottom 2.5% of expressions were clipped to further reduce noise. Finally, we transformed the resulting bulk gene expression values into the log space and mapped those for each gene to the interval [0, 1] by a linear transformation.

A2 Mapping to Gene Modules and Cancer Pathways

The protein-coding gene expressions were mapped into both perturbed gene modules and cancer pathways, using the DAVID tool and external knowledge bases [18], as well as the cancer pathways in KEGG database [19]. This step compresses the high dimensional data and provides markers of cancer-related biological processes (Fig. 1 Step 1).

Gene Modules. Functionally similar genes are usually affected by a common set of somatic alterations [30] and therefore are co-expressed in the cells. These genes are believed to belong to the same “gene modules” [8, 37]. Inspired by the idea of gene modules, we fed a subset of 3,000 most informative genes out of the approximately 20,000 genes that have the largest variances into the DAVID tool for functional annotation clustering using several databases [18]. DAVID maps each gene to one or more modules. We did not force the genes to be mapped into disjunct modules because a gene may be involved in several biological functions and therefore more than one gene module. We removed gene modules that were not enriched (fold enrichment < 1.0) and kept the remaining $m_1=109$ modules (and the corresponding annotated functions), where fold enrichment is defined as the EASE score of the current module to the geometric mean of EASE scores in all modules [17]. The gene module values of all the $n=44$ samples were represented as a gene module matrix $\mathbf B _M \in \mathbb {R}^{m_1 \times n}$. The i-th gene module value in j-th sample, $\mathbf B _{M_{i,j}}$, was calculated by taking the sum of expressions of all the genes in the i-th module. Then $\mathbf B _M$ was rescaled row-wise by taking the z-scores across samples to compensate for the effect of variable module sizes.

Cancer Pathways. Although gene module representation is able to capture the variances across samples and reduce the redundancy of raw gene expressions, it has two disadvantages: First of all, lack of interpretability. Specifically, some annotations assigned by DAVID are not directly related to biological functions, and the annotations of different modules may substantially overlap. Secondly, the key perturbed cancer pathways or functions may not be always the ones that vary most across samples. For example, genes in cancer-related KEGG pathways (hsa05200) [19] are not especially enriched in the top 3,000 genes with the largest expression variances. To make better use of prior knowledge on cancer-relevant pathways, we supplemented the generic DAVID pathway sets with a KEGG “cancer pathway” representation of samples $\mathbf B _P \in \mathbb {R}^{m_2 \times n}$, where the number of cancer pathways $m_2 = 24$. The cancer-related pathways in the KEGG database are cleaner and easier to explain, more orthogonal to each other, and contain critical signaling pathways to cancer development. We extracted the 23 cancer-related pathways from the following 3 KEGG pathway sets: Pathways in cancer (hsa05200), Breast cancer (hsa05224), and Glioma (hsa05214). An additional cancer pathway RET pathway was added, since it was found to be recurrently gained in the prior research [38]. See y-axis of Fig. 3d for the complete list of 24 cancer pathways. We considered all the $\sim$20,000 protein-coding genes other than top 3,000 genes. The following mapping of cancer pathways and transformation to z-scores were similar to that we did to map the gene modules.

Until this step, the raw gene expressions of n samples were transformed into the compressed gene module/pathway representation of samples $\mathbf B = \begin{bmatrix} \mathbf B _M^\intercal , \mathbf B _P^\intercal \end{bmatrix}^\intercal \in \mathbb {R}^{m \times n}$, where $m = m_1 + m_2$. The gene module representation $\mathbf B _M$ serves for accurately deconvolving and unmixing the cell communities, while the pathway representation $\mathbf B _P$ serves as markers/probes and for interpretation purposes.

A3 Deconvolution of Bulk Data

1.1 A3.1 Non-convexity of Deconvolution Problem

Theorem 1

The deconvolution problem Eqs. (1–3) below is not convex:

$$\begin{aligned} \min _\mathbf{C , \mathbf F } \quad&f(\mathbf C , \mathbf F ) = \left\| \mathbf B - \mathbf C {} \mathbf F \right\| _\mathrm{Fr}^2,&\end{aligned}$$

(7)

$$\begin{aligned} \text {s.t.} \quad&\mathbf F _{l j} \ge 0,&l = 1, ..., k, \ j = 1,...,n, \end{aligned}$$

(8)

$$\begin{aligned}&\sum \nolimits _{l=1}^k \mathbf F _{l j} = 1,&j = 1,...,n. \end{aligned}$$

(9)

Proof

If the problem is convex, we should have: $\forall \lambda \in (0,1)$, and $\forall \mathbf C _x, \mathbf C _y, \mathbf F _x, \mathbf F _y$ in the feasible domain, the following inequality always holds:

$$\begin{aligned} \lambda f(\mathbf C _x, \mathbf F _x) + (1-\lambda ) f(\mathbf C _y, \mathbf F _y) \ge f(\lambda \mathbf C _x + (1-\lambda ) \mathbf C _y, \lambda \mathbf F _x + (1-\lambda ) \mathbf F _y). \end{aligned}$$

(10)

However, for the following setting:

$$\begin{aligned}&\mathbf B = \begin{bmatrix} -1.38 &{} 0.92 \\ 1.03 &{} -0.15 \end{bmatrix},&\end{aligned}$$

(11)

$$\begin{aligned}&\mathbf C _x = \begin{bmatrix} -1.74 &{} 2.21 \\ 1.00 &{} -3.97 \end{bmatrix},&\mathbf C _y = \begin{bmatrix} 1.03 &{} -0.46 \\ -3.13 &{} 0.16 \end{bmatrix}, \end{aligned}$$

(12)

$$\begin{aligned}&\mathbf F _x = \begin{bmatrix} 0.83 &{} 0.32 \\ 0.17 &{} 0.68 \end{bmatrix},&\mathbf F _y = \begin{bmatrix} 0.09 &{} 0.34 \\ 0.91 &{} 0.66 \end{bmatrix}, \end{aligned}$$

(13)

and $\lambda = 0.5$, we have

$$\begin{aligned} \lambda f(\mathbf C _x, \mathbf F _x) + (1-\lambda ) f(\mathbf C _y, \mathbf F _y) =4.86 < 11.74 = f(\lambda \mathbf C _x + (1-\lambda ) \mathbf C _y, \lambda \mathbf F _x + (1-\lambda ) \mathbf F _y). \end{aligned}$$

(14)

This is contradictory to Eq. (10).    $\square$

1.2 A3.2 Architecture Specifications of NND

In the NND architecture, $\left| \mathbf X \right|$ applies element-wise absolute value, $\text {cwn}\left( \mathbf X \right)$ column-wisely normalizes $\mathbf X$, so that each column of the output sums up to 1. The two operations of Eq. (5) naturally rephrase and remove the two constraints in Eqs. (2–3), and meanwhile fit the framework of neural networks. An alternative to the absolute value operation $\left| \mathbf X \right|$ might be rectified linear unit $\text {ReLU}(\mathbf X )=\max \left( \mathbf 0 , \mathbf X \right)$. However, this activation function is unstable and leads to inferior performance in our case, since $\mathbf X _{lj}$ will be fixed to zero once it becomes negative and will lose the chance to get updated in the following iterations. One may also want to replace the column-wise normalization $\text {cwn}\left( \mathbf X \right)$ with softmax operation $\text {softmax}(\mathbf X )$. However, the nonlinearity introduced by softmax actually changes the original optimization problem Eqs. (1–3) and the fitted $\mathbf F$ is therefore not sparse.

1.3 A3.3 Hyperparameters of NND

We used an Adam optimizer with default momentum parameters and learning rate of $1\times 10^{-5}$ [20]. The mini-batch technique is not required since the data size in our application is small enough not to require it ($\mathbf B \in \mathbb {R}^{m\times n}$, $m=133$, $n=44$). The training was run until convergence, when the relative decrease of training loss is smaller than $\epsilon =1 \times 10^{-10}$ every 20,000 iterations.

1.4 A3.4 Fitting Ability of NND

One might be suspicious whether the neural network fits precisely in practice, since it is based on a simple gradient descent optimization. To validate the fitting ability of NND, we plotted the PCA of original samples $\mathbf B$ and the fitted samples $\hat{\mathbf{B }} = \mathbf C {} \mathbf F$ (Fig. 4). One can easily see that NND provides good model fits to the data.

Fig. 4.

PCA of pathway representation $\mathbf B$ and nnMF fitted $\hat{\mathbf{B }}$. Each dot represents the pathway values of a sample $\mathbf B _{\cdot j}$ or fitted $\hat{\mathbf{B }}_{\cdot j}$. The first two PCA dimensions of original data and fitted data are almost in the same positions, which indicates that NND is able to fit precisely in our application. The number of components is set to be $k=5$ here.

Fig. 5.

Distribution of elements in fraction matrix $\mathbf F ^{\star }$. Since each column of $\mathbf F$ is forced to sum up to be one, a Laplacian prior is applied to the elements of matrix $\mathbf F$. This leads to the sparsity of $\mathbf F ^{\star }$: 24 out of its 220 elements ($k\times n = 5\times 44$) are zeros (threshold set to $2.5 \times 10^{-2}$).

1.5 A3.5 Sparsity of NND Results

See Fig. 5 for distribution of fraction matrix in NND deconvolution results.

1.6 A3.6 Cross-Validation of NND

In each fold of the CV, we used $\hat{\mathbf{B }}=\mathbf C {} \mathbf F$ to only fit some randomly selected elements of $\mathbf B$, and then the test error was calculated using the other elements of $\mathbf B$. This was implemented by introducing two additional mask matrices $\mathbf M _{\text {train}}, \mathbf M _{\text {test}} \in \{0, 1\}^{m \times n}$, which are in the same shape of $\mathbf B$, and $\mathbf M _{\text {train}}+\mathbf M _{\text {test}}=1^{m \times n}$. During the training time, with the same constraints in Eq. (5), the optimization goal is:

$$\begin{aligned} \min _\mathbf{C , \mathbf F _{\text {par}}} \left\| \mathbf M _{\text {train}} \odot \left( \mathbf B - \mathbf C \mathbf F \right) \right\| _\mathrm{Fr}^2 \end{aligned}$$

(15)

where $\mathbf X \odot \mathbf Y$ is the Hadamard (element-wise) product. At the time of evaluation, given optimized $\mathbf C ^{\star }$, $\mathbf F _{\text {par}}^{\star }$, and therefore optimized $\mathbf F ^{\star }=\text {cwn}\left( \left| \mathbf F _{\text {par}}^{\star } \right| \right)$ for the optimization problem during training, the test error was calculated on the test set: $\left\| \mathbf M _{\text {test}} \odot \left( \mathbf B - \mathbf C ^{\star } \mathbf F ^{\star } \right) \right\| _\mathrm{Fr}^2$. We used 20-fold cross-validation on the NND, so in each fold 95% positions of $\mathbf M _{\text {train}}$ and 5% positions of $\mathbf M _{\text {test}}$ were 1s.

A4 Derivation of Quadractic Programming, $\mathbf P (\mathcal {W})$, and $\mathbf q (\mathcal {W}, \mathbf c )$

Recall Sect. 2.5, for the phylogeny $\mathcal {G} = (\mathcal {V}, \mathcal {E})$, the Steiner nodes are indexed as $\mathcal {V}_S = \{1,2,...,k-2\}$ ($\left| \mathcal {V}_S \right| =k-2$), the extant nodes are indexed as $\mathcal {V}_C = \{k-1,k,...,2k-2\}$ ($\left| \mathcal {V}_C \right| =k$). The i-th pathway values of Steiner nodes are denoted as $\mathbf x = \left[ x_1, x_2,...,x_{k-2} \right] ^\intercal \in \mathbb {R}^{k-2}$, and values of extant nodes as $\mathbf y = \left[ y_{k-1}, y_{k}, ..., y_{2k-2} \right] ^\intercal \in \mathbb {R}^{k}$. Since we consider each pathway dimension separately here, the subscript i for $\mathbf x$ and $\mathbf y$ is omitted for brevity. The weight of edge $(u, v) \in \mathcal {E}$ connecting nodes u and v is $w_{u v}$ ($1 \le u < v \le 2k-2$). Denote $\mathcal {W} = \{ w_{u v} \ |\ (u,v) \in \mathcal {E} \}$. The inference of the i-th element in the pathway vector of the Steiner nodes can be formulated as minimizing the elastic potential energy $U(\mathbf x , \mathbf y ; \mathcal {W})$ shown below:

$$\begin{aligned} \min _{ \mathbf x } \quad U(\mathbf x , \mathbf y ; \; \mathcal {W}) = \sum _{\begin{array}{c} (u, v) \in \mathcal {E} \\ v \le k-2 \end{array}} \frac{1}{2} w_{u v} ( x_u - x_v )^2 + \sum _{\begin{array}{c} (u, v) \in \mathcal {E} \\ v \ge k-1 \end{array}} \frac{1}{2} w_{u v} ( x_u - y_v )^2, \end{aligned}$$

(16)

Theorem 2

Equation (16) can be further rephrased as a quadratic programming problem:

$$\begin{aligned} \min _\mathbf{x } \quad \frac{1}{2} \mathbf x ^\intercal \mathbf P (\mathcal {W}) \mathbf x + \mathbf q (\mathcal {W}, \mathbf y )^\intercal \mathbf x , \end{aligned}$$

(17)

where $\mathbf P (\mathcal {W})$ is a function that takes as input edge weights $\mathcal {W}$ and outputs a matrix $\mathbf P \in \mathbb {R}^{(k-2) \times (k-2)}$, $\mathbf q (\mathcal {W}, \mathbf y )$ is a function that takes as input edge weights $\mathcal {W}$ and vector $\mathbf y$ and outputs a vector $\mathbf q \in \mathbb {R}^{k-2}$.

Proof

Based on Eq. (16), $U(\mathbf x , \mathbf y ; \mathcal {W}) \ge 0$. Each term inside the first summation ($v \le k-2$) can be written as:

$$\begin{aligned} \frac{1}{2} w_{uv} (x_u - x_v)^2 = \frac{1}{2} \mathbf x ^\intercal \mathbf P (w_{uv})\mathbf x , \end{aligned}$$

(18)

where

(19)

Each term ($v \ge k-1$) inside the second summation can be rephrased as:

$$\begin{aligned} \frac{1}{2} w_{uv} (x_u - y_v)^2 = \frac{1}{2} \mathbf x ^\intercal \mathbf P (w_{uv}) \mathbf x + \mathbf q (w_{uv}, y_v)^\intercal \mathbf x + C(w_{uv},y_v), \end{aligned}$$

(20)

where

(21)

and $C(w_{uv},y_v) = \frac{1}{2}w_{uv}y_v^2$ is independent of $\mathbf x$. Therefore the optimization in Eq. (16) can be calculated and written as below:

$$\begin{aligned} \min _\mathbf{x } \quad&\sum _{\begin{array}{c} (u, v) \in \mathcal {E}\\ v \le k-2 \end{array}} \frac{1}{2} \mathbf x ^\intercal \mathbf P (w_{uv})\mathbf x + \sum _{\begin{array}{c} (u, v) \in \mathcal {E}\\ v \ge k-1 \end{array}} \left( \frac{1}{2} \mathbf x ^\intercal \mathbf P (w_{uv}) \mathbf x + \mathbf q (w_{uv}, y_v)^\intercal \mathbf x \right) , \end{aligned}$$

(22)

$$\begin{aligned} \Leftrightarrow \min _\mathbf{x } \quad&\frac{1}{2} \mathbf x ^\intercal \left( \sum _{ \begin{array}{c} (u, v) \in \mathcal {E}\\ v \le k-2 \end{array} } \mathbf P (w_{uv}) +\sum _{ \begin{array}{c} (u, v) \in \mathcal {E}\\ v \ge k-1 \end{array} } \mathbf P (w_{uv}) \right) \mathbf x + \sum _{ \begin{array}{c} (u, v) \in \mathcal {E}\\ v \ge k-1 \end{array} } \mathbf q (w_{uv}, y_v)^\intercal \mathbf x , \end{aligned}$$

(23)

$$\begin{aligned} \Leftrightarrow \min _\mathbf{x } \quad&\frac{1}{2} \mathbf x ^\intercal \mathbf P (\mathcal {W}) \mathbf x + \mathbf q (\mathcal {W}, \mathbf y )^\intercal \mathbf x . \end{aligned}$$

(24)

   $\square$

Remark 1

The optimal $\mathbf x ^\star$ of the Eq. (16), or the solution to the quadratic programming problem Eq. (17) can be solved by setting the gradient to be $\mathbf 0$:

$$\begin{aligned} \mathbf P (\mathcal {W}) \mathbf x ^\star + \mathbf q (\mathcal {W}, \mathbf y ) = \mathbf 0 . \end{aligned}$$

(25)

Therefore,

$$\begin{aligned} \mathbf x ^\star = - \mathbf P (\mathcal {W})^{-1} \mathbf q (\mathcal {W}, \mathbf y ). \end{aligned}$$

(26)

Remark 2

Based on the proof, we can derive how to calculate the matrix $\mathbf P (\mathcal {W})$ and vector $\mathbf q (\mathcal {W}, \mathbf y )$.

Initialize the matrix and vector with zeros:

$$\begin{aligned} \mathbf P \leftarrow 0^{(k-2) \times (k-2)}, \quad \mathbf q \leftarrow 0^{k-2}. \end{aligned}$$

(27)

For each edge $(u, v) \in \mathcal {E}$ with weight $w_{u v}$, there are two possibilities of nodes u and v: First, if both of them are Steiner nodes ($u \le k-2$, $v \le k-2$), we update $\mathbf P$ and keep $\mathbf q$ the same:

$$\begin{aligned} \mathbf P _{u u} \leftarrow \mathbf P _{u u} + w_{u v}, \; \mathbf P _{v v} \leftarrow \mathbf P _{v v} + w_{u v}, \; \mathbf P _{u v} \leftarrow \mathbf P _{u v} - w_{u v}, \; \mathbf P _{v u} \leftarrow \mathbf P _{v u} - w_{u v}. \end{aligned}$$

(28)

Second, if u is Steiner node and v is an extant node ($u \le k-2$, $v \ge k-1$), we update both $\mathbf P$ and $\mathbf q$:

$$\begin{aligned} \mathbf P _{u u} \leftarrow \mathbf P _{u u} + w_{u v}, \quad \mathbf q _u \leftarrow \mathbf q _u - y_v \cdot w_{u v}. \end{aligned}$$

(29)

We apply the same procedure to all dimension of pathways $i = 1, 2, ..., m$ to get the full pathway values for each Steiner node.

A5 Differentially Expressed Cancer Pathways

Table 2 provides a list of the identified differentially expressed cancer pathways.

Table 2. Differentially expressed cancer pathways between primary and metastatic samples (FDR < 0.05).

A6 Portions of Cell Communities in BrM Patients

Figure 6 shows the inferred cell community portions across the BrM samples. The figure displays, for each patient, the proportion of each community in the primary and the metastatic sample.

A7 Perturbed Cancer Pathways Along Phylogenies

Tables 3, 4, 5 and 6 provide a full list of perturbed pathways across the phylogenies for Case 1, 2, 3, and 4 in Fig. 3e.

Fig. 6.

Classification of BrM patients based on the consisted cell subcommunities in matched samples. There are 12 subcases of the 4 cases mentioned in Sect. 3.2. Specifically, there are 9 specific cases (Case 1a-i) in Case 1. Most patients (7) have all the five cell communities in both primary and metastatic samples (Case 1i). A few patients (4) have all communities in metastasis samples and all clones but community C3|P in primary samples. The element $\mathbf F _{lj}$ is taken as 0 when it is smaller than a threshold $2.5 \times 10^{-2}$, and therefore the l-th community is missing in the j-th sample.

Table 3. Perturbed pathways during the evolution of cell communities in primary and metastatic tumors (Fig. 3e Case 1). The top five perturbed pathways whose gain or loss greater than 1.0 along each edge of phylogeny are shown. Clinically actionable perturbed cancer pathways during metastasis are shown in boldface, i.e., ErbB, RET, and PI3K-Akt [4, 31, 38].

Table 4. Perturbed pathways during the evolution of cell communities in primary and metastatic tumors (Fig. 3e Case 2). The top five perturbed pathways whose gain or loss greater than 1.0 along each edge of phylogeny are shown.

Table 5. Perturbed pathways during the evolution of cell communities in primary and metastatic tumors (Fig. 3e Case 3). The top five perturbed pathways whose gain or loss greater than 1.0 along each edge of phylogeny are shown.

Table 6. Perturbed pathways during the evolution of cell communities in primary and metastatic tumors (Fig. 3e Case 4). The top five perturbed pathways whose gain or loss greater than 1.0 along each edge of phylogeny are shown.