9. Signals of Death—Post-Diagnostic Single Gene Expression Trajectories in Breast Cancer—A Proof of Concept

Recruitment for the prospective Norwegian Women and Cancer (NOWAC) study started in 1991 (Lund et al., 2008). Women were randomly sampled from the Norwegian population register in Statistics Norway. The women were mailed a letter of invitation and a questionnaire. Follow-up was based on linkage to the Cancer Registry of Norway and the register of deaths were based on the unique national birth number given to all Norwegian inhabitants. Repeat questionnaires were mailed with intervals of some years. In the years 2002–2006, women were invited to participate in a subcohort, the NOWAC Postgenome cohort study; for further details see Dumeaux et al., 2008. The main purpose was to establish a biobank suitable for analyses of functional genomics, in particular transcriptomics. Random samples of NOWAC women were drawn in weekly batches of 500 women until 50 000 women had responded positively. Women were invited to fill in another questionnaire and donate a blood sample at a health-care institution such as a GP’s office. The blood samples were sent overnight to the institute by special post for biological samples. The tube contained a protective buffer that preserved the mRNA in the blood (PAX gene blood RNA system), allowing frozen storage over time and optimizing sensitivity of the analysis.

The present analysis used a subsample of the NOWAC postgenome biobank participants. Women who had both filled in a questionnaire in 1996–1998 and had given a Postgenome blood sample were eligible, a total of 31 101 women. Since collection of blood was at random without knowledge of disease status, the procedure gave a uniform distribution of gene expression measurement over time.

In 2010, breast cancer cases diagnosed between 1996 and 2006 were identified through a linkage to the national cancer registry. An age-matched control was drawn at random from the same batch of 500 women. A total of 394 incident breast cancer cases were identified. Those rendered non-eligible were six technical outliers, seven cases with unknown metastases, seven cases with another incidence of breast cancer before blood collection, ten controls diagnosed with cancer before blood donation, and one who emigrated, leaving 363 case-control pairs for the present analyses.

A linkage to the register of vital status in Norway gave a complete follow-up after blood donation until the end of the study on 31.12.2014, or death or emigration. Causes of death according to different strata of metastatic/invasive cancer at time of diagnoses are given below.

In order to update changes in clinical stage or a second breast cancer and to remove controls with an incidence of cancer, another linkage was performed in 2018 with the Cancer Registry of Norway. For six women with metastases and ten cases with another incidence of breast cancer, the updated information was used to change the start of follow-up.

Of the 363 case-control pairs, 85 were omitted since the follow-up time for the cases that are observed before 18 months from diagnosis are heavily influenced by the treatment. We therefore first analyze a dataset of 39 cases who died from cancers and compare them with 239 cases who did not die of cancer, i.e. a total dataset of 278 case-controls. Later, we reduce this to a dataset consisting of 23 cases with metastatic breast cancer who died of breast cancer and 79 cases with metastatic breast cancer who did not die of cancer; see Table 9.1.

We first focus on identifying genes with a different time development. Let X _{c , g} be the difference in log2 gene expression for case-control pair c = 1,2, ..., M for gene g = 1,2, ... N _g . Further, let t _c be the time of observation relative to diagnosis for the case in the case-control pair c. We assume X _{c , g} ~N(f _{s ( c ) ,g} (t _c ), σ _g ) where σ _g is the standard deviation and s(c) is the stratum of case c. We estimate the function f _{s , g} (t) by taking an average of the observations X _{c , g} from stratum s(c) in an interval that includes the n nearest observations in time, i.e. the n/2 observations with largest t _c but t _c  < t and the n/2 observations with smallest t _c but t _c  > t. The number n is a tradeoff between precision and resolution. It should be large enough that the estimate in an interval should not depend on a single data point and at least smaller then M/4 in order to get resolution in time. If there is a large difference in the time development between the two strata, the test statistic or area V _g  = |f _{a , g}  – f _{b , g} | = ∫|f _{a , g} (t) – f _{b , g} (t)|dt describing the area between the curves, will be large where the two strata are denoted a and b, respectively. This estimate is the sum of the absolute value of the differences in average gene expression between the two strata in equally spaced time points assuming the controls have similar values. The integral is evaluated in a time interval where there are observations from both strata.

We make N _g , independent hypotheses, i.e. one hypothesis for each gene:

H0A: For gene g, the time development of X _{c , g} is independent of the stratum s(c), i.e. f _{a , g}  = f _{b , g} .

For each gene g, we compare the observed V _{g , o} with the variable V from a simulated distribution where we use a standardization of the same variables X _{c , g} for all the genes simultaneously, but where we randomize the strata s(c) of the cases. We maintain the observations for each gene and the number of observations from each stratum. From the N _u simulations, we estimate the probability distribution g(v) = P(V > v) that is independent of the genes. Based on this, we find a p-value p _g  = p(V > V _{g , o} ) = (k + 1)/(N _u N _g  + 1) for each gene g if k of the N _u N _g simulations have V > V _{g , o} . We correct for multiple testing using the (Benjamini & Hochberg, 1995).

We estimated the functions f _{s , g} with a moving average, where the window size is one-quarter of the respective datasets, i.e. 9 and 59 points, respectively. These functions were evaluated in regularly spaced points, making it easy to evaluate the functions when the observations for each stratum changes position in time. The integral was evaluated in the largest interval such that there were data points from the two strata before and after the interval making the interval equal to (547, 2255) days after diagnosis. The method was applied on a dataset with N _g  = 8400 genes. The analysis is performed for standardized gene expressions for each gene

where the standard deviation σ _g is taken over the case-controls pairs for each gene. This normalization is necessary in order to compare area between the curves since we want to focus on the differences in time development and not in the mean value and the variance. The results were based on simulations with N _u  = 1000 realizations.

We will also make a weaker hypothesis where we analyze all the genes simultaneously:

H0B: For all genes, the time development of X _{c , g} is independent of the stratum s(c), i.e. for all genes f _{a , g}  = f _{b , g} .

Note that we only make one hypothesis here. We perform the same N _u simulations as in the hypothesis test for H0A, but we use the test statistics V _{(1), o}  > V _{(2), o}  > ··· which is the V _{g , o} variables for all the genes that are sorted in decreasing order. From the simulation, we find the probability for the ordered variables g _m (v) = P(V _{( m )} > v) for m = 1,2, ..., and the p-value for the hypothesis test p _{( m )} = P(V _{( m )} > V _{( m ), o}  = (k + 1)/(N _u  + 1) if k of the N _u simulations have V _{( m )} > V _{( m ), o} . Here, we have many highly correlated test statistics V _{( m ), o} for m = 1,2, ..., for testing the same hypothesis H0B. The integer m is chosen by the user. m = 1 means that we are only interested in the most extreme gene and m = 10 means that we are interested in the 10 most extreme genes. This method is most interesting for 3 < m < 50, i.e. where no single gene is significant, but several/many genes have deviating values and where we avoid the multiple testing problem. Ideally, m should be decided before the data is analyzed, but this is not as critical as when alternative test statistics are independent of each other.

It is also possible to use the same technique in order to predict the stratum of a case. The idea is to find out if the observations X _{c , g} for g = 1,2, ..., N _g is closest to f _{a , g} (t _c ) or f _{b , g} (t _c ) for the genes with lowest p-values in the hypothesis test H0A above. Our ambition is only to find the quality of the prediction, not to make a diagnosis for each case. Hence, we make the following hypothesis:

H0C: The prediction P _{a , c} , that the case c belongs to stratum a, is independent of the stratum s(c).

We test the hypothesis using cross-validation. The case-control pairs are divided into the D ₁, D ₂, ... D _Nd groups, which are described in more detail further down. For each of the pairs c ϵ D _k we find

where f _{a , g , k} (t _c ) is the estimated gene expression for gene g and stratum a at time t _c , i.e. the time of observation X _{c , g} based on all the data except the data in D _k . This is based on the assumption that X _{c , g} ~N(f _{s ( c ) ,g ,k}(t _c ), σ _g). The weight w _g may be set equal to , possibly modified based on the correlation between the gene expressions for different genes and how significant this gene is for the prediction. How important gene g is for the prediction is estimated from p _{g , k} , the p-value for hypothesis test HOA estimated from all the data except D _k . The prediction that the observation X _{c , g} is from stratum a is then

This model gives probabilities that are approximately uniform in (0,1), see Figure 9.1. If we had assumed X _{c , g} ~N(f _{s ( c ) ,g,k} (t _c ), σ _g ) independently for each gene g, then most P _{c , a} would be close to 0 or 1, which does not correspond to our ignorance in the classification. We use the test statistics

which is the L₁ distance between the prediction for stratum a and the indicator for stratum a. We may randomize P _{c , a} between the observations and find a distribution for S. The p-value for the hypothesis test H0C is found from the distribution for S, i.e. p = P(S < S _o ).

We used cross-validation and therefore needed to divide the dataset into separate groups. The 39 case-control pairs where the case died of cancer and 239 case-control pairs where the case did not die of cancer were divided into 13 groups, D ₁, D ₂, ... D ₁₃. The data in each stratum was divided into three time periods for each of the two strata with an (almost) equal number of observations. Each of the 13 groups had (almost) the same number of observations from each stratum in each of the three time periods. For each group k we find the values p _{g , k} from the hypothesis test H0A based on all the data except the data in D _k based on 1000 randomizations of the strata s(c).

Results from testing the H0A hypothesis are shown in Table 9.2. The function of the top 10 is shown in Chapter 10.

The q-values are the FDR-corrected p-values. From 8400 genes, 733 genes had q-values below the given threshold for hypothesis test H0A (97 with q<0.01 and 636 with q<0.05). The result shows that many genes have a different time development between the two strata. The reduced dataset with only metastatic breast cancer is too small to get significant results on this test. Figure 9.2 shows the functions f _died,g (t) and f _survived,g (t) separately for each of the 12 most significant genes of the 8400 genes (p<0.000001). The test statistics is the area between the pair of curves.

As shown for most genes, the gene expression increases. Noticeably, f _survived,g (t) is almost constant and close to 0 in the entire period while f _died,g (t) is closer to 1 or –1 in the period (1000,1500) days and then for many genes closer to 0 after 1500 days. The normalization (1) implies that the data for each gene have average value 0 and standard deviation 1 in order to compare data between genes. Since the stratum that survived is much larger, it is natural that the average of these curves is smoother and close to 0. The statistical test shows that deviation in averages value between the strata is significant for many genes. The p-value depends on whether there is a systematic difference in level or time variation of the gene expression, not the size of the difference in average value between the strata since this is removed in the standardization.

We also want to test all the genes simultaneously. Since we only make one test, it is easier to reject the hypothesis for a smaller dataset. First, we test hypothesis H0B on the dataset with 39 and 239 case-control pairs. There is only one hypothesis, but we have many different test statistics, one for each of m ordered V_(m) test statistics for the area between the two curves. The different test statistics indicate whether there is a strong difference in the time development in one or a few genes compared to a smaller difference in many genes. The test for each of the ordered variables is highly correlated. Table 9.3 shows the p-values from the H0B.

Notice that we get very significant results and that many genes have a different time development between the two strata.

This test is also performed on the reduced dataset with only metastatic breast cancers. There are only 23 and 79 case-control pairs in the two strata (Table 9.1), those with metastases who died of breast cancer and those who did not die of cancer, respectively. We still get significant p-values, but much larger values than in the larger data set with both metastatic and invasive cancer; see Table 9.4. The differently ordered variables are highly correlated and give typically p-values between 3–6%.

We also want to test whether it is possible to predict the stratum of each case by testing hypothesis H0C. The 13 different datasets leaving out one of the groups D _k give a slightly different rankings of the importance of the different genes. The mean correlation between the rankings of the genes for these 13 datasets is 0.85. Table 9.4 shows that there is a large overlap in the most important genes in the 14 different datasets when we include the ranking using all the data. On average, four of the five genes with the lowest p-value when using the entire dataset were among the 10 smallest p-values in the reduced datasets. We have marked the five genes with the smallest p-values when using all the genes with colors. Notice that many of the same genes have small p-values for the different subsets.

We have tested different predictions methods, i.e. different choices of the weights w _{g , k} . The different choices give highly correlated probabilities. We have found out that w _{g , k}  = 1/p _{g , k} for the 50 genes g with smallest p _{g , k} value for each group is a quite robust choice. Figure 9.3 shows the predicted probabilities for each of the 278 case-control pairs after time of follow-up. Ideally, we wanted all the 39 red and yellow circles to be equal to 1 and the remaining circles equal to 0.

Based on these variables, we find P _{c , a} , S _o and the p-value p = P(S < S _o ) based on the 10 000 randomization of P _{c , a} We find the p-value less than 0.004 indicating that the prediction is far from random. Table 9.5 gives another presentation of the prediction based on whether P _{c , a}  > 0.3 or not.

Increasing the thresholds from 0.3 will decrease both the number of true positive and the number of false positive. The threshold 0.3 is chosen as a balance between false positive and false negative. This gives a sensitivity equal 0.56 and specificity equal 0.69. The mean prediction value for the 39 cases who died is 0.32 and the mean prediction value for the 239 cases who survived is 0.23. The predictions are also shown in the boxplot in Figure 9.4.

Notice that the invasive cases who died and the metastatic cases who survived have a relatively similar prediction which is between the prediction of the metastatic cases who died and the invasive cases who survived.

No unifying theory exists for human carcinogenesis; the number of proposals is many (Vineis, Schatzkin, & Potter, 2010). To date, most mechanistic or pathways analyses have been experimental in-vitro or animal studies. With the increasing knowledge about human carcinogenesis in tumor tissues or in blood at time of diagnosis, some disturbing facts about the validity of the animal models for human carcinogenesis have been brought up. First, the biology of mice and men is comparatively different (Mak, Evaniew, & Ghert, 2014; Anisimov, Ukraintseva, & Yashin, 2005), and a controversial Nature editorial (“Of men, not mice”, 2013) advocated the need for human functional studies. Similarly, the translational value of mouse models in oncology drug development was recently questioned (Gould, Junittila, & de Sauvage, 2015). While cancer can be manufactured in mice quite easily, these models do not necessarily apply to humans (Mak et al., 2014). Consequently, an increasing number of studies use functional genomics as biomarkers, looking both at the exposure relationship and the outcome. While interesting, this approach lacks the distinct focus on the time-dependent process of carcinogenesis. Few, if any, prospective studies have been designed for longitudinal analyses of functional genomics related to the processes of carcinogenesis and metastasis.

The interpretation of these genes points towards important changes in genes known to be affected at breast cancer, and in addition some more general ones.

During the different laboratory steps, several decisions had to be taken on level of noise and the use of specific distribution of noise. Further, since a gene maybe not expressed in all individuals, the percentage of cases or controls with sufficient signals had to be decided. The stronger the criteria moving towards hundred percent, the harder the exclusion.

The strength of the study is the unique biobank created with the purpose of gene expression analysis in peripheral blood. This gave a unique opportunity to study the immune response since the mRNA in blood came from immune cells. This opened for the view that the carcinogenic process not only included exposures to carcinogens, but also has an important counterforce in the immune system. This has been known for more than a hundred years, and today documented through the new immune therapies.

The design has been population-based with a complete follow-up on cancer incidence, emigration and death based on linkage to national registers using the unique national birth number given to all residents in Norway from 1960. In addition, we had access to updated information on metastases and second breast cancers in the time between inclusion and blood donation. This somewhat reduced the noise from carcinogenic processes hidden at the time of diagnosis.