A molecular staging model for accurately dating the endometrial biopsy
Other
OA: gold
CC0
⤵ 3 in-corpus citations
AI-generated summary
This study developed and validated a molecular staging model based on global gene expression that precisely dates the endometrial cycle, revealing synchronized daily gene changes and enabling new analyses of endometrial data.
One-sentence paraphrase of the abstract; not a substitute for reading it. No clinical advice. How this works
AI-generated deep summary
This study developed and validated a “molecular staging model” to determine human endometrial cycle stage more precisely than methods based on luteinising hormone surges, hormone measurements, ultrasound, last menstrual period, or subjective histopathology. Using RNA-seq data from endometrial biopsies taken at diagnostic laparoscopy for suspected endometriosis (96 samples for model development and 236 samples in the main study, with cycle stages initially classified by routine pathology), the authors fit splines to expression of 20,067 genes to estimate post-ovulatory day and then extend the approach across the whole cycle using seven pathology stages, reassigning proliferative samples into early/mid/late based on gene-expression patterns. The model revealed synchronized daily gene-expression changes across more than 3400 endometrial genes, with the strongest shifts in the secretory phase, and enabled reinterpretation of published endometrial transcriptome datasets; a noted caveat is that the model relies on accurate, agreement-based pathology staging as the reference. Relevance to endometriosis: the study was performed on women with suspected endometriosis undergoing diagnostic laparoscopy and explicitly positions accurate endometrial molecular staging as important for studying endometriosis, while it also reports that nearly all women experience endometrial-related disorders including endometriosis and adenomyosis. This paper is centrally about endometriosis — it includes endometriosis status in the study cohort and frames the staging model as foundational for understanding endometriosis-related endometrial dysfunction.
Read from the paper's body, not the abstract. Not a substitute for reading the paper. No clinical advice. How this works
Abstract
Natural variability in menstrual cycle length, coupled with rapid changes in endometrial gene expression, makes it difficult to accurately define and compare different stages of the endometrial cycle. Here we develop and validate a method for precisely determining endometrial cycle stage based on global gene expression. Our 'molecular staging model' reveals significant and remarkably synchronised daily changes in expression for over 3400 endometrial genes throughout the cycle, with the most dramatic changes occurring during the secretory phase. Our study significantly extends existing data on the endometrial transcriptome, and for the first time enables identification of differentially expressed endometrial genes with increasing age and different ethnicities. It also allows reinterpretation of all endometrial RNA-seq and array data that has been published to date. Our molecular staging model will significantly advance understanding of endometrial-related disorders that affect nearly all women at some stage of their lives, such as heavy menstrual bleeding, endometriosis, adenomyosis, and recurrent implantation failure.
Full text
98,081 characters
· extracted from
oa-pdf
· 13 sections
· click to expand
Methods
currently in use for estimating endometrial cycle stage
have limitations. Endocrine methods measuring the luteinising hor-
mone (LH) surge or peripheral blood estrogen and progesterone are
indirect and do not allow for variability over time in endometrial
response. Ultrasound scans to detect follicle size and/or ovulation do
not provide an obligatory correlation with endometrial development.
Recording the commencement of last menstrual period (LMP) gives an
accuratefix on a major endometrial event, but as a singlefixed point in
the cycle is of limited use for accurately comparing different stages of
cycles of variable length. Histopathology of the endometrium is the
most direct measure of endometrial stage and normalcy
14, although
this is a subjective method with inherent inaccuracy even among
experts
15. Although significant advances have been made using endo-
metrial gene expression to determine cycle stage, particularly in the
mid-luteal phase around the time of embryo implantation 4,16,17,t h e s e
Methods
do not cover the whole cycle.
A more precise method for normalising endometrial gene
expression across the menstrual cycle will provide a major contribu-
tion to understanding endometrial function and provide foundational
information to investigate the pathophysiology of common gynaeco-
logical conditions such as heavy menstrual bleeding, recurrent
implantation failure and endometriosis.
Therefore, thefirst aim of this study was to develop and validate a
new method for accurately determining menstrual cycle stage based
on changing endometrial gene expression. The second aim was to
demonstrate the functional utility of the new method by using nor-
malised endometrial gene expression data to identify genes that
change expression most rapidly across the menstrual cycle, as well as
investigate the effects of increasing age and ancestry on differential
gene expression in the endometrium. We have previously demon-
strated strong genetic effects on endometrial gene expression with
some evidence for genetic regulation of gene expression in a men-
strual cycle stage-speci fic manner
18,19. However, to date no-one has
identified differentially expressed endometrial genes between women
of different ancestries, despite well-established differences in genetic
makeup.
Results
Subject Details
The median age of all subjects (study 1, 2 and Illumina HT-12 validation
study) at time of endometrial biopsy was 33 years (range 18–49). Of the
total of 358 subjects, 214 had con firmed endometriosis, 131 did not
have endometriosis and in 13 endometriosis status was unknown.
Similarly, 167 had had a prior clinical pregnancy, 183 had never been
pregnant, and pregnancy status information was unavailable for the
remaining 8.
The average age at time of endometrial biopsy of the 236 sub-
jects in Study 1 from which the final molecular model was developed
was 31.1 years (range 18 –49). All these women provided endometrial
biopsies at the time of diagnostic laparoscopy for suspected endo-
metriosis, with the primary symptom for investigation being pelvic
pain. 168 (71%) had endometriosis and 68 (29%) did not. Of the
236 subjects, 60 (25%) had had a successful prior pregnancy (Sup-
plementary Table 1). Fertility intention information was available for
136 of the women (Supplementary Table 2). Only 28 women repor-
ted problems conceiving (defined as trying for more than 12 months
to conceive), of whom 10 went on to have successful pregnancies,
and only 6 out of 96 subjects reported pregnancy loss due to mis-
carriage, although this number could have been higher due to
missing data. All subjects reported regular menstrual cycles and
normal endometrium as assessed by at least one experienced
pathologist.
Analysis 1: Development of the ‘molecular staging model’ to
assign cycle stage for secretory stage samples only
Splines were fitted to RNA-seq expression data for each of 20,067
genes from 96 endometrial samples where 2 or 3 independent
pathology reports agreed to within 2 post-ovulatory days (Fig. 1,
panel 1). For each endometrial sample, an estimated post-ovulatory
day (POD) was obtained using the day which minimised mean squared
error (MSE) between the observed expression and the expected
expression across all genes. Examples of MSE plots are shown in the
Fig. 1, panel 2. There was a strong correlation between the POD cycle
time calculated from the lowest MSE value and the average of the
pathology estimates (r = 0.9297) (Fig. 1, panel 3). To illustrate that
larger, less precise, units of time can be used to estimate cycle time
using the same method, an additional model was built using the
pathology-assigned 3 secretory cycle stages (i.e., early-, mid-, and late-
secretory). The cycle time estimated from the 3 stages model showed a
strong correlation to the cycle time estimated from the 14-day POD
model (r = 0.9807) (Fig. 1, panel 4).
Analysis 2: Molecular staging model using 7 pathology stages for
the whole cycle with RNA-seq and array expression data
In Analysis 2 we modelled RNA-seq expression data from all 236 sam-
ples collected in Study 1. These samples had been classified by routine
pathology into 1 of 7 cycle stages. Because the majority of the pro-
liferative phase samples were not assigned as early, mid or late by the
pathologist, we re-assigned all samples labelled as proliferative into
early, mid, and late byfitting a penalised cubic regression spline (k = 3)
using gene expression data from samples classified by the pathologists
as menstrual, proliferative, and early secretory (Fig. 2a). Then a pro-
liferative time point was estimated from the minimised MSE between
the observed expression and the expected expression across all genes
(Fig. 2b). The proliferative samples were then split into equal sized
groups of early, mid, and late using this time point (Fig.2c). A penalised
cyclic cubic regression spline (k = 8) wasfit for all 20,067 genes using
the 7 stages of the menstrual cycle, which included the re-assigned
early, mid and late proliferative samples (Fig. 2d). Each endometrial
sample was then assigned a‘day’ or ‘model time’ using the time which
minimised the MSE between the observed expression data for all genes
and their corresponding gene models (Fig. 2e). ‘Model time’ is a rela-
tive timepoint in the cycle and does not correspond to a real day.
Under the assumption that all 236 women were approximately uni-
formly distributed across the menstrual cycle, the data were trans-
formed so that the distance in time between each sample was identical
Article https://doi.org/10.1038/s41467-023-41979-z
Nature Communications| (2023) 14:6222 2
(Fig. 2f). This ranked all the samples in order from the start to the end
of the cycle, removing the need for cycle stages or an idealised 28-day
cycle. At this point the x-axis was changed to show the percentage of
the way through the menstrual cycle that each sample was. The new
time points were also compared to the pathology-derived cycle stages
to get an approximation how the model time corresponds to stages in
the menstrual cycle (Fig.2g ) .G e n ec u r v e sw e r et h e nr efitted using the
newly derived cycle times for each sample with a penalized cyclic cubic
regression spline (k = 30) (Fig. 2h). For visualisation purposes, nor-
malisation of gene expression for cycle stage was then derived by
subtracting the expected expression from the observed expression
(i.e., calculating the residuals) and re-adding the mean (Fig.2i).
Validation of the molecular staging model
Various validation studies were undertaken using the molecular sta-
ging model. As an initial check, data from Analysis 1 using POD to
develop the secretory model was plotted against secretory stage data
from the final molecular staging model generated in Analysis 2
(Fig.3a). A second comparison confirmed that using only 3 cycle stages
(early, mid and late secretory from only 1 pathologist) gave similar
Results
to having more frequent daily POD information from 2 or 3
independent pathologists (Fig.1, panel 4). To assess the repeatability
of the molecular staging model method, Analyses 1 and 2 were repe-
ated using Illumina HT-12 data and the results compared for the 198
samples that had both RNA-seq and Illumina HT-12 data (Fig.3b). There
was a high level of correlation in cycle stage determination using data
from the 2 different gene expression platforms, with slightly more
variation being seen in the mid-proliferative phase. Additionally, vali-
dation using unsupervised methods with initial groupings based on a
PCA plot (Supplementary Fig 1) and not using pathology dating
information at all also yielded a high level of concordance between the
molecular staging model and the validation model (Supplementary
Fig 2). The correlation between the two models was 0.989 and the
mean absolute difference between estimated sample cycle times was
1.67%. Peripheral blood estradiol and progesterone levels were not
used to help determine cycle stage and could therefore be considered
as an independent variable. Estradiol ( N =1 5 9 ) a n d p r o g e s t e r o n e
(N = 187) data were plotted against molecular staging model cycle
stage and showed typical expected menstrual cycle distributions
(Fig. 3c, d).
Reanalysis of published data
The molecular staging model was used to re-analyse 3 published
endometrial gene expression datasets available on GEO (GSE65099,
endometrial samples from GSE141549, and GSE180485). These
datasets were chosen because they contained endometrium from
natural cycles with attached estimates of cycle stage and plots. We
first produced a principal component analysis (PCA) plot using our
own RNA-seq dataset ( N = 266) with cycle stage as determined by
the molecular staging model (Fig. 4a). This PCA plot has a char-
acteristic pattern with samples clustering according to cycle stage
as determined using the molecular staging model, with no outliers.
We then generated PCA plots from each of the 3 published datasets.
The PCA plot using data from GSE141549
20 is shown in Fig. 4b, with
samples labelled as per information in GEO as menstrual, pro-
liferative, secretory and unknown. There is mixing of proliferative
samples with menstrual and secretory ones within the PCA plot
when using the GEO assigned labels. The same PCA plot using data
from GSE141549 but with cycle stage assigned by our molecular
staging model has minimal overlap between different cycle stages
(Fig. 4c), demonstrating that the molecular staging model accu-
rately aligns with PCA analyses of endometrial gene expression data
Repeat steps 1 and 2
using the 3 secretory
stages instead of POD
Validation 2
●
●
●
●●
●
●
●
●●
●
●●
● ●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
r =0.98075
6
7
51 0
Cycle Time from POD Secretory Model
Cycle Time from 3 Stages Model
Step 2:
Estimate post-ovulation day
of 96 samples
1.0
1.5
2.0
1 2 3 4 5 6 7 8 9 1 01 11 21 31 4
Post−Ovulation Day
Mean Square Error
Estimate
Model
Pathology
Sample X210475
P2
0.5
1.0
1.5
123456789 1 0 1 1 1 2 1 3 1 4
Post−Ovulation Day
Mean Square Error
Estimate
Model
Pathology
Sample X210470
P2P1
1.0
1.5
2.0
123456789 1 0 1 1 1 2 1 3 1 4
Post−Ovulation Day
Mean Square Error
Estimate
Model
Pathology
Sample X210244
P2 P1
Validation
●●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
r = 0.9297
5
10
51 0
Average Post−Ovulation Day from Pathology
Post−Ovulation Day from Model (5−Fold CV)●●
●
●
●
●
●
●
●
●● ●
●
●
●
●
●●
●
●
●●
●
●
●
●
●
●●
●
● ●
●
●
●
●
●
●
●
●
● ●● ●
●
●
●
●
●
●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
Radj
2 = 0.8107
2
4
6
8
10
51 0
Post−Ovulation Day
log2 (CPM)
ENSG00000106688
●
●
●
●
●
● ●
●
●
●● ●●
●
●
●
●
●
●●
●
●
●
●
●● ●
●
●●
●
●
● ● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●●●●●●●
●
●
●
●●
●
●●
●
●
Radj
2 =0.65534
5
6
7
51 0
Post−Ovulation Day
log2 (CPM)
ENSG00000188313
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
● ●
●
●
●
●
● ●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●●●●●
●●●
●
Radj
2 =0.1728
−4
−2
0
2
4
51 0
Post−Ovulation Day
log2 (CPM)
ENSG00000163464
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●●
●
●
●
●
●
●
● ●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●●
● ●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
● ●
●
Radj
2 = 0.0977
−0.5
0.0
0.5
1.0
1.5
51 0
Post−Ovulation Day
log2 (CPM)
ENSG00000260852
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ● ●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●● ●●●●●
Radj
2 = 0.309
2.8
3.2
3.6
4.0
4.4
51 0
Post−Ovulation Day
log2 (CPM)
ENSG00000103047
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●●
●
● ●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
Radj
2 = 0.482
1
2
3
4
5
51 0
Post−Ovulation Day
log2 (CPM)
ENSG00000167984
Step 1:
Model gene expression of 20,067 genes
1. 2. 3.
4.
P1
Fig. 1 | Analysis 1: Development of the molecular staging model to assign cycle
stage for secretory stage samples only. Panel 1. Examples of regression spline
fitting to expression data for individual genes from 96 endometrial samples taken
between post-ovulatory days (POD) 1–14. Splines were fitted to a total of 20,067
genes. The error bands surrounding the regression lines represent the 95% con-
fidence interval. Panel 2. Plots showing post-ovulatory time that gives lowest Mean
Squared Error (MSE) using spline data for all 20,067 probes for 3 different
endometrial samples (solid line). Dotted lines show POD estimates from 2 inde-
pendent evaluations by experienced pathologists. Panel 3. Correlation between
average POD from 2 or 3 independent pathology evaluations and the POD time at
which the lowest MSE occurred. Panel 4. Correlation between the estimated POD
from the POD secretory model and the estimated cycle time from the 3 stages
secretory model. Source data are provided as a Source Data file.
Article https://doi.org/10.1038/s41467-023-41979-z
Nature Communications| (2023) 14:6222 3
across the menstrual cycle. In a similar fashion, but with a smaller
dataset from GSE65099 21, samples reported as LH + 6 to LH + 10 do
not group in a consistent fashion by PCA (Fig. 4d). When the same
samples are assigned cycle stage times by the molecular staging
model, the same PCA analysis shows consistent grouping according
to cycle stage for all samples, with 2 outlying samples on the PCA
plot being reassigned as proliferative and not secretory (Fig. 4e). We
repeated the PCA comparison approach using 36 RNA-seq endo-
metrial data sets deposited in GSE180485 (Supplementary Fig 3).
These samples are from a study called EndoTime to determine
whether the accuracy of timing of luteal phase endometrial biopsies
based on urinary ovulation testing could be improved by measuring
the expression of a small number of genes using a continuous, non-
categorical modelling approach
17. We produced identical PCA plots
with the first showing samples labelled with percentage cycle times
derived from our molecular model, and the second using the
EndoTime model eLH+ estimates. While our percentage data shows
a steady progression through the cycle in concordance with the PCA
plot, the eLH+ days 4-5 data points are widely spread (i.e., these
samples have signi ficantly different molecular pro files despite the
EndoTime model calling them as similar), and several of the eLH+
days 8 –10 samples are grouped together suggesting that this ver-
sion of the EndoTime model cannot reliably discriminate between
these post LH surge days.
Changes in Endometrial Gene Expression with Increasing Age
and Different Ancestries
Using our RNA-seq data ( n = 266, 20,067 genes analysed) with
menstrual cycle staging calculated using the molecular staging
model, a total of 60 endometrial genes showed signi ficant changes
in expression with increasing age. Examples of 2 signi ficant genes
are shown in Fig. 5a. Re-running the age analysis using the original 7
cycle stage pathology data instead of the staging from the mole-
cular staging model reduced the number of age-related signi ficant
genes from 60 to 32, providing evidence that the molecular staging
model provides a superior approach for identifying differentially
expressed genes. To further explore the effects of aging on endo-
metrial gene expression, an additional n = 87 Illumina microarray
endometrial samples from GSE141549 were analysed and combined
with our RNA-seq differential expression results as part of a meta-
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●●●●●
●
●
●
●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●
●
●
●
●●
●●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
● ●
●
●
●
●
●●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●● ●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●
●
Radj
2 = 0.7396
4
5
6
7
8
9
12345
Pathology Cycle Stage
log2 (CPM)
●●
●●
●
●
●
●
●●
●
● ●● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●●
●
●
●●
●
●
● ●
●
●●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●● ● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●●
●
● ●●
●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●●
●
●
●
●
● ●
●
●
● ●
●
●
●
●
●
●
●
●
1
2
3
4
5
12345
Pathology Cycle Stage
Model Estimate of Cycle Stage
●● ●
●
●
●
●
●●
●
●
●●
●
●
●
●
●●
●
●
●
●
●● ●
●
●
●
●
●
●
●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
● ●
●
●●
●
●
●
●●
● ●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●●
1
2
3
4
5
12345
Reclassified Pathology Cycle Stage
Model Estimate of Cycle Stage
●
●
●
●
●
●●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●● ●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●●● ●
●●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●●
●
●
●
●
●
●
●
●
●
●●
●
●
●●
●
●●
●●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
● ●
●●●
●
●
●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●●
●●●
●●
●
●
●
●
●
●●●●●● ●
●
●●
●
●
●●●●●●●●●●●●
●●●●●●●●●●●●●●●●●●●
●●
Radj
2 = 0.80580.0
2.5
5.0
7.5
10.0
1234567
Reclassified Pathology Cycle Stage
log2 (CPM)
●
●
●●
● ●
●
●
●●
●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
● ● ●●●
●
●
●
●
●●
●
●
●
●●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●●
●
● ●
●
●
●
● ●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●●
●
●
● ●
●
●
●
●
●●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
1
2
3
4
5
6
7
1234567
Reclassified Pathology Cycle Stage
Model Estimate of Cycle Stage ●
●
●●
●●
●
●
● ●
●
●
●● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●●●● ●
●
●
●
●
●●
●
●
●
● ●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●●
●
●●
●
●
●
● ●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●●
●
●
● ●
●
●
●
●
●●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
1
2
3
4
5
6
7
02 5 5 0 7 5 1 0 0
Model Cycle Time
Model Estimate of Cycle Stage
● ●●● ●
● ● ●● ●● ●●● ● ●● ●●● ●●● ●●● ●● ● ●●● ●● ●● ●●
●●
● ●● ●● ● ● ●
●●● ●● ● ●●
●●●● ●● ●
● ●●
● ● ● ●● ●●● ●●● ● ●● ● ●●●
●● ●● ●●
●●● ● ●● ●
●●● ●● ●●●
●
●●● ●● ● ●●● ●● ●●● ●
● ●
●● ●●
● ●● ●●
● ● ●●
●
●●
● ● ●● ●●● ●●●
● ●●
●● ●●
●
●●●
●
●● ●●
●
●
●
●
●
● ●● ●
●
●●● ●●
●
●
●
●
●
●
● ●
●
●
● ●
●
●
●●● ●
●
●
●
●● ●
●
●
●●
●
●
●
●
●
●
●
●
●●
●● ●●
●
●
●
●
●
● ●
●
●
●
●
●
●
●●
●●
●
7
6
5
4
3
2
1
0 1 02 03 04 05 06 07 08 09 0 1 0 0
Model Cycle Time
Pathology Cycle Stage
secretory
●
●
●
● ●●
● ●
● ●
●
●
●● ●
●
●
●
●
● ●
●
●
●
●
●
●●
●
●
●● ●● ●
●
●
● ●
●
● ●
●
●
● ●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●● ●
●
●
●
●
● ●
●
●
●
●●
●
●
●●
●
●
●
●
● ●
●
●
●
●
●
●●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●●
●
●
● ●
●●
●
●
●
●●
●
●
●
●●
●
●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●●
●
●
●
●
● ●
●
●
● ●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●●
●
●
●
●●● ●
●
●
● ●
●
●
●●
●
●
●
●
●
●●●●●●
●●●●
● ●●
●●
●●
●
●
●
●
●
● ●
●
●
Radj
2 = 0.8944
5
6
7
8
9
02 5 5 0 7 5 1 0 0
Model Cycle Time
log2(CPM)
●
●
●
●
●●
●
●● ●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
● ●●
●
●
●
● ●● ●
●
●●
●
●
●
●
●
●
● ●
● ●
●
●
●
●
●
● ●
● ●
●
●
●
●● ●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●●
●
●
●
●
●
●
●
●
●
● ●
●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
● ●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●● ●●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●●
●
●●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
● ●
●
●
●
● ●
●
●
●
●
●
●●
●
●
●
●
●
●
●●
●
● ●
●● ●●
●
●
●
●
●
●●
●
●
●
●
●
●
●●
●
●
●
●
●●
●
●
●
●
●
●
●
●
● ●●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
5
6
7
8
9
02 5 5 0 7 5 1 0 0
Model Cycle Time
Cycle Normalised log2(CPM)
a
d e f
b c
g ih
menstrual proliferative
Fig. 2 | Analysis 2: Development of the molecular staging model to assign cycle
stage for the whole menstrual cycle.Figure 2a. Example of a spline curvefitted to
menstrual, combined proliferative and early secretory expression data. The error
bands surrounding the regression line represent the 95% confidence interval.
Figure 2b. Menstrual, proliferative and early secretory samples with estimated
molecular cycle stage calculated from minimum mean squared error data. Figure
2c. Combined proliferative samples reassigned into early, mid and late proliferative
groups containing equal numbers. Figure 2d. Example of a spline curve fitted to
data from all 7 stages of the cycle using reclassified proliferative cycle stage
information and pathology-derived menstrual and secretory staging. The error
bands surrounding the regression line represent the 95% confidence interval.
Figure 2e. Comparison of cycle time within the 7 cycle stages estimated from the
model and pathology-derived cycle stage. Figure 2f. Under the assumption that 236
women underwent surgery at random stages of the menstrual cycle, data from
236 samples were transformed to be uniformly spaced along the x-axis on a
0–100 scale. This transformation allows each sample to be identified as being a
percentage of the way through the menstrual cycle. Figure 2g. Plot showing rela-
tionship between pathology staging and the percentage of cycle from the mole-
cular staging model. Menstrual is 0–8%, proliferative is 8–58% and secretory is
58–100% of the molecular staging model cycle respectively. Figure 2h, i. Expression
data for ENSG00000187231 replotted using derived‘percentage’ cycle times and
then normalised across the menstrual cycle. The error bands surrounding the
regression line represent the 95% confidence interval. Source data are provided as a
Source Data file.
Article https://doi.org/10.1038/s41467-023-41979-z
Nature Communications| (2023) 14:6222 4
analysis. Considering only genes in both datasets, this reduced the
number of genes tested to 12,868, which still included 32 of the
60 signi ficant genes from our original dataset. 65 signi ficant genes
were found in the GSE141549 data when analysed on its own. How-
ever, when the 2 data sets were combined ( n = 353), 206 signi ficant
genes were identi fied across the whole menstrual cycle (Supple-
mentary Table 3). We then split the samples into 3 cycle stages;
menstrual, proliferative and secretory (equivalent to 0 –8%, 8 –58%
and 58 –100% of the molecular staging model cycle respectively)
and analysed each stage of the cycle separately (Fig. 5b). Of note,
nearly all (218/222 [98%]) of the genes showing signi ficant changes
with age were found in samples taken in the secretory phase of the
cycle (Supplementary Data 1). A gene ontology enrichment analysis
was run using the 218 genes from secretory samples that changed
significantly with age (Supplementary Data 2). The top biological
processes enriched with upregulated genes were related to axo-
nemes, cilia and microtubules while the top processes enriched with
downregulated genes were related to blood vessels, endothelial
cells and angiogenesis.
Ancestry of subjects as defined by a previous study
19 was used to
look for differences in endometrial gene expression using pairwise
comparisons of each ancestry group. In our Australian population the
majority of subjects were of European ancestry, however, despite small
numbers in other groups signi ficant differences in gene expression
were identified between the groups (Supplementary Fig 4, Supple-
mentary Table 4 & Supplementary Data 3).
0
25
50
75
100
Cycle Time
●
●
●
●
●
●
●
●
● ●●
●●●●
●
●
●
●
●
●
●●
●
●● ●
●●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●●
●
●
●●
●●
●
●●
●
●
●
●
●
●
●
●● ●
●●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●●
●
● ●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●●
●
●
●●
●●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
● ●
●
●
●●
●
●
●
●
●
●
●
●
●
●●
●
●
●●
● ●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●● ●
●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
−100
0
100
−100 0 100
PC1
PC2
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●●
−60
−30
0
30
PC2
−50 −25 0 25
PC1
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
−50
0
50
PC2
−50 0 50 100
PC1
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●●
−60
−30
0
30
PC2
−50 −25 0 25
PC1
0
25
50
75
100
Cycle Time
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
−50
0
50
PC2
−50 0 50 100
PC1
0
25
50
75
100
Cycle Time
PROLIFERATIVE
MENSTRUAL
SECRETO
R
Y
a b c
d e
Cycle Stage
●
●
●
●
Menstruation
Proliferative
Secretory
Unknown
LH+ Day
●
●
●
●
●
LH+6
LH+7
LH+8
LH+9
LH+10
Fig. 4 | Reanalysis of published data. Figure 4a. PCA plot using 266 endometrial
samples coloured to show the molecular staging model cycle stage. The PCA plot
has a characteristic pattern with the samples in cyclic pattern corresponding to
approximate cycle stage order. Figure 4b. PCA plot using Illumina microarray data
from GSE141549
20, with samples identified as menstrual, proliferative, or secretory
by the authors. Note significant mixing of proliferative samples with menstrual and
secretory ones. Figure 4c. The same PCA plot using data from GSE141549 with cycle
stage assigned by the molecular staging model. Note minimal overlap between
different molecular staging model cycle stages. Figure 4d. PCA plot using RNA-seq
data from GSE65099
21 with samples identified as LH + 6 to LH + 10 as reported in the
study. Figure 4e. The same PCA plot using data from GSE65099 with cycle stage
assigned by the molecular staging model. Note reassignment of 2 outlying samples
on the PCA plot as proliferative and not secretory. Source data are provided as a
Source Data file.
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
r = 0.9790
5
10Cycle Time from POD Secretory Model
50 60 70 80 90 100
Model Cycle Time
●
●
●●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●●●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
r = 0.9744
0
25
50
75
100Cycle Time from Microarray Model
02 5 5 0 7 5 1 0 0
Cycle Time from RNA−Seq Model
●●●
●
●● ●●● ●●●●● ●●●● ●●●●●● ● ●● ●● ● ●●● ●● ● ●● ●●● ●● ●● ●●● ●● ●●●● ● ●● ● ●● ●●●● ●● ●● ● ●● ●● ●●● ●
●●●● ●
●●●●●●
●●
●
●
●
●●
●●
●
●
●
●
●
●●●
●
●
●●
●
●
●
●
● ●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
● ●
●
●
●●● ●
●
● ●
●
●●
● ●●● ●● ● ●●
●
●●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
0
5
10
15
20
25
02 5 5 0 7 5 1 0 0
Model Cycle Time
P4 ng/mL
●●●● ●
●
●●●
●
●●●
●
●
●●● ●
●
●
●
●● ●
●
●
●● ●
●
●●
●
● ●
●● ●
●
● ●
●
●● ●
●
●
●
●
●
●
●
● ● ●●
● ●
●
●●
●
●
●
●
●
●
●
●● ●
●
●
●
● ●
●
●
●
● ●
●
●●●
●●●●
●
●
●
●
●
●
●●
● ●
●●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
● ●
●
●
●●● ●●
●
●
●
●
●
●
●●●
●
●
● ●
●
●●● ●●
●
●
●
● ●●
0
100
200
300
02 5 5 0 7 5 1 0 0
Model Cycle Time
Estradiol pg/mL
a b c
d
Fig. 3 | Validation of the NGS Molecular Staging Model. Figure 3a. For endo-
metrial secretory samples, POD predicted from the secretory molecular staging
model was compared against the full molecular staging model. Note, POD 1 is
approximately 58% of the way through the cycle. Figure 3b. Illumina HT-12 micro-
array data were also available from 198 endometrial samples used to generate the
molecular staging model. A validation study was run comparing NGS vs Illumina
data. The molecular staging model results showed a strong correlation when
comparing the 2 different gene expression platforms. Figure 3c, d. Peripheral blood
progesterone (n = 187) and estrogen (n = 159) data plotted against the molecular
staging model cycle stage showed expected typical menstrual cycle-stage dis-
tribution. Source data are provided as a Source Data file.
Article https://doi.org/10.1038/s41467-023-41979-z
Nature Communications| (2023) 14:6222 5
Differential gene expression across the cycle using the mole-
cular staging model
To investigate changing gene expression across the cycle, all samples
with RNA-seq data from Study 1 (n = 236) were ranked in chronological
order from start to end of the molecular staging model cycle. A‘sliding
window’ approach was then used to compare differential gene
expression (DGE) between samples 1–8 versus samples 9–16, followed
by samples 2 –9v e r s u s1 0–17, then 3 –10 vs 11 –18 and so on for all
samples across the menstrual cycle. Group sizes were arbitrarily set at
8 because this represents 3.4% of the 236 samples or approximately
1 day assuming a mean cycle length of 28 days. Moderated t-tests were
used to identify differentially expressed genes withP < 0.05 following
multiple testing correction, at each window. Using adjusted P values,
488 unique genes significantly changed expression during menstrua-
tion, 44 during the proliferative phase, and 2921 during the secretory
phase. Peak times of rapid change in gene expression occurred during
menstruation (3% of the way through the cycle), late proliferative
(51%), POD3 (66%), POD5 (71%), POD11 (94%) and POD13 (98%) (Fig.6a).
Examples of 12 endometrial genes showing significant and very rapid
changes in expression across different stages of the menstrual cycle
are provided in Fig. 6b.
The original Endometrial Receptivity Analysis (ERA) publication
identified 238 genes that show major changes in expression before,
during and after the time of embryo implantation at POD 3 –7
22.O f
these 238 genes, 207 were recognised in our NGS data, and 70% of
these (145/207) changed expression significantly between cycle times
66 ± 2% and 76 ± 2% (POD 3–7). Supplementary Fig 5 shows the 6 most
significantly down-regulated genes and the 6 most signi ficantly up-
regulated ERA genes that we identified.
Discussion
We have developed and validated a method for accurately determining
endometrial cycle stage based on global gene expression. We did this
by generating mathematically de fined curves fitted to RNA-seq
expression data from 236 endometrial samples for each of 20,067
genes. From these curves we then found the time of the cycle that
minimised a loss function (mean squared error) which gives the bestfit
for all genes simultaneously for any individual endometrial sample. By
placing the 236 samples in chronological order from the start to end of
the menstrual cycle and converting the x-axis to percentage, we cre-
ated a method for both de fining how far through the cycle any given
sample is, as well as being able to normalise gene expression for cycle
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
● ● ●
● ●
●
●
●
●
●
●
● ●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●●● ●●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
● ●
●
● ●
●
●
● ●
●
●
●
●
● ●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●● ●
●
●
●
●
●
●
●
●
● ●
●
●●
●●
●
●
●
●
●
●
●
●
●●
●
●
●
● ●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●●
●●
●
●
●
●
● ●
●
●
●
●
●
●
●● ●
● ●
●
●
●
●
●
●
● ● ●●●●
● ●
●
●
●● ●
● ●●●● ●●●●● ●●●●● ●●
● ●●●●●●●●●●
padj < 0.001
4
5
6
20 30 40 50
Age
Normalised log2(CPM)
ENSG00000180543
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●●
● ●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●●
●
●
●
●●
●
●
●
● ●
●
●
●
●
●
●● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
● ●
●
●
● ●
●
●
●
●
●
●
●
●●
●
●
●●
●
● ●
●
●
●
●
●
●
●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●●
●
●
●
●● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●●
● ●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
● ●●
●● ●
●
●
●
●
● ●
●●●●●●●●
●●●●●●●●●●
●●
●●●●
●●●●●●●●
padj < 0.001
−4
0
4
8
12
20 30 40 50
Age
Normalised log2(CPM)
ENSG00000129988
●
●
●
●●
●● ●
●
●
●
●
●
●
●
●
●
●
●
● ●● ●● ●
padj = 1
●
● ●
●
●
●
●
●●
●
●●●●
●
●
●
●
●
●
● ●●
●
●
●
●
● ●
●
● ●
●●
●
●
●
● ●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●●
●
●
●
●
●
●● ●
●
●
●
●●
●
●
● ●
●
● ●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●●
●
●
●
●
●●
●
●
●●
●●
●
●
●
●
●
●●
●
● ●
●
●
●
● ●● ●●● ●●●●●● ●●●●●●●
● ●●● ●●
●●●
●●
●●
●
padj = 1
●
●
●
●
●
●
●
●
●
●
●
● ●
● ●
●
●
●
●
●
●
●
● ●
●
●
●
●●●
●
●
●
●●
●
●
●
●
●
●
●●
●
●
●
●
●
●● ●
●
●●
● ●
●
●
●●
●
●
●
●
●
●
●●
●
●●
●
●
●
●
● ●
●
●
●
●
● ●
●● ●
●
●
●●●
●
●
●
●
●
● ●
● ●
●
●●●
●
●
●
●
●
●
●
●
●●●
●
●
●
●
●
● ●
●
●
●●
●
●
●
●
●
●●●●
●● ● ●●●● ●
padj = 0.010
Menstrual Proliferative Secretory
20 30 40 50 20 30 40 50 20 30 40 50
3
4
5
Age
Normalised log2(CPM)
ENSG00000147113
● ●
●
●●
●
● ●● ●● ●
●
● ●
●
●
●
● ●
●● ●●●●●●● ●● ●●
padj = 1
●
● ●● ●
●
● ●●
●
●●
●
● ●
●
●
●
● ●
●
●● ●
● ● ●● ●
●
●
● ●
●●
●
●●
●●●
●●
●
● ●●
● ●
● ●
●
●
●
● ●
● ●
●
● ●
●●
●
●
●
●
●
● ●●
●
●
●
●
●
●●
●
● ●
● ●
●● ●
●
●
●
●
● ●
●
●
●
● ●
●
●
●●
●
●●
● ●●● ●●
●● ●
●
●
●
●●●
●
●
●
●
●●●●●●●●●● ●● ●● ●● ●●●●● ●●●● ●●●●●●●●●●●● ●●●● ●●●●
●
●
●●●●●●●● ●●● ●●
padj = 1
●
●
●
●
●
● ●●●
●
●
●
●
●
●
●
●
● ●● ● ●
● ●
●
●
●
●
●●
●
●
●●●
●
●●
●
●
● ●
●
●●
●
●
●
●●
●●
●●
●
●
●
●
●● ●
●
● ● ●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
● ● ● ●
●
●
●
●●
●●
● ●● ● ●●
●
●
● ●
●●●
● ●●
●
●
●
●● ●
●
●
●
●
●● ●● ●
● ●●●●●●●● ●
● ●●●●●●●●●● ●●●●●●● ●●●●●● ● ●●●●●●
●
●●●●●●●●●● ●●●● ●●●●●●●●●
padj = 0.013
Menstrual Proliferative Secretory
20 30 40 50 20 30 40 50 20 30 40 50
−8
−4
0
4
Age
Normalised log2(CPM)
ENSG00000174175
a
b
Fig. 5 | The impact of age on endometrial gene expression. Figure 5a. Two
examples from the 60 endometrial genes that change expression significantly with
increasing age. Statistical analysis used two-sided empirical Bayes moderated
t-tests implemented in limma35 with corrections for multiple comparisons per-
formed using the Benjamini-Hochberg method (ENSG00000180543: t = 6.025,
p =5 . 8 2 e - 9 ,p
adj = 5.84e-5; ENSG00000164778: t =−5.886, p =1 . 2 9 e - 8 ,padj =8 . 2 2 e -
5). Expression data were plotted following normalisation for changing expression
across the menstrual cycle (n = 266 RNA-seq samples, 87 samples from
GSE14154920). Figure 5b. Expression data for 2 genes plotted separately for men-
strual, proliferative, and secretory samples (n = 266). Statistical analysis used two-
sided empirical Bayes moderated t-tests implemented in limma with corrections for
multiple comparisons performed using the Benjamini-Hochberg method (secre-
tory phase ENSG00000147113: t =−5.074, p =1 . 4 3 e - 6 ,padj = 0.010; secretory phase
ENSG00000174175: t =−4.429, p =2 . 0 9 e - 5 ,padj = 0.013). Source data are provided
as a Source Data file.
Article https://doi.org/10.1038/s41467-023-41979-z
Nature Communications| (2023) 14:6222 6
stage so that gene expression differences between any 2 samples can
be compared.
Our work has immediate translational application. Being able to
accurately identify differentially expressed genes from relatively small
numbers of samples will substantially increase the prospect of identi-
fying genes linked to different endometrial pathologies, thus helping
to identify biomarkers for diagnosis and potential targets for therapy.
Similarly, detailed knowledge of normal gene expression pro files
throughout the endometrial cycle will assist in identifying ‘off-target’
effects during clinical trials of new therapeutics. Accurate cycle staging
is also critical for synchronisation of frozen-thaw embryo replacement
during in vitro fertilization treatment. Longer-term, improved under-
standing of endometrial disorders such as heavy menstrual bleeding,
endometriosis and adenomyosis will lead to improved outcomes for
many women.
Our molecular model for staging the endometrium creates many
new research opportunities in benign gynaecology. There is a com-
pelling argument for more research in thisfield given the considerable
financial and quality of life burdens that endometrial-related disorders
have on patients and their families. The total direct and indirect costs
of abnormal uterine bleeding in the UK and USA were reported more
than 10 years ago as greater than £1.2 billion and $37 billion,
respectively
1. More recently in the USA, heavy menstrual bleeding
(HMB) with uterine fibroids had a mean all-cause total cost of $16,762
per patient, while HMB only was $11,135 compared to a control cohort
at $6,691
23. In 2011, in the UK, it was estimated that in England and
Wales approximately 1.5 million women were affected by HMB
annually. HMB is the fourth most common reason for referral to
gynaecology services in the UK and approximately 20% of 1.2 million
referrals to specialist gynaecology services concern women with HMB
in the UK
1.I na d d i t i o nt ot h e financial impact, women with
menstruation-related symptoms have lower scores on several domains
of quality of life such as general health and physical, mental, social and
occupational functioning during their periods
24. A recent review of the
economic burden of endometriosis determined the overall direct
medical cost range to be from US$1459 to US$20,239 per patient per
year, and the indirect costs between US$4572 and US$14,079
2.
Endometrial cycle stage classi fication is currently based on his-
topathological criteriafirst published over 70 years ago14.T h eo n g o i n g
need for endometrial cycle staging is evidenced by the fact that this
1950 paper has received over 3,500 citations25, the most for any study
in obstetrics and gynaecology. Accuracy of pathology dating has well
established limitations. It has been reported that accuracy in dating the
same POD between different pathologists using Noyes criteria was
0
250
500
750
1000
0 1 02 03 04 05 06 07 08 09 0 1 0 0
Model Time
Number of Significant Genes
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●●
●●
●
●●
●
●●
●●●
●
●
●
●
●
●●
●
●
●
●●
●
●
●●
●
●
●
●
●
●●
●
●●
●
●
●●
●●●
●
●
●
●
●
●
●●
●
●●●●●
●
●
●
●
●
●
●●●
●
●
●
●
●●
●
●
●
●
●
●●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●●
●
●●
●●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●●●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●
●
●●●
●●●●
●
●
●●
●
●●
●
●
●
●
●
●
−5
0
5
0 2 04 06 08 0 1 0 0
Model Cycle Time
log2 (CPM)
left−right determination factor 1
LEFTY1
●
●
●●
●
●
●
●
●
●
●
●
●●
●
●●●
●
●●●●●●●●●
●
●●●●●●
●
●●●●●●●●●●●●●●●●●●●●●●
●
●●●
●
●●●●
●
●●●
●
●●●
●
●●●
●
●●●●
●
●●●
●
●●●●●●●●●●
●
●●●
●●
●
●
●●●
●
●●●●
●●
●●●●●●●●●●●●●●●●●●●●●●●
●
●●●●●●●●●●●●
●
●●●●
●●
●●●●●●
●
●
●
●
●
●●
●
●
●●
●
●
●
●
●●●
●
●
●
●
●
●
●
●
●
●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●●
●
●
●
●
−5.0
−2.5
0.0
2.5
02 0 4 0 6 0 8 0 1 0 0
Model Cycle Time
log2 (CPM)
long intergenic non−protein coding RNA 2274
LINC02274
●
●●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●●●●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●
●
●
●
●●●
●●
●●●
●
●
●●
●●●
●●
●●
●●●●
●
●
●●●
●
●
●
●
●●
●
●
●
●
●●
●●
●●
●●
●●
●
●●●●●
●
●
●
●●
●
●
●
●
●
●
●
●●●
●
●
●
●
●●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●●
●
●
●
●
●
●●
●
●
●
●
●
●●
●
●
●●
●
●
●
●
●
●
●
●
●
●●
●●
●
●
●
●
●
●
●●
●●
●
●
●
●●
●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
2.5
5.0
7.5
10.0
02 0 4 0 6 0 8 0 1 0 0
Model Cycle Time
log2 (CPM)
DEPP1 autophagy regulator
DEPP1
●
●
●●
●●
●
●
●
●
●●
●●
●●●●
●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●●●
●
●
●●
●●
●
●
●
●
●
●●
●●
●
●
●
●
●
●●
●●
●●●
●
●
●
●
●
●
●●●●
●
●
●
●●
●
●●
●
●
●●
●
●
●●●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●●
●●●
●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●●●
●●●
●
●●●
●●●
●
●
●●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●
−5.0
−2.5
0.0
2.5
5.0
7.5
02 0 4 0 6 0 8 0 1 0 0
Model Cycle Time
log2 (CPM)
cytochrome P450 family 3 subfamily A member 5
CYP3A5
●●●●●●●●●●●●●●
●
●●
●
●
●●●●●●●●●●●●●●●●●●
●
●
●●●●●●●●
●
●●●●
●
●
●●
●●●●●●
●
●●●●●●●●●●●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●●
●
●
●
●
●
●
●
●●
●●
●
●●
●
●
●
●
●
●
●●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●
●
●
●
●
●
●
●
●●
●
●
●
●●●●
●
●●●
●
●
●
●
●
●●●
●
●●●●●●
●●
●
●
●
●●●●●●●●●●●●●●●●●●●●●●●●
●
●
●
●
−5.0
−2.5
0.0
2.5
5.0
02 0 4 0 6 0 8 0 1 0 0
Model Cycle Time
log2 (CPM)
3−hydroxy−3−methylglutaryl−CoA synthase 2
HMGCS2
●
●
●
●
●●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●●
●
●●
●
●●
●
●
●●
●
●
●●
●●
●
●
●
●
●●
●●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●●●●
●●
●
●
●
●●
●
●●
●
●
●
●
●●
●●
●●
●
●●
●●●●
●
●●
●
●
●
●
●●
●
●●
●
●●
●
●
●●●●●
●
●
●●
●●●
●●●●
●●
●
●
●
●
●●
●●
●●
●
●
●●
●●
●
●
●
●
●●
●●●
●●
●●
●
●
●
●●●●
●●
●
●●
●●
●
●
●
●
●
●
●
●
●●●
●
●
●
●
●●
●
●●
●
●
●
●
●
●
●
−5
0
5
02 0 4 0 6 0 8 0 1 0 0
Model Cycle Time
log2 (CPM)
phospholipase A2 group IVF
PLA2G4F
●
●●
●
●●
●
●●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●●
●
●
●
●●●
●●●
●●●●●●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●
●
●
●●
●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●●●
●
●
●
●
●
●
●
●
●●
●
●
●●
●●
●
●
●
●
●
●●
●
●
●
●
●
●●●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●●
●
●
●
●
●
●●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
−4
−2
0
2
4
02 0 4 0 6 0 8 0 1 0 0
Model Cycle Time
log2 (CPM)
coagulation factor II, thrombin
F2
●●●
●●●
●
●
●
●
●
●
●
●●
●●
●
●
●
●●
●
●●
●
●
●
●
●
●
●
●●●
●
●
●
●
●
●
●
●
●
●
●●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●●
●
●●●
●●
●
●●
●●
●●
●
●
●
●
●●
●
●
●
●
●
●
●
●●
●
●●
●
●
●
●●●●●
●
●●
●
●
●
●
●
●
●
●
●●
●●
●
●
●
●
●
●●●
●
●
●
●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●●
●
●
●●
●
●
●
●
●
●●
●●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●●●●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●●●
●
●
●
●
●
●●
●●
●
●
●
●
●●
●
●
●
●
●
●●●
●
●●
●
●
●
●
●
●
●
●
−2.5
0.0
2.5
5.0
02 0 4 0 6 0 8 0 1 0 0
Model Cycle Time
log2 (CPM)
solute carrier family 24 member 4
SLC24A4
●
●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●●
●●
●
●
●●
●
●
●
●●●
●
●
●
●
●
●
●
●●
●
●
●●
●
●●●●
●
●●●
●
●
●
●●
●●
●
●
●
●
●
●
●
●●
●●●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●●
●●
●
●
●
●
●
●
●
●●
●●
●
●
●
●
●●●
●
●●
●
●
●●●
●
●●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●●●
●
●
●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●●
●●
●
●
●
●
●
●
●
●
●
●●
●
●
●●
●
●●
●
●
●
●
−5.0
−2.5
0.0
2.5
5.0
02 0 4 0 6 0 8 0 1 0 0
Model Cycle Time
log2 (CPM)
serpin family A member 4
SERPINA4
●
●●●●
●
●
●
●
●●
●●
●
●●
●
●
●
●●
●
●
●
●
●●
●
●●
●
●
●●
●
●●
●
●
●
●
●●
●●●
●●
●
●
●
●
●
●
●●
●
●
●●
●
●●
●
●
●●
●
●
●●
●
●
●
●
●
●●
●
●
●●
●●
●
●
●
●
●
●●
●
●
●
●
●●●
●
●
●
●
●●●
●
●
●●●●●●●●
●
●●
●
●●●
●
●
●
●
●
●●
●●●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●●
●
●●
●
●●
●
●
●
●
●
●●
●
●●●●
●
●
●
●
●
●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●●●●●
●
●
●●
●
●●
●●●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●●●●●
●
●
4
6
8
02 0 4 0 6 0 8 0 1 0 0
Model Cycle Time
log2 (CPM)
OCIA domain containing 2
OCIAD2
●
●
●
●
●
●
●
●●●●
●
●●
●
●
●
●
●
●
●
●
●
●●●●●
●●●●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●●●
●
●
●●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●●●
●
●
●
●
●
●●●
●
●●
●
●●
●
●
●
●
●●
●
●
●
●
●●●
●●
●
●
●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●●
●
●
●
●
●●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●●●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●●●●●
●
●●
●
●
●
●
●
●
●●
−5.0
−2.5
0.0
2.5
02 0 4 0 6 0 8 0 1 0 0
Model Cycle Time
log2 (CPM)
shisa family member 3
SHISA3
●●
●
●
●
●
●
●
●
●●
●
●
●
●●
●
●●
●
●
●
●●
●●
●
●
●
●
●
●
●
●
●●●
●
●
●
●
●
●
●
●
●
●●
●
●●
●
●
●
●
●
●
●
●●
●●
●●
●
●●
●
●●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●●●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●●
●●
●
●
●●
●
●
●●
●
●
●
●●
●
●
●
●
●
●
●●
●
●
●
●
●
●●●●●
●
●
●
●●
●●
●
●
●●
●●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●●
●
●●
●
●
●
●
●
●
●
●●●●
●
●
●
●●
●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●●●
●
●
●●●●
●
●
●
●
●
−3
0
3
6
02 0 4 0 6 0 8 0 1 0 0
Model Cycle Time
log2 (CPM)
LHFPL tetraspan subfamily member 3
LHFPL3
a
b
Fig. 6 | Changing gene expression across the cycle. Figure 6a. Genes from RNA-
seq analysis that significantly change expression (Padj <0.05) over 3.4% of the cycle
(approximately equal to a 24-hour window) at different stages of the menstrual
cycle. Statistical analysis used two-sided empirical Bayes moderated t-tests, tested
relative to a threshold fold change of 1.2, implemented in limma
35 with corrections
for multiple comparisons performed using the Benjamini-Hochberg method. Using
adjusted P values, 488 unique genes significantly change expression during
menstruation, 44 during the proliferative phase, and 2921 during the secretory
phase. Peak times of rapid change in gene expression approximately correspond to
menstrual (3% of the way through the cycle), late proliferative (51%), POD3 (66%),
POD5 (71%), POD11 (94%) and POD13 (98%). Figure 6b. Examples of 12 genes that
change expression significantly at different times across the cycle. Source data are
provided as a Source Data file.
Article https://doi.org/10.1038/s41467-023-41979-z
Nature Communications| (2023) 14:6222 7
poor, ranging from 18 to 40%, rising to 60–68% for POD dating within
1 day and to 65–81% for dating within 2 days15. The same study reported
that accuracy and interobserver reproducibility were unaffected by
refresher training, suggesting the limits of histological dating have
been reached. By contrast, our model reliably identi fies signi ficant
differences in gene expression from day to day through the cycle (see
Fig. 6)a n dc a na s s i g nat i m ep o i n to nas c a l em o r er efined than a single
day. In addition, our model accurately determines endometrial stage
regardless of cycle length and does not rely on measurement of
external events such as changes in circulating hormones orfixed time
points such as ultrasound detection of ovulation or first day of men-
strual bleeding. How accurate these results are in terms of cycle stage
is more difficult to ascertain. The comparison of NGS versus Illumina
data from the same samples (Fig. 3b) shows that at times of the cycle
where gene expression is changing rapidly, such as menstruation and
implantation (see Fig.6a), there is very high concordance between the
2 different methods. This agreement suggests that results at these
stages of the cycle are accurate to within less than 3-4% or approxi-
mately 1 day. At other times, such as during most of the proliferative
phase, there is less agreement, suggesting an accuracy of 7–11% or 2-3
days. By contrast, pathology dating during the proliferative phase is
limited to early, mid or late. The use of percentage to describe distance
through the cycle removes dif ficulties in comparing different length
cycles. The human menstrual cycle is inherently variable in overall
length and in the relative lengths of the menstrual, proliferative and
secretory phases
10–13. This variability makes accurate comparisons
between endometrial samples challenging. For example, 10 days post
last menstrual period (LMP) can be proliferative endometrium in one
woman and secretory in another. Similarly, POD6 can be pre-receptive
in one case and post-receptive in another. By de fining cycle stage in
terms of the percentage of the way through the menstrual cycle that
the endometrium is (which still allows precise mapping to events such
as menstruation and ovulation), it is possible to accurately synchronise
all endometrial samples to their own well-de fined point in the cycle.
From our data, menstrual is 0 –8% of the cycle, proliferative is 8 –58%
and secretory is 58–100%.
Another advantage of the molecular staging model is that it
generates mathematically de fined curves (splines) for expression
levels of all genes across the cycle. This allows gene expression from
different stages of the cycle to be readily normalised for comparison
studies. It also provides information on expression profiles for every
gene and can be used across different gene array platforms. As part of
the validation studies undertaken for the molecular staging model, we
used the same methods to derive a similar staging model using Illu-
mina HT-12 data from 198 samples that also had RNA-seq data. This
e x e r c i s ed e m o n s t r a t e dah i g hd e g r e eo fa g r e e m e n ti nt h et w om o d e l s’
cycle stage between the RNA-seq and Illumina HT-12 data sets, con-
firming that the molecular staging model method gives robust out-
comes when used with different transcriptomics platforms.
Endometrial expression data in GEO and other databases can now
be precisely staged and reinterp reted using this methodology. To
support this claim, we have reinterpreted 3 published data sets and
generated PCA plots to demonstrate that the molecular model pro-
vides accurate timing for all samples in each data set. We show good
agreement between the molecular staging model and most of the
published staging but identify several samples in each published data
set that are misclassi fied. With relatively small sample sizes, a few
misclassified samples can severely reduce the ability to identify dif-
ferentially expressed genes. Based on these results, we suggest the
need for caution in interpreting endometrial gene expression studies
with small sample sizes and/or where the authors have not been able to
incorporate accurate cycle staging into their analyses.
This day-to-day volatility in gene expression during the secretory
phase provides a compelling explanation for the consistent lack of
agreement between studies comparing gene expression where
comparison groups contain all stages of the secretory phase com-
bined. Use of the molecular staging model to accurately assign cycle
stage and account for these rapid changes in gene expression data will
significantly enhance our ability to compare endometrial gene
expression between experimental groups and identify real differences
due to phenotype or treatment.
The molecular staging model reveals remarkably synchronised
and dramatic changes in expression for many genes at different stages
of the cycle. The model delivers highly accurate cycle staging because
most of the variance in endometrial gene expression is due to cycle
stage, and not other factors. Over 50% of the gene models had an
adjusted R-squared value greater than 0.19 and 30% were over 0.36,
indicating that the proportion of variance in gene expression
explained by modelling cycle time is considerable. Based on this var-
iance, we have been able to identify 3453 genes that change expression
significantly within a 3.4% window of the endometrial cycle (which
approximately equates to a 24-hour period in a 28-day cycle). Of these,
488 were during menstruation, 44 during the proliferative phase, and
2921 during the secretory phase. These changes are remarkably con-
sistent between samples and synchronised between the 3453 genes
within each sample (see for example Fig.6 and Supplementary Fig 5). It
is these cycle-driven changes that provide the statistical strength and
remarkable chronological consistency of our modelling across the
cycle. If very few genes changed expression across the cycle, most of
the expression curves would be relatively flat and there would be no
obvious global minimum in the mean squared error plots. Instead, the
highly dynamic patterns of gene expression ensure that the mean
squared error plots give robust results for predicting a time in the
cycle. It is also important to reiterate that the mean squared error
solution is simultaneous for all 20,067 genes, giving con fidence that
significant changes in expression for individual genes are real.
Because most of the variance in endometrial gene expression is
linked to cycle stage, cycle stage determination using the molecular
staging model is unlikely to be in fluenced by other variables that
impact on gene expression during the menstrual cycle. This is because
even if a pathological condition changed the expression 100 ’so f
genes, it will only minimally influence the accuracy of the model which
is based on significant changes in the expression of several 1000’so f
genes. The overwhelming contribution of cycle stage variance to the
accuracy of the molecular model provides the rationale for the gen-
eralisability of the model to endometrial samples taken from women
with a wide range of uterine phenotypes and pathologies during the
menstrual cycle.
The Endometrial Receptivity Analysis (ERA) is a clinical test used
to determine if the endometrium is at the optimal stage of the men-
strual cycle for embryo implantation. The test was originally based on a
group of 238 genes that change expression rapidly around the time of
implantation
22, thus allowing the relative expression pro files to be
equated with ‘pre-receptive’, ‘receptive’ and ‘post-receptive’ endome-
trium. We correctly predicted that using the molecular staging model
to identify endometrial genes with rapidly changing expression pro-
files would also identify many of the ERA genes. As well as identifying
70% of the ERA genes we also found many hundreds of additional
genes that change expression signi ficantly around the time of
implantation, reinforcing how complex endometrial function is at this
critical stage of the menstrual cycle.
Using the molecular staging model to normalise gene expression
across the cycle we were able to identify, for the first time, numerous
endometrial genes that change expression significantly with increasing
age (Supplementary Table 3). Importantly, when we performed a meta-
analysis by adding a published dataset from another laboratory and
analysing it using our method, the number of signi ficant genes was
increased; a result that supports our original finding. An interesting
finding from this study was that when examining menstrual, pro-
liferative and secretory stages separately, nearly all the genes that
Article https://doi.org/10.1038/s41467-023-41979-z
Nature Communications| (2023) 14:6222 8
showed DGE with increasing age did so in endometrial samples taken
in the secretory phase. This suggests that it is the endometrial
response to progesterone that changes with increasing age, or at least
genes with biological actions related to secretory phase events such as
implantation, early pregnancy or the late secretory transition to
menstruation. One explanation for not identifying significant changes
in gene expression during the menstrual phase is low subject numbers
(total N = 28 with 0 genes signi ficant). But for the proliferative phase
from N = 165 subjects we found only 4 signi ficant genes while for the
secretory phase withN = 160 subjects there were 218 significant genes.
These data highlight the fundamental difference between proliferative
and secretory endometrium and bring in to question the biological
relevance of comparing these two very different tissues.
The most signi ficant up and down regulated GO pathways with
increasing age have strong biological plausibility (Supplementary
Data 2). Upregulated pathways include axoneme assembly, micro-
tubule bundle formation and cilium assembly. Microtubules are pre-
sent in nearly all cells and play a prominent role, along with ciliated
cells, in well characterised changes to the endometrial luminal epi-
thelium during the secretory phase. Although the functional relevance
is not clear, the transition from luminal microvilli to pinopodes around
the expected time of embryo implantation is thought to be a key
morphological marker for endometrial receptivity
26.M o r er e c e n te v i -
dence has found an increased incidence of luminal epithelial micro-
tubule abnormalities in the endometrium of women with
adenomyosis, possibly due to increased inflammation
27. Blood vessel
development and regulation of endothelial cell proliferation were the
top 2 down regulated GO biological pathways. The endometrium is a
site of rapid cyclical angiogenesis, vascular development and break-
down, with spiral arterioles developing their classical shape and other
endometrial vessels growing to match the increasing thickness of the
endometrium
28. Age-related alterations in endometrial angiogenic or
vascular potential could be expected to have profound effects on
endometrial function, and in particular fertility and subsequent
placentation.
Future endometrial DGE studies should account for age when
performing statistical modelling. It seems likely that use of the mole-
cular staging model will identify other gynaecological conditions that
influence endometrial DGE and hence need to be accounted for in
subsequent analyses.
There are limitations to the current study that need to be con-
sidered. The 236 endometrial samples used to develop the menstrual
cycle staging model in analysis 2 were nearly all taken from women
undergoing diagnostic laparoscopy to investigate pelvic pain and
suspected endometriosis (IVF subjects with secretory stage endome-
trium that were used in analysis 1 were not included in thefinal model).
Of the 236 women, 173 were diagnosed with endometriosis. Whether
this population can be considered as having normal endometrium is
open for debate. However, as reported in Supplementary Tables 1 & 2,
60 of the 236 women had delivered 1 or more live births, with 38 of
these 60 also being diagnosed with endometriosis. The ability to carry
a normal pregnancy to a successful outcome is evidence that the
endometrium was functionally‘normal’ in these women, and suggests
that if there are pathological differences in the endometrium of
women with endometriosis, they are either not universal, or not severe
enough to preclude normal endometrial function. Further, the 140
women who had not been pregnant were on average 6.2 years younger
than those that had had a successful pregnancy (29.2 vs 35.4 years of
age), suggesting that many had not yet tried to conceive and could well
have ‘normal’ endometrium. While it is possible that gynaecological
pathologies that affect endometrial gene expression may exist in some
of our 236 subjects, the conclusion from this work, as discussed earlier,
is that the predominant driver of changing endometrial gene expres-
sion is menstrual cycle stage, and n ot phenotypic factors. If pheno-
types that do influence gene expression are identified in the future, for
example age and ethnicity in the current work, then these can be
corrected for in the model and this will reduce variability and increase
sensitivity of future studies.
The relative lengths of the menstrual, proliferative and secretory
phases in the molecular staging model are determined by the dis-
tribution of women at each stage of the menstrual cycle. We made the
assumption that our population of 236 women were randomly dis-
tributed across the cycle although it is possible that this was not the
case. Hence some phases in the molecular staging model could be
compressed or extended over time compared to the normal popula-
tion. This does not alter the validity of the model for the relative
position of the samples in chronological order across the cycle but
does mean that mapping of molecular staging model onto cycle stages
may update as larger studies are performed.
In conclusion, we have developed and validated a method for
accurately determining endometrial cycle stage based on global gene
expression. Our ‘molecular staging model ’ reveals signi ficant and
remarkably synchronised daily changes in expression for over 3400
endometrial genes at different stages of the cycle, with most change
occurring during the secretory phase. These major day-to-day differ-
ences in endometrial gene expression provide a compelling explana-
tion for the failure of studies that lack accurate cycle staging to reach
consensus on genes of interest. Our study supports selected previous
findings and significantly extends existing data. Using the molecular
staging model to normalise expression data we can demonstrate sig-
nificant changes in endometrial gene expression with increasing age.
The molecular staging model provides a wealth of new data on
endometrial gene expression and establishes a method for investi-
gating the role of the endometrium in critical biological events such as
uterine receptivity for embryo implantation as well as gynaecological
pathologies such as endometriosis and menstrual disorders.
Methods
Subject recruitment and tissue collection
This work was approved by the Human Research Ethics Committee of
the Royal Women’s Hospital, Melbourne, Australia (Projects 11–24 and
16–43 for Study 1) and Melbourne IVF (Project 13/17 for Study 2), and all
subjects gave written informed consent. Study subjects did not receive
compensation for participating. A total of 358 endometrial samples
were collected for this study, comprising 264 samples taken at the time
of diagnostic surgery from women reporting pelvic pain with sus-
pected endometriosis ( ‘Study 1 ’) and 94 samples from individuals
undergoing IVF (‘Study 2’).
Some of the 358 subjects have had data published as part of
previous studies investigating genetic regulation of endometrial gene
transcription. Specifically, 123 Illumina Human HT-12 v4.0 samples
were included in a 2017 study
18, 229 Illumina Human HT-12
v4.0 samples in a 2018 study 29, and 169 & 206 RNA sequencing sam-
ples in studies from 202019,30.
Endometrial tissue samples (collected by curette or Pipelle
biopsy) were obtained for gene expression analysis, along with blood
samples for DNA extraction and hormone assays, patient ques-
tionnaires, past and present clinical histories, pathologyfindings and
surgical notes. All subjects were premenopausal and free from hor-
mone treatment at the time of biopsy. Endometrial tissue samples
were split and either stored in RNA later (Life Technologies, Grand
Island, NY, USA) at 4
oC before being stored at -80 oC for total RNA
extraction, or formalin fix e da n dp r o c e s s e dr o u t i n e l yf o rh i s t o l o g i c a l
assessment.
Histological dating of endometrium
All 358 endometrial samples were routinely evaluated14 by at least one
experienced pathologist and allocated to one of seven menstrual cycle
stages. Menstrual cycle stage de finitions, assuming a standardised
28 day cycle, and numbers of samples in each group were as follows:
Article https://doi.org/10.1038/s41467-023-41979-z
Nature Communications| (2023) 14:6222 9
Stage 1 = menstrual (n =1 8 ,d a y s1–4), 2 = early proliferative (n = 5, days
5–7), 3 = mid proliferative ( n =1 0 4 , d a y s 8–11), 4 = late proliferative
(n =2 9 ,d a y s1 2–15, includes‘interval’), 5 = early secretory (n =6 4 ,d a y s
16–19, or post ovulation days 2 –5), 6 = mid secretory ( n =7 6 , d a y s
20–23, or post ovulation days 6 –9), 7 = late secretory ( n =4 0 , d a y s
24–28, or post ovulation days 10 –14). Twenty-two biopsies taken by
Pipelle did not have adequate tissue for reliable pathology reporting. If
the pathology report crossed two definitions (e.g. mid/late secretory),
then the later cycle stage was used. The exception to this was in the
secretory phase where post-ovulatory day (POD) range crossed 2 cycle
stages with the majority of the days in the earlier stage, i.e. POD 4–6=
early secretory. If the pathology report only recorded‘proliferative’ or
‘secretory’, then mid-proliferative or mid-secretory was assigned,
resulting in elevated numbers in these 2 cycle stages. A subset of
secretory stage samples (n = 164) underwent additional evaluation and
were assigned an individual post-ovulatory day by a further 1 or 2
pathologists working independently of each other and the previous
assessment. Endometrial samples for which the pathologist reported
any abnormalities or evidence of exogenous hormone had already
been excluded from the study.
RNA-seq and gene expression array data preparation
Of the 358 endometrial samples, 290 had Illumina Human HT-12 v4.0
performed and 266 underwent RNA-seq (198 samples had both tech-
niques performed). Total RNA was isolated from endometrial samples
using the Allprep DNA/RNA Mini Kit (Qiagen, CA) as per the manu-
facturer’si n s t r u c t i o n s
29. Brie fly, RNA quality was checked using a
Bioanalyzer 2100 (Agilent Technologies, CA) and RNA concentration
was measured using a NanoDropND-6000 (Thermo Fisher Scientific,
USA). All samples were high quality with an RNA integrity number
greater than 8. Expression profiles in endometrial tissue were gener-
ated by hybridizing 750 ng of cRNA to Illumina Human HT-12 v4.0
Beadchips.
RNA samples were treated with Turbo DNA-free kit (Thermo
Fisher Scientific, USA) prior to RNA-seq library generation
19. Stranded
RNA-seq libraries were prepared using the Illumina TruSeq Stranded
Total RNA Gold protocol which includes ribosomal depletion
(Illumina, USA).
RNA-seq data processing
Raw sequencing reads were quality checked using FastQC v0.11.7 and
MultiQC v1.6. Low quality reads and contaminating HiSeq Illumina
adapter sequences were trimmed using Trimmomatic v0.36
31.T r i m -
med reads were aligned against the human reference genome
(Ensembl Homo sapiens GRCh38 release 84) using HISAT2 v2.0.5
32.
Transcript assembly was performed using StringTie v1.3.1 and the
EnsemblHomo sapiensGRCh38 release 91 reference annotation. Reads
mapping to each known transcript were quanti fied in StringTie
33 to
generate transcript-, exon- and intron-level expression matrices for
each individual. Raw gene count matrices were also produced using a
Python script provided by StringTie.
Normalisation of RNA-Seq and array expression values
Genes expressed at a low level by RNA-seq, i.e. genes with counts per
million (CPM) 80% of the samples, were removed. Raw gene
counts were normalized for composition bias and total raw reads
(library size) using the Trimmed Mean of M (TMM) method in the
edgeR R package
34. Normalized counts were converted to CPM and
log2 transformed (log2-CPM). Batch effects from sequencing were
removed using the ComBat function from the sva R package (Leek
et al., 2020). To load and normalise the microarray data, the R
packages limma
35 and lumi36 were used. Background correction and
robust spline normalization (RSN) produced logged values of probe
intensity and ComBat was used to remove microarray batch effects.
For probes to be included in the array analysis, annotation probe
quality was required to be “Good” or “Perfect” and detection p-
value < 0.05 in at least 20% of samples.
Genotyping
For determination of ancestry, DNA samples from each of the 358
individuals were genotyped on HumanCoreExome or In finium Psy-
chArray chips (Illumina, USA)
19. Quality control (QC) was performed in
PLINK29. Following QC, a total of 282,625 SNPs (hg19) were phased
using Shapelt V2 and taken forward to imputation using the haplotype
Reference
consortium reference panel (version r1.1 2016) on the
Michigan Imputation Server. SNPs with low imputation quality
(R
2 5%, minor allele frequency (MAF) < 1×10–4,a n d
Hardy–Weinberg equilibrium<1×10−6 after imputation were removed.
The remaining SNP positions were lifted over to the Ensembl genome
build 38 (GRCh38) using CrossMap v.0.2.8. SNPs failing to lift-over
were assigned to their new GRCh38 position manually based on
dbSNP151 GRCh38 patch release 7 (GRCh38.p7), leaving 6,230,993
SNPs for further analysis.
Hormone assays
Estradiol and progesterone concentrations were measured in bloods
taken at the time of endometrial sampling. Some of these hormone
data have been published previously
37. An additional 28 bloods were
assayed for progesterone (Serum P was tested on the Roche Cobas
e601 immunoanalyser, utilising electrochemiluminescence (ECLIA).
The lower limit of detection was 0.06 ng/mL. The inter-assay at a target
mean of 1.4 ng/mL returned a CV% of 3.7. The intra-assay at a target
mean of 1.5 ng/mL returned a CV% of 6.5). This gave a total of 187
progesterone results and 159 estradiol results that could be plotted
against the molecular staging model cycle stage.
Analysis 1
Development of the‘molecular staging model’ to assign cycle stage for
secretory stage samples only: All secretory stage samples with RNA-seq
data from Study 1 and Study 2 where 2 or 3 independent pathology
assessments agreed on the post-ovulatory day (POD) to within 2 days
(n = 96 of a possible 180 secretory stage samples) were selected for
analysis 1. For each gene, batch-corrected expression was fitt oP O D
using a penalised cubic regression spline using 3 knots implemented
with the generalised additive model (gam) function from the mgcv R
package
38. Each curve was used to obtain the expected expression
value for each gene for any given day. For each sample, an estimated
POD was obtained using the day which minimised the mean squared
error (MSE) between the observed expression and the expected
expression across all genes.
Alternatively, this procedure can be described minimisingd in the
loss function:
LdðÞ =
X
g2G
yg /C0 f g dðÞ
/C16/C17 2
ð1Þ
where d is the POD constrained between 1 and 14 days, g is a gene in
gene set G,y g is the observed expression of gene g,a n d fg(d) is the
spline function that describes the expected expression of gene g for
day d. Additionally, K-fold cross-validation where K = 5 was performed
to ensure the model was not overfitting.
To illustrate that using larger, less precise, units of time can be
used to estimate cycle time using the same method, an additional
model was built using the pathology-assigned 3 secretory cycle stages
(i.e. stages 5, 6, and 7 corresponding to early-, mid-, and late-secretory
respectively) instead of the pathology-assigned POD (i.e. 1 –14 days).
Using the RNA-seq batch-corrected expression data, each gene wasfit
using the same penalised cubic regression spline (k = 3) as a function of
stage, and curves were generated for each gene. Cycle time was esti-
mated for each sample by selecting the time point in the stage (from
Article https://doi.org/10.1038/s41467-023-41979-z
Nature Communications| (2023) 14:6222 10
4.5 to 7.5) that minimised the MSE between the observed expression
and the expected expression across all genes. As validation, cycle time
estimated from the 3 stages model was compared to the cycle time
estimated from the 14-day POD model.
Analysis 2
Molecular staging model using 7 pathology stages for the whole cycle
with RNA-seq and array expression data: The method for developing
the POD of cycle prediction model was replicated with some mod-
ifications using all samples from Study 1 ( N = 236 for NGS) classi fied
into 7 cycle stages by histopathology. Because the majority of the
proliferative phase samples were not assigned as early, mid or late by
the pathologist, we re-assigned all proliferative samples into early, mid
or late by fit t i n gap e n a l i s e dc u b i cr e g r e s s i o ns p l i n e( k = 3 )u s i n g
menstrual, proliferative, and early secretory gene expression data. A
proliferative time point was estimated using the time-point which
minimised the mean squared error between the observed expression
and the expected expression across all genes. The proliferative sam-
ples were then split into equal groups of early, mid, and late using this
time point. This approach assumes that patients presenting for surgery
in a public hospital system will approximate a uniform distribution
across the menstrual cycle, so the number of early, mid, and late
proliferative samples will be approximately equal. Once the pro-
liferative samples were assigned to early, mid and late stages, a pena-
lised cyclic cubic regression spline (k = 8) was fit using the 7 stages.
Each endometrial sample is then assigned a‘day’ from the model using
the time which minimises the mean squared error between the
observed expression data for all genes and their corresponding gene
models. This is equivalent to minimising the loss function in Eq. ( 1),
except now d can represent any timepoint in the cycle. This ‘day’ is a
relative timepoint in the cycle and does not correspond to a real day.
Continuing with the assumption that all 236 women were approxi-
mately uniformly distributed across the cycle, the data are trans-
formed so that the distance in time between each sample is identical.
This process in effect ranks all the samples in order from the start to
the end of the cycle, and no longer relies on assigning days from an
idealised 28-day cycle. The x-axis was then scored from 0–100 so that
the individual scores for each sample represented the percentage of
the way through the menstrual cycle that each sample was.
The gene curves were then re fitted using the new derived cycle
times for each sample with a penalized cyclic cubic regression spline
(k = 30). For visualisation purposes, normalisation of gene expression
for cycle stage was then derived by subtracting the expected expres-
sion from the observed expression (i.e. calculating the residuals) and
re-adding the mean.
Validation of the molecular staging model
Several studies were conducted t o assess how endometrial gene
expression data that were generated using the molecular staging
model performed.
The first validation study as part of Analysis 1 was to compare
Results
from the secretory stage that had daily POD pathology, and
hence substantially more accurate cycle stage information, with results
that only used 3 pathology stages across the secretory phase. The
similar results from the 2 different pathology inputs validated the use
of pathology data dividing the endometrium into 7 cycle stages to
develop the molecular staging model across the whole menstrual cycle.
The second validation study was to repeat Analyses 1 and 2 using
Illumina HT-12 data and then compare results for the 198 out of
358 samples that had both RNA-seq and Illumina HT-12 data.
The third validation study was to repeat Analysis 2 using unsu-
pervised methods (i.e. without the use of pathology cycle stage
information). K-means clustering (k = 5) was performed using thefirst
2 dimensions after dimensionality reduction principal component
analysis (PCA) (Supplementary Fig 1). Clusters were ordered by their
distance between each other, and for each gene, expression wasfitted
using a penalised cyclic cubic regression spline (k = 6) using the
ordered cluster number as a proxy for cycle time. Identical to Analysis
2, for each sample, an optimal cycle time was obtained by calculating
the time that minimised the MSE, and subsequently all times were
transformed to be uniformly distributed on a scale from 0 to 100.
Cycle time estimates from the validation model were compared to the
Analysis 2 cycle times, however, since time 0 for the validation model
was not expected to correspond to the start of menstruation, time was
offset by the earliest menstrual sample (identified by Analysis 2) before
assessing performance (Supplementary Fig 2).
The fourth validation study was to determine whether indepen-
dent endocrine data supported the molecular staging model. For a
subset of samples where endocrine data were available, peripheral
blood estradiol (n = 159) and progesterone (n = 187) levels were plotted
against the cycle stage from the molecular staging model.
The molecular staging model was also used to re-analyse 3 pub-
lished datasets available in GEO (GEO DataSets ID: GSE65099
21,
GSE14154920 and GSE18048517). Principal component analysis (PCA)
plots from these data sets were replotted with the unaltered cycle
stage from the original publication and our molecular staging model
cycle stage for comparison.
Application of the molecular staging model
The molecular staging model for gene expression across the menstrual
cycle was applied to 3 questions: Does (1) age or (2) ancestry have any
influence on endometrial gene expression, and (3) at what stages of the
cycle does gene expression change most rapidly? Differential expres-
sion analysis was performed with the predicted cycle time as an
additional factor in the linear model, modelled as a regression spline.
Empirical Bayes moderated t-tests, implemented in the R limma
package, were used to assess if genes were differentially expressed.
Age of patient at biopsy was analysed as a continuous variable
using all subjects with RNA-seq data ( N = 266). Since differential
expression effects due to age were identified, subsequent differential
gene expression (DGE) analyses included age as a factor when fitting
the linear model. To further explore the effects of age on endometrial
gene expression n = 87 subjects from GSE141549 were also analysed,
and a meta-analysis was performed using the weighted Fisher ’s
Method
for combining p-values implemented in the metapro R
package
39, where weights were proportional to each study ’ss a m p l e
size. Ensembl ID’s from the current data were matched with the Illu-
mina probe ID’s from GSE141549 resulting in 12,868 genes in common
between the 2 data sets. If multiple probes matched to the same
Ensembl ID, the probe with the greatest mean expression was used.
Analyses were run for the whole menstrual cycle, and separately for the
menstrual, proliferative and secretory phases. After multiple hypoth-
esis correction, genes with false discovery rate (FDR) corrected
P < 0.05 were considered to be differentially expressed. Gene ontology
enrichment analysis was performed using the clusterPro filer R
package
40.
Ancestry of subjects as defined by a previous study19 was used to
look for differences in endometrial gene expression using pairwise
comparisons of each ancestry group. Genetic ancestry of subjects was
assigned using principal component analysis (PCA). Briefly, genotype
data from participants was merged with the 1000 Genomes P3v5
Reference
data using markers common to both cohorts. Population
clusters were determined using the first 5 PCs and were annotated
according to the five 1000 genomes super populations (European,
Eastern Asian, Admixed America, African, Southern Asian).
T oi n v e s t i g a t ec h a n g i n gg e n ee x p r e s s i o na c r o s st h ec y c l ea l l
samples with RNA-seq data from Study 1 ( n = 236) were ranked in
chronological order from start to end of the cycle. A ‘sliding window’
approach was then used to compare DGE between samples 1–8v e r s u s
samples 9–16, followed by samples 2 –9 versus 10 –17, then 3 –10 vs
Article https://doi.org/10.1038/s41467-023-41979-z
Nature Communications| (2023) 14:6222 11
11–18 and so on for all samples across the menstrual cycle. Group sizes
were arbitrarily set at 8 because this represents 3.4% of the 236 samples
or approximately 1 day assuming a mean cycle length of 28 days.
Moderated t-tests with a fold-change cut-off of 1.2, as implemented in
limma’s treat function, were used to identify differentially expressed
genes with P
adj <0 . 0 5a te a c hw i n d o w .
To validate a subset of genes we identi fied as having rapidly
changing expression around the time of embryo implantation, we ran a
comparison with genes identified in the original Endometrial Recep-
tivity Analysis (ERA) publication
22. The ERA identi fied 238 genes that
show major changes in expression before, during and after the time of
embryo implantation at POD 3–7. We identified all endometrial genes
that showed significant changes in expression within any given 24 hr
period at the same time of the menstrual cycle to confirm that our list
contained a high proportion of the ERA genes.
Reporting summary
Further information on research design is available in the Nature
Portfolio Reporting Summary linked to this article.
Data availability
The RNA-seq data generated in this study have been deposited in the
Gene Expression Omnibus database under accession number
GSE234354 and the Illumina HT-12 microarray data under accession
number GSE234368. Additional public datasets used are available on
GEO with accession numbersGSE65099 (PMID: 26418742),GSE141549
(PMID: 32859947) and GSE180485 (PMID: 35092277). Homo sapiens
Reference
genome GRCh38 and genome annotation were obtained
from Ensembl (release 91). Source data are provided with this paper.
Code availability
The code for our model is available as an R package at https://github.
com/jessicachung/endest (DOI: 10.5281/zenodo.8321573) and an R
Shiny application is available at https://github.com/jessicachung/
endspect. Analysis scripts for this manuscript can be found at
https://github.com/jessicachung/endo_model_paper.E n d e s ti sp u b -
lished pursuant to the terms located in the applicable repository at
github.com which terms permit reproduction, publication and adap-
tation of endest solely for non-commercial purposes. Publication of
endest at github.com is not subject to the publication terms applicable
to Nature publications.
References
1. Jain, V. et al. Uterine bleeding: how understanding endometrial
physiology underpins menstrual health.Nat. Rev. Endocrinol. 18,
290–308 (2022).
2. Darba, J. & Marsa, A. Economic Implications of Endometriosis: A
Review. Pharmacoeconomics40, 1143–1158 (2022).
3. Hoppes, E. et al. Global research and learning agenda for building
evidence on contraceptive-induced menstrual changes for
research, product development, policies, and programs.Gates
Open Res 6, 49 (2022).
4. Ponnampalam, A. P. et al. Molecular classi fication of human endo-
metrial cycle stages by transcriptional profiling. Mol. Hum. Reprod.
10,8 7 9–893 (2004).
5. Kao, L. C. et al. Expression pro fil i n go fe n d o m e t r i u mf r o mw o m e n
with endometriosis reveals candidate genes for disease-based
implantation failure and infertility.Endocrinology144,
2870–2881 (2003).
6. Girling, J. E. et al. Differential Gene Expression in Menstrual Endo-
metrium From Women With Self-Reported Heavy Menstrual
Bleeding. Reprod. Sci. 24,2 8–46 (2017).
7. Koot, Y. E. et al. An endometrial gene expression signature accu-
rately predicts recurrent implantation failure after IVF.Sci. Rep. 6,
19411 (2016).
8. Aghajanova, L. et al. Stanniocalcin-1 expression in normal human
endometrium and dysregulation in endometriosis.Fertil. Steril.106,
681–691.e1 (2016).
9. Aghajanova, L. et al. Effects of noncavity-distorting fibroids on
endometrial gene expression and function.Biol. Reprod. 97,
564–576 (2017).
10. Najmabadi, S. et al. Menstrual bleeding, cycle length, and follicular
and luteal phase lengths in women without known subfertility: A
pooled analysis of three cohorts.Paediatr. Perinat. Epidemiol.34,
318–327 (2020).
11. Soumpasis, I., Grace, B. & Johnson, S. Real-life insights on men-
strual cycles and ovulation using big data.Hum. Reprod. Open
2020, hoaa011 (2020).
12. Bull, J. R. et al. Real-world menstrual cycle characteristics of more
than 600,000 menstrual cycles.NPJ Digit Med 2,8 3( 2 0 1 9 ) .
13. Tatsumi, T. et al. Age-Dependent and Seasonal Changes in Men-
strual Cycle Length and Body Temperature Based on Big Data.
Obstet Gynecol. 136, 666–674 (2020).
1 4 . N o y e s ,R .W . ,H e r t i g ,A .T .&R o c k ,J .D a t i n gt h eE n d o m e t r i a lB i o p s y .
Fertil. Steril.
1,3 –25 (1950).
15. Duggan, M. A. et al. The accuracy and interobserver reproducibility
of endometrial dating.Pathology 33,2 9 2–297 (2001).
16. Ruiz-Alonso, M. et al. Endometrial Receptivity Analysis (ERA): data
versus opinions.Hum. Reprod. Open 2021, hoab011 (2021).
17. Lipecki, J. et al. EndoTime: non-categorical timing estimates for
luteal endometrium.Hum. Reprod. 37,7 4 7–761 (2022).
18. Fung, J. N. et al. The genetic regulation of transcription in human
endometrial tissue.Hum. Reprod. 32,8 9 3–904 (2017).
19. Mortlock, S. et al. Tissue speci fic regulation of transcription in
e n d o m e t r i u ma n da s s o c i a t i o nw i t hd i s e a s e .Hum. Reprod. 35,
377–393 (2020).
2 0 . G a b r i e l ,M .e ta l .Ar e l a t i o n a ld a tabase to identify differentially
expressed genes in the endometrium and endometriosis lesions.
Sci. Data 7, 284 (2020).
2 1 . L u c a s ,E .S .e ta l .L o s so fE n d o m e t rial Plasticity in Recurrent Preg-
nancy Loss. Stem Cells 34,3 4 6–356 (2016).
22. Diaz-Gimeno, P. et al. A genomic diagnostic tool for human endo-
metrial receptivity based on the transcriptomic signature.Fertil.
Steril. 95,5 0–6 0 ,( 2 0 1 1 ) .6 0 . e 1 - 1 5 .
23. Wang, A. et al. Health Care Costs and Treatment Patterns Asso-
ciated with Uterine Fibroids and Heavy Menstrual Bleeding: A
Claims Analysis.J. Women’sH e a l t h( L a r c h m t . )31,8 5 6–863 (2022).
24. Schoep, M. E. et al. Productivit y loss due to menstruation-related
symptoms: a nationwide cross-sectional survey among 32 748
women. BMJ Open 9, e026186 (2019).
25. Niederberger, C. et al. 25 historic papers: an ASRM 75th birthday
gift from Fertility and Sterility.Fertil. Steril. 112,e 2–e27 (2019).
26. Quinn, K. E. et al. Pinopodes: Recent advancements, current per-
spectives, and future directions.Mol. Cell Endocrinol. 501,
110644 (2020).
27. Khan, K. N. et al. An axonemal alteration in apical endometria of
human adenomyosis.Hum. Reprod. 36,1 5 7 4–1589 (2021).
28. Rogers, P. A. et al. Endometrial angiogenesis, vascular maturation,
and lymphangiogenesis.Reprod. Sci. 16,1 4 7
–151 (2009).
29. Fung, J. N. et al. Genetic regulation of disease risk and endometrial
gene expression highlights potential target genes for endome-
triosis and polycystic ovarian syndrome.Sci. Rep. 8, 11424 (2018).
30. Holdsworth-Carson, S. J. et al. Obesity does not alter endometrial
gene expression in women with endometriosis.Reprod. Biomed.
Online 41, 113–118 (2020).
31. Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trim-
mer for Illumina sequence data. Bioinformatics30,
2114–2120 (2014).
32. Kim, D. et al. Graph-based genome alignment and genotyping with
HISAT2 and HISAT-genotype.Nat. Biotechnol.37,9 0 7–915 (2019).
Article https://doi.org/10.1038/s41467-023-41979-z
Nature Communications| (2023) 14:6222 12
33. Kovaka, S. et al. Transcriptome assembly from long-read RNA-seq
alignments with StringTie2.Genome Biol. 20,2 7 8( 2 0 1 9 ) .
34. Robinson, M. D. & Oshlack, A. A sc aling normalization method for
differential expressionanalysis of RNA-seq data. Genome Biol. 11,
R25 (2010).
35. Ritchie, M. E. et al. limma powers differential expression analyses
for RNA-sequencing and microarray studies.Nucleic Acids Res. 43,
e47 (2015).
36. Du, P., Kibbe, W. A. & Lin, S. M. lumi: a pipeline for processing
Illumina microarray.Bioinformatics24,1 5 4 7–1548 (2008).
37. Marla, S. et al. Genetic risk factors for endometriosis near estrogen
receptor 1 and coexpression of genes in this region in endometrium.
Mol. Hum. Reprod. 27, gaaa082 (2021).
38. Wood, S. N., Pya, N. & Säfken, B. Smoothing Parameter and Model
Selection for General Smooth Models.J. Am. Stat. Assoc. 111,
1548–1563 (2016).
39. Yoon, S. et al. Powerful p-value combination methods to detect
incomplete association.Sci. Rep. 11,6 9 8 0( 2 0 2 1 ) .
40. Yu, G. et al. clusterPro filer: an R package for comparing biological
themes among gene clusters.OMICS 16,2 8 4–287 (2012).
Acknowledgements
We acknowledge the work of Research Nurses Ranita Charitra, Tracy
Middleton and Irene Bell who recruited all the participants in this study
and Leonie Cann who provided laboratory support. Thank you to the
surgeons and theatre booking staff at the Royal Women’s Hospital and
Staff at Melbourne IVF. We also acknowledge and thank the many
women who provided tissues and clinical data that were used in this
study. Research reported in this publication was supported in part by
National Health and Medical Research Council (NHMRC) project grants
GNT1012245 and GNT1105321 (PAWR, GWM, JEG, SJH), NHMRC Fel-
lowship GNT1078399 and Investigators Grant GNT1177194 (GWM),
NHMRC Medical Postgraduate Scholarship No. 1055814 and RANZCOG
Fellows’ Clinical Research Scholarship (2013) (WTT), funds from Mel-
bourne IVF, and the Medical Research Future Fund MRF1199715 (PAWR,
MH, SJH, JFD).
Author contributions
WTT, JC and PAWR were responsible for the conception and design of
the study. PAWR, WTT, SJH, JFD, MH and JEG acquired the tissues along
with clinical and surgical data. HR, SB and VO provided expert patho-
logical assessment of endometrial biopsies. JC and CS provided bioin-
formatics and statistical input. RK, JNF, SM and GWM undertook all the
gene expression assays and initial QC of the expression data. All authors
have read, provided input and approved the final manuscript.
Competing interests
The method described in this study is the subject of a patentfiled at the
Australian Patent office (21/6/23 (PCT/AU2023/050559).), by the Uni-
versity of Melbourne with WTT, JC and PAWR as inventors. The remaining
authors declare no competing interests.
Additional information
Supplementary informationThe online version contains
supplementary material available at
https://doi.org/10.1038/s41467-023-41979-z.
Correspondenceand requests for materials should be addressed to P.
A. W. Rogers.
Peer review informationNature Communicationsthanks Sascha Ott and
the other, anonymous, reviewer(s) for their contribution to the peer
review of this work. A peer review file is available.
Reprints and permissions informationis available at
http://www.nature.com/reprints
Publisher’s note Springer Nature remains neutral with regard to jur-
isdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons
Attribution 4.0 International License, which permits use, sharing,
adaptation, distribution and reproduction in any medium or format, as
long as you give appropriate credit to the original author(s) and the
source, provide a link to the Creative Commons licence, and indicate if
changes were made. The images or other third party material in this
article are included in the article’s Creative Commons licence, unless
indicated otherwise in a credit line to the material. If material is not
included in the article’s Creative Commons licence and your intended
use is not permitted by statutory regulation or exceeds the permitted
use, you will need to obtain permission directly from the copyright
holder. To view a copy of this licence, visithttp://creativecommons.org/
licenses/by/4.0/.
© The Author(s) 2023
Article https://doi.org/10.1038/s41467-023-41979-z
Nature Communications| (2023) 14:6222 13
Text is read by the "Ask this paper" AI Q&A widget below. Extraction quality varies by source — PMC NXML preserves structure cleanly, OA-HTML may include some navigation residue, and OA-PDF can have broken hyphenation. The publisher copy (via DOI) is the canonical version.
My notes (saved in your browser only)
Ask this paper
Answers must be backed by verbatim quotes from this paper's full text. Hallucinated quotes are dropped automatically; if no verbatim passage answers the question, we say so. How this works
Condition tags
MeSH descriptors
Citation neighborhood
Papers in the corpus that this work cites (lower rings, blue) and that cite this one (upper rings, green). Dot size scales with the paper's in-corpus citation count — bigger dot = more influential within the endo/adeno field. Click a dot to open that paper. [ expand to 2 hops ] — adds papers reached through this work's immediate citers/citees. Heavier; up to 60 extra dots.
References (41)
- An axonemal alteration in apical endometria of human adenomyosis via openalex
- A relational database to identify differentially expressed genes in the endometrium and endometriosis lesions via openalex
- Dating the Endometrial Biopsy via openalex
- Differential Gene Expression in Menstrual Endometrium From Women With Self-Reported Heavy Menstrual Bleeding via openalex
- Effects of noncavity-distorting fibroids on endometrial gene expression and function via openalex
- Endometrial angiogenesis, vascular maturation, and lymphangiogenesis via openalex
- Expression profiling of endometrium from women with endometriosis reveals candidate genes for disease-based implantation failure and infertility via openalex
- Genetic risk factors for endometriosis near estrogen receptor 1 and coexpression of genes in this region in endometrium via openalex
- Obesity does not alter endometrial gene expression in women with endometriosis via openalex
- Stanniocalcin-1 expression in normal human endometrium and dysregulation in endometriosis via openalex
- Tissue specific regulation of transcription in endometrium and association with disease via openalex
- Uterine bleeding: how understanding endometrial physiology underpins menstrual health via openalex
- W2593643174 via openalex
- W2884664850 via openalex
- W2954816350 via openalex
- W2964849163 via openalex
- W2970585464 via openalex
- W2980696159 via openalex
- W2980707312 via openalex
- W2984136994 via openalex
- W3007415140 via openalex
- W3016770067 via openalex
- W3086103533 via openalex
- W3137679227 via openalex
- W3154390398 via openalex
- W3203033666 via openalex
- W4205159168 via openalex
- W4224295413 via openalex
- W4308385954 via openalex
- W2995527562 via openalex
- W2035618305 via openalex
- W2069902854 via openalex
- W2106608030 via openalex
- W2107018762 via openalex
- W2110494659 via openalex
- W2131271579 via openalex
- W2146512944 via openalex
- W2152689989 via openalex
- W2159922867 via openalex
- W2400783197 via openalex
- W2586065198 via openalex
Cited by (4)
- NRF2 signalling in endometriosis: A hormone-independent axis for fertility-preserving therapeutic innovation 2026
- Endometriosis and adenomyosis unveiled through single-cell glasses 2025
- Deconvolution of bulk endometrial tissue identifies cell type proportions and expression signatures associated with endometrial function and disease 2025
- Regulation of RNA splicing in endometrial tissue and its association with endometriosis 2025
Source provenance
- europepmc
- last seen: 2026-10-13T06:10:10.580696+00:00
- openalex
- last seen: 2026-06-04T00:00:01.174412+00:00
- pubmed
- last seen: 2026-10-13T06:09:08.380565+00:00
License: CC0
· commercial use OK