Local in Time Statistics for detecting weak gene expression signals in blood – illustrated for prediction of metastases in breast cancer in the NOWAC Post-genome Cohort

Abstract

Background:

Functional genomics in a processual analysis cover the time-dependent changes in transcriptomics and epigenetics before diagnosis of a disease, reflecting the changes in both life style and disease processes. The aim of this paper is to explore the dynamic, time-dependent mechanisms of the metastatic processes, using blood transcriptomics and including time in a continuous manner. For achieving this goal, we have developed new statistical methods based on statistics that are local in time.

Methods:

The new statistical method, Local in Time Statistics (LITS), is based on calculating statistics in moving windows and randomization. The method has been tested for the analysis of a dataset that collectively provides information on the blood transcriptome up to 8 years before breast cancer diagnosis. The dataset from the Norwegian Women and Cancer (NOWAC) Post-genome Cohort consists of 467 case-control pairs matched on birth year and time of blood sampling. The data for a pair are the difference in log₂ gene expression between the case and control. The stratified analyses are based on important biological differences like metastatic versus non-metastatic cancer, and the mode of cancer detection, ie, screening-detected cancers versus clinically detected cancers. The dataset was used for examining whether the gene expression profile varies between cases and controls, with time, or between cases with and without metastases.

Results:

The null hypotheses of no differences between cases and controls, no time-dependent changes, and no differences between different strata were all rejected. 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 predictor was not very sensitive to the number of genes included.

Conclusion:

Using a new statistical method, LITS, we have demonstrated time-dependent changes of the blood transcriptome up to 8 years before breast cancer diagnosis.

Keywords::

• processual analysis
• transcriptomics
• prediction
• breast cancer
• blood
• Local in Time Statistics

Acknowledgments

We are thankful to and impressed by the women who donated blood for this cancer research project. Bente Augdal, Merete Albertsen, and Knut Hansen were responsible for all infrastructure and administrative issues. We thank Clara-Cecilie Günther for preprocessing the data. This study was supported by a grant from the European Research Council (ERC-AdG 232997 TICE). The funders had no role in the design of the study; in the collection, analyses, and interpretation of the data; in the writing of the manuscript; or in the decision to submit for publication. Some of the data in this article are from the Cancer Registry of Norway. The Cancer Registry of Norway is not responsible for the analysis or interpretation of the data presented. Microarray service was provided by the Genomics Core Facility, Norwegian University of Science and technology, and NMC – a national technology platform supported by the functional genomics program (FUGE) of the Research Council of Norway. The data will be stored at the European Genome-phenome Archive (EGA, https://www.ebi.ac.uk/ega/),³³ where it will be accessible on request.

Author contributions

All authors contributed toward data analysis, drafting and revising the paper and agree to be accountable for all aspects of the work.

Disclosure

The authors report no conflicts of interest in this work.

Supplementary material

Method

Adjusting for the batch effect

Here, we give a short description of the ComBat method developed by Johnson et al¹ for estimating the batch effects and how to use these estimates for adjusting for the batch effects when computing sample means and standard deviations.

The log₂ gene expression value Y_ijg for gene g and sample j from batch i is modeled as

$Y_{ijg}=a_g+X{\beta }_g+{\gamma }_{ig}+{\delta }_{ig}{\epsilon }_{ijg}\text{and}{\epsilon }_{ijk}\sim \text{Normal}\left(0,{\sigma }^2\right),$

where

• α_g is the overall gene expression,

• X is a design matrix for sample conditions,

• β_g is the vector of regression coefficients corresponding to X,

• γ_ig is the additive batch effect, and

• δ_ig is the multiplicative batch effect.

The batch-adjusted data Y_ijg can then be computed as

$Y_{ijg}^{*}=\frac{Y_{ijg}-{\widehat{a}}_g-X{\widehat{\beta }}_g-{\widehat{\gamma }}_{ig}}{{\widehat{\delta }}_{ig}}+{\widehat{a}}_g+X{\widehat{\beta }}_g.$

The estimates of the parameters α_g, β_g, γ_ig, and δ_ig are computed using an empirical Bayes method.¹ Note that in the implementation of the method, the batch-adjusted data Y_ijg are computed as $Y_{ijg}^{*}=\frac{\frac{Y_{ijg}-{\widehat{\alpha }}_g-X{\widehat{\beta }}_g}{\widehat{\sigma }}-{\widehat{\gamma }}_{ig},}{{\widehat{\delta }}_{ig}}\widehat{\sigma }+{\widehat{\alpha }}_g+X{\widehat{\beta }}_g$, where ${{\widehat{\gamma }}^{\prime }}_{ig}=\frac{{\widehat{\gamma }}_{ig}}{\widehat{\sigma }}$ is the parameter that is estimated instead of ${\widehat{\gamma }}_{ig}$.

Both the expectation and the variance of a gene for the cases can vary both with time and stratum. We therefore cannot use the ComBat method for batch adjusting the dataset that consists of differences in log₂ gene expression between cases and controls. Instead we will use ComBat to estimate the batch effects ${\widehat{\gamma }}_{ig}$ and ${\widehat{\delta }}_{ig}$ from a dataset that includes only the log₂ gene expressions for the controls.

Log₂ gene expression data that are adjusted for the additive batch effect γ_ig, but not for the multiplicative batch effect δ_ig, can then be computed as

${Y^{\prime }}_{ijg}=Y_{ijg}-{\widehat{\gamma }}_{ig}={\mu }_{Gg}+{\widehat{\delta }}_{ig}{\epsilon }_{ijg}\text{where}{\epsilon }_{ijg}\sim \text{Normal}\left(0,{\sigma }_G^2\right)\text{for}\text{group G}\text{.}$

For case-control pair c from batch i with sample j₁ as control (from group G1) and sample j₂ as case (from group G2), we have the log₂-expression difference of X_g,c.

$X_{g,c}=Y_{i{j}_{2}g}-Y_{i{j}_{1}g}={Y^{\prime }}_{i{j}_{2}g}-{Y^{\prime }}_{i{j}_{1}g}={\mu }_g+{\widehat{\delta }}_{ig}{\epsilon }_{g,c}\text{where}{\epsilon }_{g,c}\sim \text{Normal}\left(0,{\sigma }^2\right)$

We observe that X_g,c is adjusted for the additive batch effect γ_ig, but not for the multiplicative batch effect δ_ig.

We compute the estimate of μ_g, ${\widehat{\mu }}_g$ as the weighted average of X_g,c, where the weights are $\frac{1}{{\widehat{\delta }}_{ig}}$, and we compute the estimate of ${\sigma }^2,{\widehat{\sigma }}^2,\text{as}\frac{1}{n-1}{\displaystyle \sum _{c=1}^n{\left(\frac{X_{g,c}-{\widehat{\mu }}_g}{{\widehat{\delta }}_{ig}}\right)}^2}$. We will compare estimated sample means and standard deviations between genes. For each gene, we therefore multiply the estimates of ${\widehat{\delta }}_{ig}$ by a constant so that for this gene $\frac{1}{B}{\displaystyle \sum _{i=1}^B{\widehat{\delta }}_{ig}}=1$, where B is the number of batches/runs.

Figure S1 Boxplots illustrating how the score used in the predictor depends on the number of genes included in the score.

Notes: The score has been normalized by dividing with the number of genes included in the score. The score for the cases with metastases should be positive (lower panel), while the scores for the cases without metastases should be negative (upper panel). (A) Scores for case-control pairs around 6 months from the screening group. (B) Scores for case-control pairs around 2 years and 6 months from the clinical group.