Abstract
COVID-19 vaccines have been instrumental tools in the fight against
SARS-CoV-2 helping to reduce disease severity and mortality. At the
same time, just like any other therapeutic, COVID-19 vaccines were
associated with adverse events. Women have reported menstrual cycle
irregularity after receiving COVID-19 vaccines, and this led to renewed
fears concerning COVID-19 vaccines and their effects on fertility.
Herein we devised an informatics workflow to explore the causal drivers
of menstrual cycle irregularity in response to vaccination with mRNA
COVID-19 vaccine BNT162b2. Our methods relied on gene expression
analysis in response to vaccination, followed by network biology
analysis to derive testable hypotheses regarding the causal links
between BNT162b2 and menstrual cycle irregularity. Five high-confidence
transcription factors were identified as causal drivers of
BNT162b2-induced menstrual irregularity, namely: IRF1, STAT1, RelA (p65
NF-kB subunit), STAT2 and IRF3. Furthermore, some biomarkers of
menstrual irregularity, including TNF, IL6R, IL6ST, LIF, BIRC3, FGF2,
ARHGDIB, RPS3, RHOU, MIF, were identified as topological genes and
predicted as causal drivers of menstrual irregularity. Our
network-based mechanism reconstruction results indicated that BNT162b2
exerted biological effects similar to those resulting from prolactin
signaling. However, these effects were short-lived and didn’t raise
concerns about long-term infertility issues. This approach can be
applied to interrogate the functional links between drugs/vaccines and
other side effects.
Subject terms: Reproductive signs and symptoms, RNA vaccines
Introduction
Corona virus disease 19 (COVID-19) pandemic, caused by the severe acute
respiratory syndrome corona virus 2 (SARS-CoV-2), is still sweeping the
world, causing more fatalities and threatening with dangerous viral
variants and more economic losses^[32]1–[33]4. The virus has infected
hundreds of millions and caused millions of deaths worldwide. COVID-19
vaccines have been instrumental tools in the fight against the virus,
and they helped reduce disease severity and mortality^[34]5–[35]10. At
the same time, COVID-19 vaccines were associated with adverse events,
just like any other therapeutic^[36]11–[37]24.
The occurrence of post-vaccine menstrual cycle disturbances has gone
unnoticed during clinical trials phase of COVID-19 vaccines. Then,
thousands of reports pointing to menstrual changes started to surface
following worldwide vaccination campaigns^[38]14,[39]16,[40]25.
Healthcare authorities dismissed these claims at the beginning and
considered them unjustified. But such reports continued to appear from
various countries around the world which led to increased vaccine
hesitancy in young women^[41]26–[42]29.
Recently, serious concerns have been raised about the effects of
COVID-19 vaccines on menstruation^[43]27–[44]31, and these fears keep
escalating. Thousands of women reported post-COVID-19-vaccine menstrual
changes to health care authorities around the world, and several
published studies indicated an association between menstrual
abnormalities and COVID-19 vaccines^[45]32–[46]34. Women feared that
menstrual changes suggest long-term adverse effects on fertility and
pregnancy, which led to hesitation against vaccination among women. In
fact, menstrual changes have been reported after receiving both mRNA
and adenovirus vectored COVID-19 vaccines which led to the speculation
that it is the nature of the immune response to vaccines, rather than
vaccine components, that led to these adverse
events^[47]13,[48]14,[49]16,[50]31,[51]35–[52]37. Furthermore, previous
reports suggested that vaccine-associated period changes occur due to
transient perturbations to the hypothalamic-pituitary-ovarian (HPO)
axis^[53]38–[54]40. In fact, changes in menstrual cycles have been
reported for non-COVID-19 vaccines including the human papillomavirus
and typhoid vaccines^[55]41,[56]42.
Herein, we devised a workflow to assess menstrual adverse events in
response to treatment with mRNA vaccine BNT162b2. Our results revealed
a causal link implicating prolactin signaling and hormone-induced
effects on the menstrual cycle and endometrium resulting from
post-vaccine gene expression perturbances. Fortunately, gene expression
perturbations were short-term and therefore are not expected to cause
long-term menstrual irregularities. The approach devised and
implemented herein can be applied to assess other vaccines and other
vaccine-induced biological effects.
Results
Systems biology findings
We undertook a systems biology approach to derive transcriptional
signatures for COVID-19 mRNA vaccines relying on BNT162b2
transcriptomics data. Our workflow is shown in Fig. [57]1. Previously,
we applied similar approaches to explore the network pharmacology of
drugs and vaccines^[58]43, as well as investigating disease
pathogenesis pathways for prioritizing biomarkers and drug
targets^[59]11,[60]12. Each study workflow is tweaked to suit the
scientific questions we are asking as well as the types of data we
have. In this study, the analysis of transcriptional raw data extracted
from [61]GSE169159 indicated that gene expression alterations on day 22
of receiving the first vaccine dose (i.e., the day after receiving the
second vaccine dose) affected genes that are known biomarkers or drug
targets for menstrual cycle disturbances. All details on deriving gene
signatures from transcriptomics data of [62]GSE169159 are described
elsewhere^[63]12. Additionally, no significant DEGs were observed on
day 28 after receiving the first vaccine dose.
Fig. 1. Informatics systems biology workflow.
Fig. 1
[64]Open in a new tab
A devised workflow for studying the mechanism(s) underlying the
biological effects of vaccines.
Vaccine Gene Signatures
Two transcriptomic gene signatures (GS) for BNT162b2 vaccine were
derived from gene expression profiling experiments in response to
treatment with vaccine ([65]GSE169159)^[66]44. The first gene signature
(GS1) consisted of 1853 differentially expressed genes (DEGs), 884
upregulated and 969 downregulated DEGs that satisfied the
prioritization criteria of a false discovery rate (FDR) ≤ 0.05, and a
log2 fold change (log[2]FC) ≥ 2.00 or ≤ −2.00. The second gene
signature (GS2) consisted of 108 DEGs, 74 upregulated and 18
downregulated DEGs that satisfied the prioritization criteria of an
FDR ≤ 0.05, and a log[2]FC ≥ 5.00 or ≤ −5.00.
In order to get a better idea about the biological significance of the
DEGs in response to treatment with BNT162b2, we applied a
bioinformatics workflow relying on both downstream enrichment analysis,
and upstream analysis for putative regulators responsible for causing
the gene expression changes observed in transcriptomics data.
Upstream regulation analysis
Upstream analysis was performed using the DEGs in GS1 (applying a
threshold of log[2]FC) and GS2 using causal reasoning^[67]45. GS1
consisted of 1853 DEGs and therefore was trimmed by our upstream
analysis algorithm to reduce the complexity of generated results. The
algorithm automatically applied log[2]FC ≥ +2.62 or ≤ −2.62 threshold
which reduced the number of DEGs to 1107 (743 upregulated and 612
downregulated); since the applied causal reasoning algorithm requires a
query list of about thousand genes on average. This analysis resulted
in the prediction of 826 activated and 480 inactivated upstream
regulators including transcription factors, kinases, phosphatases, and
microRNAs (Supplementary Table [68]1). All prioritized upstream
regulators have prediction activities with p-values ≤ 0.05 and a
calculation distance = 1–3.
To focus on upstream regulators of genes with the maximal differential
expression in response to vaccination, we applied the same analysis to
GS2 consisting of 92 DEGs (74 upregulated and 18 downregulated) that
had log[2]FC ≥ +5.00 or ≤ −5.00. Our analysis resulted in predicting
625 activated proteins 389 inactivated (Supplementary Table [69]2).
All prioritized upstream regulators were scored based on the number of
differentially expressed genes that can be reached via the shortest
paths, and the correctness of the regulation. The activity prediction
correctness is assessed based on the activation and inhibition edges
along the paths and the expected and observed directionality of fold
changes of the DEGs. It should be noted that the calculation distance
is one of the most important parameters that can distinguish direct
regulation effects from indirect effects. For example, a calculation
distance of one means that the upstream regulator is one step away from
the transcriptional event indicating that the regulation event in
direct.
Filtering upstream regulators by molecular function and distance
The changes in gene expression that we observe in response to
perturbagens (e.g., vaccines), often result from interactions between
gene regulatory regions and regulatory proteins such as transcription
factors, kinases, phosphatase, RNA molecules and others. Transcription
factors (TFs) are considered topologically more important than DEGs
especially for purposes of mechanism reconstruction; or at least
complementary to DEGs for reconstructing the biological pathways and
networks responsible for a phenotype of interest (e.g., menstrual
irregularities in response to vaccination). Although transcription
factors are not the only master regulators, but they are probably one
of the easiest to validate experimentally using in vitro assays.
Hence, we sought filtering the predicted causal regulatory proteins by
“molecular function” = “transcription factor” and “calculation
distance” = 1. A calculation distance=1, indicates direct effects on
transcription. This filtering step resulted in 13 transcription factors
(IRF1, STAT1, STAT2, RELA, IRF9, SPI, NFKB1, IRF3, IRF7, BCL6, PRDM1,
GATA3) for DEGs in GS1 with log[2]FC ≥ +2.62 or ≤ −2.62 (Supplementary
Table [70]3). Applying the same filters on the causal regulatory
proteins predicted for GS2 resulted in the prediction of five
transcription factors (IRF1, STAT1, RELA, IRF3 and STAT 2), with a
calculation distance = 1 (Supplementary Table [71]4). The overlapping
TFs resulting for the previous two causal reasoning analyses are listed
in Table [72]1. Such hits can be considered higher confidence
transcription factors for causing the observed phenotype.
Table 1.
Top commonly predicted upstream regulatory hubs at a distance of 1 for
DEGs in GS1 and GS2.
Molecular entity Molecular function Gene Predicted activity^a
Correct/total^b Prediction p-value^c Distance^d Evidence^e
IRF1 Transcription factor IRF1 +
60/62*
17/17**
4.24E-16^f
7.63E-06^g
1 Yes^[73]59,[74]60
STAT1 Transcription factor STAT1 +
64/70*
15/16**
1.22E-13^f
2.59E-04^g
1 Yes^[75]118–[76]122
RelA (p65 NF-kB subunit) Transcription factor RELA +
59/78*
14/15**
3.21E-06^f
4.88E-04^g
1 Yes^[77]123–[78]132
STAT2 Transcription factor STAT2 +
23/24*
8/8**
1.49E-06^f
3.91E-03^g
1 Yes^[79]117,[80]133–[81]137
IRF3 Transcription factor IRF3 +
19/22*
9/9**
4.28E-04^f
1.95E-03^g
1 Yes^[82]138–[83]148
[84]Open in a new tab
^aPredicted activity of the key hub by causal reasoning is denoted by –
if the hub is inhibited, and denoted by + if the hub is activated.
^bCorrect/total network predictions show the number of genes in the
dataset predicted correctly over the total number of genes in the
causal reasoning network.
^cCalculation distance from the upstream regulatory key hubs and
downstream genes.
For example, distances of 2 and 3 identify key hubs that are distant
key hubs, while a distance of 1 identify closest one-step away
transcriptional factors.
^dp-value of the predicted protein activity calculated using the
polynomial test.
^eEvidence for an existing experimental link between the transcription
factor and menstrual irregularity.
^fp-value of the predicted protein activity calculated using the
polynomial test for GS1.
^gp-value of the predicted protein activity calculated using the
polynomial test for GS2.
*For GS1.
**For GS2.
Two of the five higher confidence transcription factors (IRF1 and IRF3)
belong to the interferon regulatory factors, and other two
transcription factors (STAT1 and STAT2) belong to the signal
transducers and transcription activators which mediate cellular
responses to interferons. The fifth transcription factor, RelA, belongs
to the Rel homology domain/immunoglobulin-like fold and is a regulator
of NF-kB activity. As a validation of these findings, we reported
examples of supporting evidence in the biomedical literature linking
these five predicted higher-confidence transcription factors to
menstrual cycle irregularities (Table [85]1). In fact, supporting
studies cited in Table [86]1 brought our attention to significant
interactions between these transcription factors and prolactin/PRL
gene.
Among all predicted transcription factors for GS1 and GS2, IRF1 had the
smallest activity prediction p-values values (4.237E-16 for GS1, and
7.63E-06 for GS2). Therefore, IRF’s causal reasoning network shown in
Fig. [87]2a, b. This network serves as an example of the causal
reasoning networks we relied in this work. Additional networks for
STAT1, RELA, IRF3 and STAT 2 are provided in Supplementary Material
(Supplementary Fig. [88]1–[89]4).
Fig. 2. Causal reasoning networks.
[90]Fig. 2
[91]Open in a new tab
a Causal reasoning network of highest confidence transcription factor
IRF1 using DEGs in GS1. b Causal reasoning network of highest
confidence transcription factor IRF1 using DEGs in GS2. Gene expression
changes are shown in green and red sectors around each molecule.
Increased expression value corresponds to the green sector which size
increases clockwise around the molecule icon. Decreased expression
value corresponds to the red sector which size increases
counterclockwise. Supportive data panel contains over and
under-expressed genes from the experimental data set which support a
hypothesis that IRF1 is in a predicted predominant “active” state.
Conflicting data panel contains over and under-expressed genes from the
experimental data set which are discordant with the hypothesis that
IRF1 is in predicted predominant “active” state.
Identifying important RNA molecules as upstream regulators
Many RNA molecules, including microRNAs and long non-coding RNAs, haven
been predicted as upstream molecular regulators that worth further
analysis by experimental scientists. DEGs in GS1, with log[2]FC ≥ +2.62
or ≤ −2.62 and FDR ≤ 0.05, led to the prioritization of 182 RNA
molecules including miR-502-5p, miR-345-5p, miR-548x-3p, miR-548x-3p,
miR-935, miR-383-5p, miR-450a-5p, miR-450a-5p, miR-4674 and miR-3941
which topped the list based on their activity prediction p-values.
However, DEGs in GS2 led to the prioritization of 28 RNA molecules
including LINC02605, miR-221-5p, NBAT1, RMRP, miR-378b, MIR31HG,
HOXB-AS1, miR-514b-5p, LINC00277 and CASC9 as top hits. Additionally,
we identified 18 overlapping RNA molecules between the 182 (from the
first DEGs in GS1) and 28 (from GS2) that are considered high
confidence RNA upstream regulators. The overlapping RNA molecules were
in order of their activity prediction p-values, from smaller to larger
values, were: LINC02605, miR-221-5p, NBAT1, RMRP, miR-378b, MIR31HG,
HOXB-AS1, miR-514b-5p, LINC00277, CASC9, ZFPM2-AS1, miR-219-1-3p,
miR-3941, LINC02605, miR-210-5p, miR-3134, MIR31HG, LOC101929517,
LINC-PINT, SBF2-AS1.
Filtering upstream regulators and downstream DEGs by biomarker uses
Menstrual irregularity biomarkers were extracted from the Cortellis
Drug Discovery Intelligence (CDDI) database^[92]46 using the following
search terms: biomarker type = “gene” or “protein”; condition =
“menstrual cycle”, “menstrual abnormalities” or “premenstrual
syndrome”; biomarker role = “diagnosis”, “disease profiling”,
“prognosis” or “prognosis–risk stratification”. Finally, we retrieved
213 biomarkers in total. We also extracted 177 biomarkers for
prolactin.
Next, we checked overlaps between the 213-biomarker set and 3 gene sets
of interest: 1) DEGs with log[2]FC ≥ +2.00 or ≤ −2.00, 2) causal hubs
predicted using DEGs with log[2]FC ≥ +2.00 or ≤ −2.00, and 3) causal
hubs predicted using DEGs with log[2]FC ≥ +5 or ≤ −5. Our results
indicated that TNF was the only common gene between our biomarkers list
and all other three gene lists. Nine genes (ARHGDIB, LIF, FGF2, MIF,
IL6R, IL6ST, RHOU, BIRC3 and RPS3) were common between the biomarkers
list and two gene sets. Finally, there were 33 common genes between
menstrual cycle biomarkers and vaccine-induced DEGs and/or predicted
causal proteins (Table [93]2), and 35 gene overlaps between prolactin
biomarkers and vaccine-induced DEGs and/or predicted causal proteins
(Table [94]3). All gene overlaps are shown in Fig. [95]3a, b. We could
look at overlaps with more gene lists as a method of filtering higher
confidence “topological” genes that may be driving the menstrual
irregularities in response to COVID-19 vaccines (i.e., increased
confidence due to biological relevance).
Table 2.
List of menstrual irregularity biomarkers which overlap with BNT162b2
vaccine-induced DEGs and/or predicted causal proteins.
# Gene Confidence Score Activity Activity Prediction p-value GS1
Activity Prediction p-value GS2 DEG
1 TNF 4 + 3.80E-11 2.96E-03 Yes
2 IL6R 3 + 3.08E-12 1.37E-04 No
3 IL6ST 3 + 3.08E-12 9.00E-05 No
4 LIF 3 + 5.24E-12 2.59E-04 No
5 BIRC3 3 + 8.89E-12 7.63E-06 No
6 FGF2 3 + 3.00E-11 1.18E-03 No
7 ARHGDIB 3 - 7.26E-11 7.63E-06 No
8 RPS3 3 + 1.25E-04 6.27E-03 No
9 RHOU 3 + 4.76E-04 3.60E-03 No
10 MIF 3 + 9.94E-03 3.60E-03 No
11 STAT4 2 + 8.48E-08 NA No
12 TEK 2 - 4.92E-06 NA No
13 CXCR4 2 + 9.06E-05 NA No
14 GAK 2 + 2.67E-04 NA No
15 ACTN1 2 + 9.62E-04 NA No
16 PGR 2 - 5.31E-03 NA No
17 MFN2 2 - 7.00E-03 NA No
18 EZH2 2 + 7.06E-03 NA No
19 AXL 2 + NA 3.69E-03 No
20 IGFBP2 2 + NA 3.81E-03 No
21 NUB1 1 NA NA NA Yes
22 ICAM1 1 NA NA NA Yes
23 PSME2 1 NA NA NA Yes
24 ADM 1 NA NA NA Yes
25 IL1B 1 NA NA NA Yes
26 HIF1A 1 NA NA NA Yes
27 GDI2 1 NA NA NA Yes
28 PHF19 1 NA NA NA Yes
29 CD1C 1 NA NA NA Yes
30 CTSW 1 NA NA NA Yes
31 KISS1R 1 NA NA NA Yes
32 DLK2 1 NA NA NA Yes
33 CCL5 1 NA NA NA Yes
[96]Open in a new tab
All molecules were ranked based on their activity prediction p-values
as well their overlap confidence score. An overlap confidence score of
3 indicates that a specific gene/protein is overlapping between the
biomarker set and 3 other gene sets, while a score of 1 indicates that
gene/protein is overlapping between the biomarker list and one other
gene set.
Table 3.
List of prolactin biomarkers which overlap with BNT162b2
vaccine-induced DEGs and/or predicted causal proteins.
# Gene Overlap confidence score Activity Activity prediction p-value
GS1 Activity prediction p-value GS2 DEG
1 TNF 3 + 3.80E-11 2.96E-03 Yes
2 MAPK14 3 + 5.46E-07 4.39E-04 Yes
3 PPARA 3 - 9.69E-06 2.97E-05 Yes
4 NFE2L2 3 - 2.32E-03 4.68E-03 Yes
5 TXNIP 3 - 3.81E-03 1.29E-03 Yes
6 IL6R 2 + 3.08E-12 1.37E-04 No
7 GNAS 2 + 9.79E-12 1.79E-05 No
8 GNAI2 2 + 1.53E-10 6.06E-05 No
9 DRD1 2 + 3.56E-09 4.88E-04 No
10 TERT 2 + 9.22E-07 5.48E-06 No
11 S100A6 2 - 1.36E-05 7.83E-05 No
12 FLT1 2 + 2.01E-05 1.63E-04 No
13 CDH1 2 - 3.28E-05 3.64E-04 No
14 MIR516A2 2 + 9.57E-04 8.45E-03 No
15 MIR516A1 2 - 9.57E-04 8.45E-03 No
16 MIF 2 + 9.94E-03 3.60E-03 No
17 VEGFA 1 + 9.58E-10 NA No
18 E2F1 1 + 1.48E-04 NA No
19 MIR576 1 - 4.01E-04 NA No
20 FTO 1 - 3.30E-03 NA No
21 AHSG 1 + 4.74E-03 NA No
22 MIR488 1 - NA 9.72E-06 No
23 CREB1 1 + NA 2.59E-04 No
24 AKT1 1 + NA 2.17E-03 No
25 NOS2 1 + NA 5.91E-03 No
26 TPT1 1 + NA 6.27E-03 No
27 CASP3 1 - NA 6.47E-03 No
28 CCL2 1 NA NA NA Yes
29 CD274 1 NA NA NA Yes
30 IFI44 1 NA NA NA Yes
31 MX2 1 NA NA NA Yes
32 IL1B 1 NA NA NA Yes
33 BMPR2 1 NA NA NA Yes
34 PARP1 1 NA NA NA Yes
35 SERPINF2 1 NA NA NA Yes
[97]Open in a new tab
All molecules were ranked based on their activity prediction p-values
as well their overlap confidence score. An overlap confidence score of
3 indicates that a specific gene/protein is overlapping between the
biomarker set and 3 other gene sets, while a score of 1 indicates that
gene/protein is overlapping between the biomarker list and one other
gene set.
Fig. 3. Overlapping DEGs, causal upstream hubs and biomarkers.
Fig. 3
[98]Open in a new tab
a Venn diagram showing overlaps between DEGs, predicted causal upstream
regulatory hubs using DEGs in GS1 and GS2, and menstrual irregularity
biomarkers. b Venn diagram showing overlaps between DEGs, predicted
causal upstream regulatory hubs using DEGs in GS1 and GS2, and
prolactin signaling biomarkers.
Pathway enrichment and interconnectivity analysis
We used pathway enrichment analysis to assess whether the identified
causal regulators work collectively to affect certain biological
pathways. Pathway enrichment results using five higher confidence
upstream causal regulators as a query gene list (IRF1, STAT1, RELA,
IRF3 and STAT 2), highlighted the prolactin signaling pathway as one of
the significantly enriched pathways with an enrichment p-value of
1.60E-05 (Fig. [99]4a). Next, we included PRL in the query list (PRL,
IRF1, STAT1, RELA, IRF3 and STAT 2) for pathway enrichment which led to
the prioritization of prolactin singling pathway as the top enriched
pathway with an enrichment p-value of 8.92E-07 (Fig. [100]4b). Similar
analyses were performed on the 33 filtered biomarker genes/proteins
(Fig. [101]4c), and biomarker proteins in addition to PRL, IRF1, STAT1,
RELA, IRF3 and STAT 2 (Fig. [102]4d).
Fig. 4. Interconnectivity between prioritized high confidence transcription
factors.
[103]Fig. 4
[104]Open in a new tab
a Direct interactions network of five higher confidence causal
transcription factors. b Direct interactions network of five higher
confidence causal transcription factors in addition to prolactin (PRL).
c Direct interactions network of 33 causal upstream regulators that are
known biomarkers for menstrual disturbances. d Direct interactions
network of 33 causal upstream regulators that are known biomarkers for
menstrual disturbances, in addition to prioritized 5 topological genes
and PRL. Thick edges correspond to confidence score
[MATH: ≥ :MATH]
0.70 (i.e., high confidence score), while the thin edge corresponds to
a confidence level
[MATH: ≤ :MATH]
0.50 (i.e., low confidence score). Nodes are color-coded using a split
pie chart coloring scheme indicating pathway/gene set contribution to
each node from the top 5 most enriched pathways/gene lists. FDR values
represent he significance of the predicted pathway. Generated based on
STRING data on 27 September 2022.
It is well known that changes in prolactin signaling can result in
menstrual cycle irregularities^[105]47–[106]51. Searching PubMed using
terms “prolactin” AND “menstrual” returned 1950 results, while search
terms “prolactin” AND “menstrual cycle” returned 1270 results, and
search terms “prolactin” AND “menstrual irregularity” returned 32
results. This is a clear indication that the studied mRNA COVID-19
vaccine BNT162b2 has the potential to cause menstrual irregularities by
inducing perturbations in the genes and/or proteins involved in
prolactin signaling pathways (i.e., through affecting the activity of
key transcription factors involved in this pathway). Prolactin
signaling pathway affects a wide range of physiological processes
ranging from reproduction and lactation to growth and development, from
endocrinology and metabolism to brain and behavior, as well as immune
regulation (Fig. [107]5a, b).
Fig. 5. Prolactin signaling pathway.
[108]Fig. 5
[109]Open in a new tab
a Prolactin signaling pathway map. A node (or object) on the map could
be a gene, protein or chemical compound. Query genes from experimental
data which intersect with pathway protein or chemical compound. Query
genes from experimental data which intersect with pathway objects are
denoted by thermometers. Thermometer 1 represents causal transcription
factors. Thermometer 2 represents DEGs in response to treatment with
vaccine, applying thresholds of log[2]FC ≥ 2.00 or ≤ −2.00, and
FDR ≤ 0.05. b Biological processes involved in prolactin signaling
pathway. The % refers to the percentage of network objects in the
pathway map. The p-value is the process prediction p-value.
Prolactin is a polypeptide hormone encoded by the PRL gene and secreted
by the anterior pituitary gland. It is known as a growth regulator for
many tissues, including cells of the immune system. It functions as a
cytokine with immunomodulatory activities. It may also play a role in
cell survival by suppressing apoptosis, and it is essential for
lactation. Chemically, prolactin’s structure is similar to those of the
growth hormone and the placental lactogen hormone, which form together
what is known as the “prolactin/growth hormone/placental lactogen”
family, and they all originated from one ancestral gene.
Signaling pathway impact analysis (SPIA)
SPIA was performed on the combined list consisting of the DEGs and the
predicted causal hubs. Enrichment analyses performed on the combined
list of DEGs, and key hubs have previously been shown to highlight more
biologically relevant results^[110]52. The top five enriched pathway
maps were: 1) immune response interferon-alpha/beta (IFN-alpha/beta)
signaling via JAK/STAT, 2) regulation of antiviral response by
SARS-CoV-2, 3) immune response antiviral actions of interferons, 4)
immune response induction of apoptosis and inhibition of proliferation
mediated by interferon-gamma (IFN-gamma), and 5) immune response
IFN-gamma in macrophages activation with impact p-values of 4.31E-25,
1.24E-11, 1.04E-10, 5.61E-10, and 1.59E-9 respectively. All predicted
pathways are highlighting the role of interferons. Hence, SPIA results
revealed a prominent role of interferon signaling in the signaling
pathways impacted by BNT162b2 vaccines.
Supporting evidence
We mined the VAERS database, PubMed using Abstract Sifter, and the
Connectivity Map to gather supporting evidence for the prioritized
hypotheses regarding vaccine-induced menstrual irregularities and
prolactin or HPO mimicking effects.
VAERS database
We searched the VAERS database for vaccine adverse events that are
relevant to menstrual irregularities. First, we extracted all COVID-19
vaccine adverse events and ranked them according to their proportions
of the overall reported adverse events for COVID-19 vaccines. VAERS
annotates menstrual irregularity under ‘menstruation irregular’, in
addition to using more specific “symptom” terms, including heavy
menstrual bleeding, dysmenorrhea, intermenstrual bleeding, amenorrhea,
postmenopausal hemorrhage, premenstrual pain, menstrual discomfort,
menstruation normal, premenstrual dysphoric disorder, menstrual cycle
management, premenstrual headache, retrograde menstruation,
premenstrual syndrome and menstrual disorder.
Our analysis indicated that none of the menstruation disturbances
listed above were among the most frequently reported adverse events for
COVID-19 vaccines that we detailed before^[111]12. An updated adverse
event report for COVID-19 vaccines is provided in Supplementary Table
[112]5. However, some cross-sectional studies reported high frequencies
of these side effects. There could be many explanations for this
discrepancy between adverse event databases and cross-sectional
studies. Women tend to neglect seeking medical attention for what they
perceive as a mild non-threatening short-term menstruation
irregularity^[113]53–[114]55. In fact, women participating in
cross-sectional studies are unlikely to report changes to periods
unless specifically asked^[115]16.
In fact, mining VAERS data for menstruation irregularities resulted in
35,386 adverse events that were not restricted to COVID-19 vaccines
(Supplementary Table [116]6). The top five vaccines that had the
highest share in these events were COVID-19 vaccine (26,714 events
comprising 85.82%), human papillomavirus recombinant vaccine (1198
events comprising 3.85%), hepatitis B vaccine (1013 comprising 3.25%),
trivalent influenza virus vaccine (581 events comprising 1.87%) and the
zoster vaccine (566 events comprising 1.82%). It is noteworthy that all
these vaccines are given later rather than the first few years of life,
permitting adverse event reporting by menstruating women.
PubMed Abstract Sifter
To increase the confidence in the prioritized causal hits, we examined
the relationship(s) between the prioritized genes and menstruation
irregularities in more depth and complexity using the PubMed Abstract
Sifter^[117]56. On the Landscape sheet we built queries that revealed
the number of articles satisfying a variety of queries related to
menstrual cycle, abnormal menstruation, vaccines, and prioritized
upstream causal regulators. These results (Supplementary Table [118]7)
revealed a sizable publication record for the relationship between
menstruation and genes affected by COVID-19 vaccines (DEGs and/or
causal hubs). The results of two queries consisting of the higher
confidence list of prioritized causal transcription factors, and
biomarkers causal hubs and DEGs are shown in Fig. [119]6a, b.
Fig. 6. Screenshots from PubMed Abstract Sifter.
[120]Fig. 6
[121]Open in a new tab
a Landscape sheet of the PubMed Abstract Sifter showing relationships,
in the form of article counts, between biological concepts highlighted
in this study. The first column “id” lists the gene symbols of
prioritized top five causal transcription factors. b Landscape sheet of
the PubMed Abstract Sifter showing relationships, in the form of
article counts, between causal biological concepts highlighted in this
study. The first column “id” lists the gene symbols of prioritized
causal genes and vaccine-induced DEGs that are known as also biomarkers
for menstrual cycle according to the CDDI database^[122]46.
The Connectivity Map (CMap)
The Connectivity Map analysis suggests that drugs capable of inducing
transcriptomics effects opposite to those induced by mRNA vaccines
could reverse vaccine side effects^[123]57,[124]58. To identify
small-molecule drugs that could prevent or reverse vaccine’s side
effects, we ranked all DEGs in response to vaccination by BNT162b2
according to their expression levels using log[2]FC values, to query
the Connectivity Map database^[125]59. The CMap query gene signature
consisted of the 50 most upregulated genes and the 50 most
downregulated genes in response to vaccination with BNT162b2. Compound
hits that produced opposite transcriptional signatures to the mRNA
vaccine BNT162b2 are listed in Table [126]4. These compounds can
reverse the transcriptomic signature of the vaccine, which will prevent
or reduce side effects. In this study, we wanted to increase the
confidence in the computational hypotheses derived from the enrichment
and network analyses described earlier.
Table 4.
Small-molecule drugs and chemical compounds that regulate gene
expression in an opposite manner to BNT162b2.
Compound Score Description Links to hypothalamic-pituitary-ovarian
function
Droxinostat −94.51 HDAC inhibitor ^[127]149
Metyrapone −93.62 Cytochrome P450 inhibitor ^[128]150,[129]151
Perospirone −92.07 Dopamine receptor antagonist ^[130]152
Nabumetone −88.8 Cyclooxygenase inhibitor ^[131]153
Salbutamol −87.35 Adrenergic receptor agonist ^[132]154,[133]155
VU-0415374-1 −86.32 Glutamate receptor modulator ^[134]156
Bromfenac −83.64 Cyclooxygenase inhibitor ^[135]153,[136]157
PF-3845 −82.22 FAAH inhibitor ^[137]158
Hexylresorcinol −80.48 Local anesthetic ^[138]159
PPT −80.07 Estrogen receptor agonist ^[139]160
[140]Open in a new tab
Discussion
This study describes the first attempt to provide a mechanistic insight
for vaccine-induced menstrual cycle irregularities. Our approach
combined the analysis of vaccine gene expression profiles with upstream
predictions of causal regulatory proteins and RNAs, and downstream
analysis of enriched biological pathways to provide a causal
mechanistic insight for vaccine-induced menstrual irregularities.
Our analysis led to the prioritization of topologically significant
genes, such as transcription factors and important enzymes (i.e.,
kinases) that were largely missed in the gene expression profiling
experiment, and therefore were not among the prioritized DEGs. We used
the ‘Causal Reasoning’ methodology to identify candidate proteins
(i.e., hypotheses) in the network that can be reached through a
pre-defined distance (i.e., maximum shortest path length) from the
DEGs. Thus, this analysis was crucial for the reconstruction networks
responsible for vaccine-induced menstrual irregularities. The top five
transcription factors, listed from highest confidence to lower
confidence based on their prediction p-values, were: IRF1, STAT1, RELA,
IRF3 and STAT 2 (Table [141]1). All were predicted to be activated, in
response to BNT162b2 vaccination, based on the directionality of
differential gene expression in GS1 and GS2.
IRF1 was ranked first as the highest confidence predicted activated
transcription factor. To assess whether changes in IRF1 activation can
affect menstrual cycle, we checked whether IRF1 is biomarker for
“menstrual cycle irregularity” but we didn’t find evidence to support
that. Next, we reviewed the biomedical literature to search for
possible links between IRF1 and the menstrual cycle. We used the
PubMed’s advanced search using query terms “IRF1” and “menstrual cycle”
and found evidence that IRF1 is upregulated by prolactin during the
secretory phase of the menstrual cycle^[142]59,[143]60. Additionally,
evidence pointed that IRF1 upregulation in the endometrium was linked
to prolactin and is localized predominantly to the granular epithelial
cells^[144]59. Network reconstruction using PLR in addition to seed
nodes IRF1, STAT1, RELA, IRF3 and STAT led to the direct interactions
network in Fig. [145]3a. Downstream enrichment analysis in biological
pathways, highlighted the prolactin signaling pathway as the most
significantly enriched pathway with the six network seeds mentioned
above.
Thus, upstream causal reasoning followed by downstream pathways
analysis highlighted a putative role for prolactin signaling in
modulating post-COVID-19-vaccine adverse events on the menstrual cycle.
Prolactin is a multi-functional molecule; it is a transcription factor
hormone, secreted from the pituitary glands, and it regulates diverse
biological functions including female menstruation^[146]61–[147]69. For
example, high prolactin levels can interfere with the production of sex
hormones including estrogen and progesterone which can further impact
menstruation regulation^[148]61–[149]69. In fact, women who experience
menstrual cycle irregularities often have higher prolactin levels than
others, a condition known as hyperprolactinemia^[150]47.
Hyperprolactinemia, is the most prevalent endocrine dysfunction of the
hypothalamic-pituitary axis in young females, accompanied with
ovulatory disorder and leading to menstrual
irregularities^[151]70,[152]71. High levels of prolactin in the body
prevent the release of (luteinizing hormone) LH and
follicle-stimulating hormone (FSH), leading to ovulation
disturbances^[153]62,[154]65,[155]69. Symptoms of hyperprolactinemia
include long or irregular cycles, anovulation, amenorrhea,
oligomenorrhea, polycystic ovarian syndrome or
amenorrhea^[156]72–[157]78. In fact, hyperprolactinemia can be caused
by some drugs, stress, and some conditions like prolactinoma
(noncancerous tumor of the pituitary gland)^[158]74,[159]79–[160]81.
All these factors were found to cause inconsistencies in menstrual
cyclicity^[161]51,[162]82–[163]84.
Although we perceive menstrual changes as adverse events,
prolactin-mimicking effects of vaccine are not necessarily a negative
consequence of vaccines. Recently, prolactin has been suggested as a
promising immunomodulator for the treatment of COVID-19^[164]85.
However, we caution that prolactin mimicking effects may worsen
auto-immune disease symptoms in patients suffering from systemic lupus
erythematosus (SLE), multiple sclerosis, rheumatoid arthritis,
psoriatic arthritis, and AIDS. Caution should be also practiced in
patients undergoing organ transplantation. Elevated PRL levels have
been reported in the previous conditions^[165]86,[166]87. Furthermore,
a recent study showed that prolactin hormones in addition to FSH and LH
of healthy vaccinated males were higher than non-vaccinated males or
COVID-19 male patients, indicating that changes in prolactin signaling
are not limited to females^[167]88. Prolactin levels were
27.86 ± 4.35 ng/L in vaccinated males, 5.35 ± 1.59 ng/L in
non-vaccinated males, and 16.65 ± 6.15 ng/L in COVID-19 male
patients^[168]88.
To identify biomolecules that are implicated in menstrual changes, or
the pathological processes that underlie the observed vaccine-induced
menstrual symptoms, we filtered all predicted causal molecules and DEGs
based on their overlaps with “menstruation irregularity”/“menstruation
abnormality” diagnostic and prognostic biomarkers found in CDDI (Fig.
[169]2a). We had four gene lists: 1) all DEGs in GS1 and GS2, 2) causal
hubs for DEGs predicted for GS1, 3) causal hubs predicted for GS2, and
4) known diagnostic biomarkers for menstrual irregularity. TNF was
identified as a high-confidence hit, i.e., a causal protein leading to
the observed changes in gene expression and the predicted prolactin
mimicking effects of the vaccine. TNF was an overlapping gene among
four gene lists: 1) a DEG (log[2]FC = 3.07), 2) a causal key hub
considering DEGs in both GS1, 3) a causal hub considering DEGs in GS2,
and 4) a diagnostic biomarker for menstrual irregularity. Furthermore,
our causal reasoning results predicted TNF-alpha activation in response
to vaccination with BNT162b2. It should be noted that all causal
predictions (Supplementary Tables [170]1 & [171]2) are based on
experimental gene expression data.
Finally, SPIA results combined the enrichment results of DEGs with the
actual amount of perturbation which highlighted the role of interferons
on the signaling pathways influenced by BNT162b2. In fact, mRNA and
vector-based COVID-19 vaccines result in the formation of neutralizing
antibodies and activation of immune cells via the release of
pro-inflammatory markers like cytokines and interferons^[172]89. There
is evidence indicating that the treatment of multiple sclerosis with
beta interferons causes menstrual irregularities associated with
increased levels of luteinizing hormone (LH) and/or
hyperprolactinemia^[173]90. Furthermore, the upregulation of
interferon-gamma perturbs calcium signaling pathways which can in turn
impact hormonal balance^[174]12.
But what is the relationship between prolactin signaling, TNF-alpha
activation and interferons? In fact, TNF-alpha activates the human
prolactin gene promoter via NF-κB signaling^[175]91. TNF-alpha
activation also stimulates the hypothalamic-pituitary-adrenal axis
while suppressing the hypothalamic-pituitary-thyroid and gonadal axes,
and growth hormone release^[176]92. Menstrual bleeding (menses) is
known to be regulated by hypothalamic and pituitary hormones, and even
the slightest changes in hormone levels, e.g., in response to
medication or stress, can result in menstrual cycle
abnormalities^[177]93. There is evidence that TNF-alpha and interleukin
1 beta (IL-1B), both are upregulated DEGs in this analysis, exert
significant inhibitory effects on the GnRH-LH system in rats^[178]94,
which may be the case in humans too. Moreover, the occurrence of
reproductive disorders in poultry is highly correlated with the HPO
axis and neuro–endocrine–immune network molecules, such as TNF-alpha
and interferon-gamma (IFN-γ, IFNG)^[179]95. Thus, integrating
enrichment and causal reasoning results with SPIA findings uncovered
causal relationships between BNT162b2-induced menstrual changes and all
the following pathways: prolactin signaling pathways, TNF-alpha
activation, interferons the
hypothalamic-pituitary-gonadal/ovarian/testicular axis. These results
agree with previous studies suggesting that stabilizing the
hypothalamic-pituitary-ovarian (HPO) axis with combined hormonal
contraception reduces the chance of experiencing vaccine-associated
menstrual changes^[180]38,[181]96.
Different lines of supporting evidence increased the confidence in the
derived causal hypothesis implicating menstrual changes with prolactin
signaling, TNF-alpha and the HPO axis. First, VAERS data showed that
post COVID-19 menstrual changes occurred in response all COVID-19
vaccines included in the databases including mRNA and vector-based
vaccines and were not tied to the vaccine platform. The menstrual
changes reported in VAERS included wide range of symptoms and were not
limited to the length of menstrual cycle or menses period. Secondly,
PubMed Abstract Sifter results highlighted 299,927 articles linking DEG
TNF to any menstrual symptoms and 141 articles linking TNF to abnormal
menses. Other high-confidence causal DEGs were progesterone receptor
(PGR) with 45,560 and 796 articles linking it to any menstrual symptoms
or abnormal menses subsequently, IL-1B with 56,099 and 43 articles
linking it to any menstrual symptoms or abnormal menses subsequently.
Finally, chemogenomics evidence from the CMap highlighted significant
links to the HPO axis per results shown is Table [182]4.
It should be noted that the transcriptomics perturbations in response
to treatment with BNT162b2 diminished on day 28 after receiving the
second vaccine dose of BNT162b2. This suggests that vaccine effects on
gonadal hormones, for females and males, and the predicted
prolactin-mimicking effects, TNF-alpha activation, and HPO signaling
changes, were temporary. However, we cannot rule out long-term effects
without clinical studies comparing vaccinated and non-vaccinated
individuals. Moreover, because our bioinformatics analysis relied on
BNT162b2 transcriptomics data, we emphasize that these findings
primarily apply to the BNT162b2 vaccines. However, data mined from
VAERS, and the biomedical literature indicated that vaccine-induced
menstrual cycle changes were reported for other COVID-19 vaccines
(e.g., mRNA-1273, and Janssen’s) and non-COVID-19 vaccines (e.g., HPV
and typhoid).
In conclusion, our integrative computational network biology approach
revealed that BNT162b2 can induce transcriptomics changes which may
induce menstrual cycle changes by several mechanisms including
prolactin-mimicking effects resulting from changes in interferon
signaling and associated hormonal imbalance particularly in the HPO
axis. This remains a high-confidence biological hypothesis supported by
different lines of computational evidence derived from transcriptomics
studies, causal reasoning analysis, downstream pathway enrichment
results, and additional supporting evidence from vaccine adverse event
databases (e.g., VAERS) and the biomedical literature. Further
experimental validation is warranted to assess whether post-vaccine
prolactin-mimicking effects are due to increased levels of prolactin or
due to other networking biology events mimicking prolactinemia. These
effects may not be restricted to COVID-19 vaccines and should be
assessed for other vaccines as well.
This study sheds the light on post-vaccine menstrual irregularity by
revealing short-term post-COVID-19 vaccine prolactin mimicking effects
resulting from the transcriptomics irregularities induced by COVID-19
vaccines. Most women associate menstruation irregularities with
infertility which is one of the leading causes of vaccine hesitancy
among females^[183]97. By providing a mechanistic insight into
post-vaccine menstrual irregularities, this study is promised to
correct misinformation about the relationship between vaccine-induced
period irregularities and infertility. Thus, it is expected to decrease
vaccine hesitancy.
This study establishes a causal relationship between COVID-19 vaccine
and menstruation regulation by highlighting perturbed gene expression
or dysregulated transcription of known diagnostic or prognostic
biomarkers for menstruation and menstruation irregularities.
Additionally, top scoring key hubs provide valuable hypotheses
explaining gene expression and can be explored further in laboratory
tests.
This study explores the causal links between COVID-19 vaccines and
menstruation regulation based on an integrative bioinformatics approach
that analyzed vaccine-induced transcriptomics irregularities.
Integrating COVID-19 vaccine transcriptomics data with menstruation
biomarkers, reinforced the selection of biologically relevant
hypotheses from an overwhelming number of statistically significant
hypotheses by increasing the confidence in computational hypotheses
predicted by several methods. The fact that our computational
hypotheses were supported by multiple lines of evidence is considered a
major strength for this study. In fact, our integrative informatics
workflow has several advantages over relying solely on conventional
enrichment analyses for identifying the biological mechanisms that
underlie vaccine side effects. Our approach integrates hypotheses
derived independently from pathway and network enrichments, causal
reasoning, SPIA, and the CMap to prioritize high confidence
computational hypotheses predicted independently by various
computational approaches and using different data types. The CMap, for
example, is considered a unique chemogenomics database capable of
connecting genes, drugs, and diseases based on genes expression
similarities between polypharmacologic drugs and studied vaccines. This
permits the prediction of vaccine side effects as well as underling
causal mechanisms based on gene expression similarities with
well-studied drugs. Finally, mining VAERS and PubMed for adverse event
reports and vaccine-relevant data, serves as a validation step for the
computationally-derived hypotheses. Thus, computational hypotheses
prioritized using our integrative informatics approach are inherently
high-confidence hypotheses with potentially improved clinical outcome.
Conversely, the applied methodologies or public databases have a few
limitations that should be pointed out. First, reports from VAERS may
not be conclusive or sufficient to establish causal relationships
adverse events and specific vaccines. Due to the voluntary nature of
VAERS reporting system, the information provided about an adverse event
can be imperfect, imprecise, coincidental, or unconfirmed, limiting the
scientific use of such reports^[184]11,[185]12. Secondly,
bioinformatics techniques relying on gene expression, pathway
over-representation and network biology have some limitations and
biases that we reviewed previously elsewhere^[186]43. Herein, the main
limitation for the generalizability of the bioinformatics results to
other COVID-19 vaccines, was the reliance on transcriptional data for
the mRNA COVID-19 vaccine BNT162b2, which was the only publicly
available COVID-19 transcriptomics data in humans at the time of
conducting this research. As a result, our bioinformatics results apply
directly to BNT162b2 or and may be extended to other COVID-19 mRNA
vaccines (e.g., mRNA-1273) since COVID-19 mRNA vaccines share common
features of the nature, strength, and timing of the immune responses as
well as similar vaccine
compositions^[187]7,[188]8,[189]12,[190]18,[191]44. The dosing regimens
of vaccines may affect the results as well^[192]30,[193]89. Our
integrative workflow can be used to assess the safety of other vaccines
using their transcriptional signatures in vaccinated individuals.
Methods
Systems biology informatics workflow
We have developed a network biology workflow to identify causal links
between COVID-19 Vaccines and menstruation irregularities. This
workflow (Fig. [194]1) incorporates three major components: (1) a
module for mining and prioritizing gene signatures representative of a
condition or a biological state; (2) a causal reasoning network module
to identify upstream regulators of gene expression and (3) a pathway
enrichment module to understand the biological processes regulated by
DEGs and predicted causal regulators of gene expression. The resulting
hypotheses are then evaluated based on existing evidence in vaccine
reporting system databases and the biomedical literature.
Vaccine-induced differential gene expression
We searched the gene expression omnibus (GEO)^[195]98–[196]101 for
transcriptional studies performed in response to treatment COVID-19
vaccines and we were able to identify one whole transcriptomics RNA-seq
dataset ([197]GSE169159) for COVID-19 vaccines in response to treatment
with BNT162b2 at different time points. Our transcriptomics data
analysis of [198]GSE169159 raw data indicated that gene expression
alterations from baseline were more prominent on day 22, which is the
day after receiving the vaccine second dose. None of the genes analyzed
at other time points (e.g., day 7, day 21, day, day 22.23, day 28)
passed the applied thresholds for the prioritization of DEGs in this
study (i.e., log[2] fold change (log[2]FC) of +2 or –2, and false
discovery rate (FDR)
[MATH: ≤ :MATH]
0.05. Therefore, we relied on differential gene expressions on day 22
for all our bioinformatics analyses.
Gene expression profiles on day 22 were used to generate two query gene
signatures to study the systems biology effects of BNT162b2: GS1 and
GS2. GS1 consisted of all differentially expressed genes
(DEGs)^[199]102–[200]104 that passed our prioritization criteria for
DEGs: 1) log[2] FC
[MATH: ≥ :MATH]
2.00 for differentially upregulated genes, and
[MATH: ≤ :MATH]
−2.00 for differentially downregulated genes; 2) FDRs
[MATH: ≤ :MATH]
0.05. GS2 consisted of all differentially expressed genes (DEGs) that
passed our prioritization criteria for DEGs: 1) log[2]FC
[MATH: ≥ :MATH]
5.00 for differentially upregulated genes, and
[MATH: ≤ :MATH]
−5.00 for differentially downregulated genes; 2) FDRs
[MATH: ≤ :MATH]
0.05. The DEGs used to derive GS1 and GS2 are provided in Supplementary
Table [201]8.
Upstream transcriptomics analysis
Causal Reasoning^[202]45,[203]105 was used to identify upstream
regulators (transcription factors, RNA molecules, kinases,
phosphatases, and others proteins) that could cause/explain the
observed post-vaccine gene expression changes in transcriptomics
experiments. We relied on MetaBase^[204]106,[205]107 as an interactions
database, and the causal reasoning algorithm implemented in Clarivate’s
Key Pathway Advisor^[206]103,[207]108. This method relies on a directed
network which is annotated with activation and inhibition edges as well
as biological mechanisms (transcription regulation). The significance
of the predictions made by a particular hypothesis is assessed using a
binomial test and calculating p-values as probabilities to get k
successes in n predictions using binomial tests with p-value = 0.50
according to the following equation:
[MATH:
p-value=nkpk(1−p)n−k :MATH]
1
Here, k is the sum of correct predictions and n is the sum of correct
and incorrect predictions.
Finally, p-values are assigned in the score matrix and hypotheses above
the p-value threshold are filtered out of the score matrix.
Downstream pathway analysis
Pathway enrichment analyses were conducted in Cytoscape version
3.9.1^[208]109 and MetaCore^TM^[209]45 to interpret the consequences of
vaccine-induced differential gene expression on biological processes.
The significance of the enrichment was determined by calculating
hypergeometric p-values^[210]110. All terms from the ontology were
ranked based on their calculated p values. Ontology terms with p-values
less than the p-value threshold 0.05 are defined as statistically
significant and therefore relevant to the studied list of genes. All
terms from the ontology were ranked according to their calculated
p-values.
Signaling pathway impact analysis (SPIA)
SPIA^[211]111,[212]112 was performed to identify the impact of the DEGs
on the activity of the enriched pathway. This method aids in the
identification of the most biologically relevant pathways and candidate
causal genes. Herein, we identified perturbed pathways in response to
vaccination by performing the enrichment analysis on the union gene
list consisting of the experimentally derived DEGs in response to
vaccination with BNT162b2, and the list of key hubs (e.g., activated,
or inhibited proteins) using causal reasoning.
Vaccine adverse events database
Raw data files were downloaded in comma-separated value (CSV) files
from the CDC website^[213]106,[214]107. CDC WONDER online search tool
was used to mine VAERS data by vaccine type and symptoms^[215]108. The
COVID-19 vaccines included in the databases were: BNT162b2, mRNA-1273
’s and Janssen’s^[216]113,[217]114.
PubMed Abstract Sifter
The advanced literature retrieval tool PubMed Abstract Sifter was used
to explore relationships between the biological concepts and molecular
concepts that play roles in this research area. The steps in using the
Abstract Sifter are described in the user guide. The tool and the user
guide are available from the US EPA and downloadable from this webpage:
[218]https://comptox.epa.gov/dashboard/downloads^[219]56,[220]115–[221]
117.
The Connectivity Map (CMap)
The Connectivity Map analysis suggests that drugs capable of inducing
transcriptomics effects opposite to those induced by mRNA vaccines
could reverse vaccine side effects^[222]57,[223]58. To identify
small-molecule drugs that could prevent or reverse vaccine’s side
effects, we ranked all DEGs in response to vaccination by BNT162b2
according to their expression levels using log[2]FC values, to query
the Connectivity Map database^[224]59. In fact, our transcriptomics
data analysis of [225]GSE169159 indicated that gene expression
alterations from baseline were more prominent on day 22, which is the
day after receiving the vaccine second dose. None of the genes analyzed
at other time points passed the applied thresholds for identifying DEGs
(i.e., log2 fold change (log2FC) of +2 or –2, and false discovery rate
(FDR)
[MATH: ≤ :MATH]
0.05. Therefore, we relied on differential gene expressions on day 22
for all our bioinformatics analyses.
Reporting summary
Further information on research design is available in the [226]Nature
Research Reporting Summary linked to this article.
Supplementary information
[227]Supplemental Information^ (2.5MB, pdf)
[228]Reporting Summary^ (72.4KB, pdf)
Acknowledgements