6. Statistics of Sparsely Sampled Curves

A central part of the statistical methodology is to examine how gene expression varies with time. By dividing the entire time period into shorter time periods and computing different kind of statistics for each such time period, we simplify the problem to examining how these statistics vary with time. We use overlapping time periods by using a moving window in time. The statistics computed for a time period are independent of time. As the distribution of the gene expressions may vary with time, the lengths of the time periods should be chosen such that we obtain as short time periods as possible. However, to obtain as good estimates as possible for each time period, there should be as many case-control pairs as possible within each time period. There is a trade-off between these two wishes. To conclude, we define M – L + 1 time periods [t ₁, t _L ], [t ₂, t _{L +1}], ..., [t _{M – L +1}, t _M ] where L is chosen such that we obtain short time periods with many case-control pairs. Typically, we let L ≈ M/4.

Compared to an approach where the time periods are not overlapping, an advantage with the moving window approach is that we are better able to identify the points in time relative to diagnosis where changes in gene expression occur.

In most hypothesis tests, we compute p-values by estimating the null distribution for the statistic of the hypothesis test by randomizing the data, i.e. interchanging covariates (time to diagnosis, case/control, etc.) between the patients. In the randomization we preserve critical properties of the genes (level of expression, complex correlation between genes, etc.) and randomize only what is connected to the changes in time, stratum or case/control status. This randomization defines the null distribution for the test statistic that is used when finding the p-value.

We can randomize the data either by randomizing the case and control in each case-control pair, by randomizing the case-control pairs between the periods, or by randomizing between the two strata within the time period. Note that all three randomization strategies maintain the correlation structure between the genes for each case-control pair. Also note that each randomization of the data leads to a different ordering of the genes if the genes are ordered according to the statistic of the hypothesis test.

The p-value of the test is set to , where N _S is the total number of randomizations and K is the number of randomizations out of N _S with a more extreme statistic than the statistic for the real data. When we test one hypothesis for each gene, we take multiple testing into account by using the Benjamin-Hochberg procedure for controlling the false discovery rate, FDR (Reiner et al. 2003).

When extracting information for a time period, we use the gene expression data in that time period and ignore information about time when we compute different kinds of statistics, identify differentially expressed genes etc. We have used the following statistics for a time period:

Let m _p,g be the sample mean and s _p,g be the sample standard deviations for the differences in log₂ gene expressions for gene g in time period p. The variable m _p,g is an approximation to f _{s ( c ) ,g} (t _c ) = E{X _g,c } and s _p,g an approximation to the standard deviation of X _gc in the interval. Let m _p,g,a (m _p,g,b ) be the sample mean and s _p,g,a (s _p,g,b ) be the sample standard deviations for the differences in log₂ gene expression for gene g in time period p for stratum a (b). We define the statistics s _{p ,( g )}, m _{p ,( g )}, and w _{p ,( g )} from these sample means and standard deviations as follows:

In some approaches, we use the sample mean m _{p , g} directly, without ranking.

We identify the significantly differentially expressed genes in the time period using the Bioconductor R-package limma, linear models for microarrays (Ritchie et al. 2015), where the response is X _g,c , i.e. the difference in log2 gene expression between the case and the control of a case-control pair.

In the next section, we describe how these statistics can be used for examining how gene expression varies over time, between strata, and between cases and controls

The objective is to be able to identify small changes that vary slowly in time and/or between strata, by using a large number of genes in each hypothesis test and predictor.

For examining whether there are differences between cases and controls, between strata or in time, we test following hypotheses using the three statistics’ standard deviation, mean and weight defined in the previous section:

H0-case-ctrl: The expectation of X _g,c is zero. This means that there is no difference between the expectations of the log₂ gene expression values for the cases and controls. It implies that f _{s ( c ) ,g} (t _c ) = 0. If the null hypothesis is false, the expectation will be different from zero for some periods and genes. We test the hypothesis by using the statistic m _{p ,( g )}.

H0-time: The distribution of X _g,c is not associated with the time to diagnosis. This means that the expectation and standard deviation of X _g,c are the same in all time periods. It implies that f _{s ( c ) ,g} (t _c ) is constant in time. If the null hypothesis is false, the standard deviation for some periods will be lower than the standard deviations for the entire time period for some genes. Also, the absolute value of the expectation for some periods will be higher than the absolute value of the expectation for the entire time period for some genes. We test the hypothesis first by using the statistic s _{p ,( g )}, and then by using the statistic m _{p ,( g )}.

H0-node: The expectation of X _g,c is not associated with stratum. This means that the expectations for the two strata are equal for all genes g and time to diagnosis t. It implies that f _{s ( c ) ,g} (t _c ) is the same for all strata. If the null hypothesis is false, the difference in expectation will be different from zero for some periods and genes. We test the hypothesis by using the statistic w _{p ,( g )}.

The null distribution of each statistic will be estimated by randomizing the data, and we compute p-values by comparing the statistic for the data to the estimated null distribution.

We will reject the H0-case-ctrl hypothesis if the hypothesis m _{p ,( g )}>0 is rejected for at least one time period p and rank g, where g belongs to a subset of the N _g ranks. In practice, we have chosen to let the subset of ranks consist of ranks between approximately 1% and 25% of the number of genes, so that the subset contains both relatively low and high ranks. This means that H0-case-ctrl is rejected based on a very large number of hypotheses, that are also highly positively correlated, and we therefore needed to adjust for multiple testing. The approaches for rejecting the H0-time and H0-node hypotheses are similar. Besides rejecting the three null hypotheses, the hypothesis tests for the statistics for each time period and rank can be used for illustrating how the p-values are associated with the time to diagnosis.

In this approach, we examine changes in gene expression over time by examining how the number of differentially expressed genes between cases and controls varies with time. The time curve in this case consists of the number of differentially expressed genes in each time period. Such time curves give an indication of when there is a large difference between cases and controls before or after diagnosis. For comparing different strata, we can compare the time curves of the strata.

When testing whether the number of differentially expressed genes are different for two strata, we use n _s , the number of genes that are differentially expressed between cases and controls in at least one time period, as test statistic for stratum s.

When comparing stratum a and stratum b, we cannot directly compare n _a and n _b if M _a , the number of case-control pairs in stratum a, is much larger than M _b , the number of case-control pairs in stratum b. To use the same number of case-controls pairs when comparing the strata, we test the following hypothesis.

H0-strata1: The number of differentially expressed genes between cases and controls is different for stratum a and b. Assume M _b  < M _a . We want to estimate the null distribution for n _a when the sample size of stratum a is M _b , and then compare this distribution to n _b . The null distribution is found using simulation by repeatedly sample M _b case-control pairs from stratum a, and then compute the number of differentially expressed genes for each sampled dataset. The p-value of the hypothesis test is computed by comparing the samples of the null distribution to n _b .

In this section we also estimate the function f _{s ( c ) ,g} (t _c ) = E{X _g,c } by taking the mean of observations of X _g,c in an interval of time. However, here we analyze the properties of the function in the entire time period at the same time instead of for each interval in time separately. Then we are able to analyze all the data from a gene instead of focusing on properties in an interval. This may be used to identify genes that are significantly different between strata and predict strata from the gene observations.

We estimate the area between two curves f _{a , g} (t) and f _{b , g} (t) by the test statistics V _g  = |f _{a , g}  – f _{b , g} | = ∫|f _{a , g} (t) – f _{b , g} (t)|dt where the two strata are denoted a and b, respectively. This is equal to the weighted sum of the absolute value of the differences in average gene expression between the two strata in each time interval, where the weight depends on the length of the time interval. The area is large if there is a large difference between the curves and we neglect whether this is due to different average value or one is increasing and the other is decreasing in time.

We test the following hypothesis:

H0-strata2: The functions f _{a , g} (t) and f _{b , g} (t) are equal. For each gene g, we compare the observed V _g with the same variable from a simulated distribution where we resample the variables X _g,c for all the genes simultaneously by randomizing the stratum s(c) for the case-control pairs. We maintain the observations for each gene and the number of observations from each stratum. We have made N _g simultaneous tests and need to use the methods for adjusting for multiple testing.

The weights w _{p ,( g )} can be estimated for each rank g from data in period p for a training dataset. The stratum of the case of a new case-control pair, i.e. a case-control pair that does not belong to the training set, can then be predicted based on the score

,

where x _(g) is the difference in log₂ gene expressions of the new case-control pair and δ _{p ,( g )} is 1 if the weight w _{p , g ʹ} is positive, and –1 otherwise, where |w _{p , g ʹ}| = w _{p ,( g )}. The n genes with highest absolute value of the weights are used for computing the score, where n is a number less than or equal to the number of genes, N _g . Large values of z indicate that the new case belongs to stratum a. If z > c, for some arbitrary threshold c we conclude that the new case belongs to stratum a, otherwise we conclude that the new case belongs to stratum b. We may set c = 0 if it is not more important to avoid false classification in one stratum relative to the other and if

,

where m _{p ,( g ), a} and m _{p ,( g ), b} are the sample means that are used when computing w _{p ,( g )}. Increasing (decreasing) c results in fewer false positives (negatives) at the cost of more false negatives (positives).

The available datasets are too small to be divided into a training and validation set. When predicting the stratum of the cases in a dataset, we should therefore use a leave-one-out or k-fold cross validation approach. When using the leave-one-out approach, we predict the stratum of case j using weights w _{p ,( g )} that have been estimated using the dataset where case-control pair j has been excluded. The k-fold cross validation approach is similar, except that we divide the dataset into k folds and predict the stratum of the cases in fold f using weights w _{p ,( g )} that have been estimated using the dataset where the case-control pairs in fold f have been excluded.

In a previous methodological study, time was categorized in three non-overlapping periods (Lund et al. 2016); see Chapter 8. The aim in that study was to show that there is signal in the data, but without showing where in time the changes in gene expression occurred or which genes were involved. The main idea used in the paper is that genes can be grouped into curve groups, each curve group corresponding to genes with a similar development over time (Figure 6.1). Based on these curve groups, we tested a set of hypotheses that determined whether there is development in gene expression levels over time, and whether this development varies among different strata. For a breast cancer dataset in the Norwegian Women and Cancer (NOWAC) postgenome cohort, the curve group analysis revealed that development of gene expression levels varied in the last years before breast cancer diagnosis, and that this development differed by lymph node status and participation in the Norwegian Breast Cancer Screening Program. The effect of the participation may be due to different treatment for the participating women representing the majority of the population.

Partly to be able to better identify the points in time relative to cancer diagnosis where changes in gene expression occur, we developed methods based on a moving window in time. This approach includes time in a more continuous manner than the approach based on curve groups. In Holden et al. 2017 we used moving windows and randomization for the same dataset as we used in Lund et al. 2016, Chapter 8. The null hypotheses of no differences between cases and controls, no time-dependent changes, and no differences between different strata were all rejected. The main conclusion of the analyses was that there are time-dependent changes of the blood transcriptome up to eight years before breast cancer diagnosis.

In Chapter 9 we examined the changes in gene expression after diagnosis for a breast cancer dataset in the Norwegian Women and Cancer (NOWAC) postgenome cohort. We stratified stage in invasive or metastatic breast cancer, and vital status in dead or alive at the end of follow-up, and observed a significant increase of differentially expressed genes among women with metastatic disease who later died both compared to invasive cases that survived (p = 0.001) and to metastatic cases that survived (p = 0.024). To illustrate the difference between strata, we made heatmaps for the most differentially expressed genes over time; see Figure 6.2 below.

We also observed a second transient increase in blood gene expression a few years after diagnosis in metastatic cases, hypothetically representing a capitulation of the immune system.

The analysis of complex curves was extended in order to identify genes that are significantly different between two strata The method tests differences in a non-parametric time development relative to time of diagnosis of the gene expressions from different strata using the area between the curves in a long time period.

The method was tested on case-control differences in log2 gene expressions in a post diagnostic time period separating between the women who survived and the women who died of breast cancer. The method clearly showed non-linear changes, with rapid transient mostly increasing fold changes, in cases who later died. Survivors had no changes. For cases that died, this transient increase was followed by a regression towards the gene expression profiles of survivors. For 9786 genes, the integrated area from 18 months to 8 years was highly significant (p < 0.00001) among women who died. There were indications of a stronger relationship in metastatic cases alone.

Figure 6.3 below shows the f _{s ( c ) ,g} (t _c ) curves for the 100 most significant genes in the period after diagnosis for the women that survived and died of breast cancer.

This systems epidemiology approach provides a proof of concept for the use of gene expression as an individualized biomarker of prognosis related to death or not. Since we had no prior knowledge of the shape of differences in gene expressions as a function of time relative to diagnosis, we needed a non-parametric model that identified all possible changes in trajectories. The aim of the study was to explore single gene expression trajectories from immune cells in blood over the first years after diagnosis as predictors of later vital status, dead or alive.

In Holden et al. 2017, we described Prediction of stratum based on local statistics, to illustrate how the predictive power of the test varies with time. In Figure 6.4 below, we observe that for screening detected cancers the probability of correct prediction of metastasis status was best in year 1 before diagnosis compared to year 3 and 4 before diagnosis for clinically detected cancers.

The mortality of breast cancer is strongly associated with parity, i.e. the number of children to whom a woman has given birth (Lund 1990). In Lund 2018, we used one of the breast cancer datasets in the Norwegian Women and Cancer (NOWAC) postgenome cohort mentioned in the previous sections to examine whether there is an association between gene expression and parity. In that study we used only one time period, which means that we ignore information about time and as a consequence will mainly be able to identify genes with expressions that do not vary with time in the years before/after diagnosis. As there is a large body of evidence demonstrating the long-lasting protective effect of each full-term pregnancy (FTP) on the development of breast cancer (BC) later in life, this is reasonable.

We used the Bioconductor R-package Limma, linear models for microarrays (Ritchie et al. 2015), to identify the genes that were influenced by parity. In the linear model, the responses were X _g,c , the differences in the log2 gene expression for each case-control pair, while we included the parity of the control and the parity of the case as covariates. In the analyses, we merged parities 1–3 and 4–6 so that the parity data consisted of three different values: 0, 1–3, and 4–6. The merging was done in order to reduce the effect of the highest parities. We identified gene sets that were influenced by parity using Limma in the same way as we did for individual genes, by using enrichment scores for gene sets instead of differences in the log2 gene expressions as responses in the linear model. The enrichment scores for gene sets were obtained from the X _g,c s using the Bioconductor R-package gene set variation analysis, GSVA (Hänzelmann et al. 2013).

We found that 756 genes showed linear trends in cancer-free controls, false discovery rate (FDR) 5%, but this was not the case for any of the genes in the breast cancer cases. Gene Set Enrichment Analysis, GSEA of immunologic gene sets, C7 collection in Molecular Signatures Database, MSigDB (GSEA MSigDB, Subramanian et al. 2005) revealed 215 significantly enriched human gene sets (FDR 5%). These marked differences in gene expression and enrichment profiles of immunologic gene sets between breast cancer cases and healthy controls suggest an important protective effect of the immune system on breast cancer risk.