Integration of Transcriptome and Metabolome Provides Unique Insights to Pathways Associated With Obese Breast Cancer Patients.

Information regarding transcriptome and metabolome has significantly contributed to identifying potential therapeutic targets for the management of a variety of cancers. Obesity has profound effects on both cancer cell transcriptome and metabolome that can affect the outcome of cancer therapy. The information regarding the potential effects of obesity on breast cancer (BC) transcriptome, metabolome, and its integration to identify novel pathways related to disease progression are still elusive. We assessed the whole blood transcriptome and serum metabolome, as circulating metabolites, of obese BC patients compared them with non-obese BC patients. In these patients' samples, 186 significant differentially expressed genes (DEGs) were identified, comprising 156 upregulated and 30 downregulated. The expressions of these gene were confirmed by qRT-PCR. Furthermore, 96 deregulated metabolites were identified as untargeted metabolomics in the same group of patients. These detected DEGs and deregulated metabolites enriched in many cellular pathways. Further investigation, by integration analysis between transcriptomics and metabolomics data at the pathway levels, revealed seven unique enriched pathways in obese BC patients when compared with non-obese BC patients, which may provide resistance for BC cells to dodge the circulating immune cells in the blood. In conclusion, this study provides information on the unique pathways altered at transcriptome and metabolome levels in obese BC patients that could provide an important tool for researchers and contribute further to knowledge on the molecular interaction between obesity and BC. Further studies are needed to confirm this and to elucidate the exact underlying mechanism for the effects of obesity on the BC initiation or/and progression.

Introduction

Breast cancer (BC), the most common malignant tumor type, was ranked top in incidence with high prevalence and mortality among females in Saudi Arabia as well as worldwide (1–3). BC is a molecularly heterogeneous, complex, and multifactorial disease with different biological and clinical characteristics (4). A number of BC-related etiological factors have been identified as hereditary, genetic factors, environmental, and lifestyle risk factors (5). Obesity poses a serious public health issue worldwide (6). In Saudi Arabia, the prevalence of obesity is 28.7% with a higher incidence among women (7). Obesity is one of the risk factors associated with the development of many types of cancer including BC. A number of studies, a few of them in Saudi Arabia, have reported an association between obesity and BC among postmenopausal women whereas the inverse relationship was reported among premenopausal women, however, this association remains unclear (8–13).

Obesity–BC molecular interaction could provide an important tool for researchers, as it may help to identify and discover new molecular fingerprints, as well as clarify molecular mechanisms involved in screening and develop therapeutic strategies for the management of BC. Notably, a few studies have focused on the molecular interaction between obesity and BC but this association still remains unclear. Recently, omics techniques, such as transcriptomics and metabolomics, have been widely used to improve understanding of the underlying biological mechanisms and biomarkers identification (14). Many studies utilized transcriptomics to investigate human diseases at the molecular level and identify the variation in transcriptomic profile in relation to diagnostic, treatment, or management and that may help to understand the mechanisms of disease initiation and progression (15–19). Moreover, many transcriptomic studies have been based on the analysis of the association between BC and obesity (18, 20). On the other hand, metabolomic investigations have been widely utilized in cancer metabolism and biomarkers identification to deduce the onset and progression of cancer (21, 22). Increasingly, studies now include measurements from multiple omics techniques rather than the single omics technology of a set of samples in early studies (14). There is currently very limited published research that has investigated the integration between deregulated transcriptomics and metabolomics data profiling in BC blood liquid biopsy. Therefore, in our study, we performed the integration between deregulated transcriptomics data and metabolomics profiling in BC patients with obesity to provide a better understanding of the biological status and shed new insights into potential molecular mechanisms, the interactions and biomarkers in the relationship between obesity and BC. Pathways and network connections were carried out to further explore the relationship between the selected metabolites and candidate transcripts. This could give considerable importance for the clinical management of BC patients and could provide an important tool for researchers as well as to increasing the knowledge on the molecular interaction between obesity and BC.

Materials and Methods

Study Subjects

We summarized the study workflow for identification and validation of signature RNAs and metabolites in obese BC patients in Figure 1. The study includes 69 newly diagnosed and before any treatment BC female subject who attended the Unit of Mammogram, Department of Radiography at King Abdulaziz University Hospital (KAUH), Jeddah, Saudi Arabia. The Unit of Biomedical Ethics, Research Committee, approved this study (Document number: HA-02-J-008). All BC patients signed the consent form. The patients' information was obtained through a standard questionnaire and the anthropometric data were collected using standard and well-established methods. The clinicopathological characteristics were obtained in collaboration with the Pathology Department at KAUH, Jeddah, Saudi Arabia.

Figure 1

Flowchart of transcriptomics and metabolomics analysis in obese vs. non-obese BC patients. BC, Breast cancer; BMI, Body mass index; DEGs, differentially expressed genes; GO, Gene ontology.

The WHO recommendations (23) were used to classify BC patients as obese [Body mass index (BMI) ≥ 30 kg/m2, n = 36] and non-obese, which include lean and overweight (BMI < 30 kg/m2, n = 33). From all patients' cohort, RNA sequencing (RNA-seq) was performed for only 10 non-obese with lowest BMI and 11 obese BC patients with highest BMI. Then the transcriptomics data were validated in all groups, non-obese (n = 33) and obese (n = 36) BC patients. However, from the RNA samples cohort, only six samples of each group were selected for the metabolomics. These selections were based on BMI differentiation as with the RNA-seq samples selection.

Blood and Serum Sample Collection and Storage

Whole blood samples were collected in PAXgene™ blood RNA tubes (PreAnalytiX, Switzerland) as well as in BD Vacutainer™ venous blood collection tubes: SST™ serum separation tubes (Fisher Scientific, USA), according to the manufacturer's instructions. Serum samples were aliquots after separated from the clotted blood. Blood samples in PAXgene™ blood RNA tubes, as well as aliquots serum samples, were stored at −80°C until used for transcriptome and metabolome analysis, respectively.

Transcriptome Analysis

RNA Extraction

Total RNA was isolated from whole blood using the PAXgeneTM blood RNA kit (Qiagen, UK). The concentration and purity of the extracted RNA were verified by DeNovix DS-11 Spectrophotometer (DeNovix, USA) and Agilent 2100 bioanalyzer measurements (Agilent Technologies, USA). The RNA samples were stored at −80°C until used.

RNA Library Preparation, Sequencing, and Differentially Expressed Genes Analysis

The next-generation sequencing technologies were used for performed RNA-seq experiment to discover the amount of RNA in a blood biological sample at a given moment by using the Next-Seq 500 system (Illumina, Singapore) as described elsewhere (24). Approximately 2 μg of total RNA was fragmented and end-repaired using the Illumina directional protocol. Complementary DNA (cDNA) sequencing libraries were constructed using Illumina ®TruSeqTM stranded total RNA sample preparation kit (Illumina, USA) according to the manufacturer's instructions. The concentration and purity of the cDNA libraries were measured (RNA integrity number score > 7.0) using Agilent 2100 bioanalyzer (Agilent Technologies, USA). The libraries were sequenced using the Next-Seq 500 platforms in single-end 150-bp mode (Illumina, Singapore) according to the manufacturer's protocol. The FASTX-Toolkit (25) was used to remove adaptor sequences and to filtered low-quality base call and low-quality reads. TopHat2 program (26) was used to mapping short filtered sequencing reads to the human genome (UCSC) in order to identify exon-exon splice junctions, and the quantified genes expression level was done using Subreads package feature counts (27). They were tested after calculating the library size and appropriate data set dispersion depending on gene expression values by using edgeR Bioconductor package (28). Differentially expressed genes (DEGs) were presented as log fold change (logFC) and p ≤ 0.05 was counted statistically significant.

Gene Ontology and Pathway Analysis of Transcriptomic Data

Gene ontology (GO) analysis (29) was carried out to determine the functions of the DEGs identified. The highly significantly DEGs (p ≤ 0.001) were uploaded into Enrichr tool (30) for analysis and organized into groups basis of cellular components, biological processes, and molecular functions. The pathway enrichment analysis was also conducted for the highly significantly DEGs to place these target genes in the pathways according to the Kyoto encyclopedia of genes and genomes database (KEGG) database (31) by used Enrichr tool.

RNAs Co-expression and Interaction Network

The co-expression network of the most highly significant DEGs (p ≤ 0.0001) was constructed to identify the potential DEGs interaction by using GeneMANIA tool (32). A co-expression network was constructed according to the correlation analysis between the DEGs associated with obesity and BC.

Validation of the Transcriptomic Data by Quantitative Real-Time PCR

Total RNA (800 ng) was reverse transcribed using QuantiTect reverse transcription kit (Qiagen, UK). The expression of a selected gene(s) was measured in duplicate in a large cohort of BC blood patients (33 non-obese and 36 obese BC patients) by quantitative real-time PCR (qRT-PCR) using IQ SYBR green mix (Bio-Rad, USA) and RPL11 as the internal control on CFX Connect™ real-time PCR detection system (Bio-Rad, USA). The primers of selected genes were designed over two different exons and the sequences are available upon request. No-reverse transcriptase controls (NRCs) and no-template controls (NTCs) were included for each primer pair. The relative expression quantification was calculated depending on the 2−ΔΔct methods (33).

Metabolome Analysis

Samples Preparation and Metabolite Extraction

Among the 21 BC patients that subjected in RNA-seq assay; six obese BC and equal numbers of non-obese BC patients were selected for metabolomics. Serum metabolites were extracted and analyzed in triplicates, by liquid chromatography-tandem mass spectrometry (LC-MS/MS), for untargeted metabolomics detection. One-hundred microliters of serum used for metabolites extraction by ice-cold Methanol: Acetonitrile: water in a ratio of (2:1:1 volume/volume), that was added to serum and mixture was vortexed, followed by incubation at −20°C for an hour, samples were then spun at 4°C for 5 min at 8,000 rpm. Supernatants analyzed by LC-MS/MS. Each sample (10 μL) was injected individually into Hypersail gold high-performance liquid chromatography (HPLC) column (150 × 4.6 mm, 5 μm) with a flow rate of 0.250 ml/min and mobile phase A 0.1% formic acid in 99.9% acetonitrile formic acid (0.1%, volume/volume) and mobile phase B is 0.1% formic acid in MilliQ. Mass Spec parameter performed as an earlier report (34). Raw data processed using the online XCMS database (35). Isotopic peaks were integrated using the CAMERA (36). Metabolites were searched using the METLIN database (37) and the pathway analysis done with the help of MetaboAnalyst 3.0 (38).

Statistical Analysis

Statistical analyses using unpaired, two-tailed t-tests were performed in GraphPad software Prism version 8.0.1 (GraphPad Software, La Jolla California USA, www.graphpad.com). The data were presented as a mean ± standard error of the mean (SEM). The level of significance was given at p ≤ 0.05.

Results

Transcripts Profiling Changes Between Obese and Non-obese BC Patients

The RNA-seq study included 21 female patients newly diagnosed with BC and before they underwent any treatment. The non-obese and obese BC patients were significantly different in the BMI, waist and hip circumference (Table 1). Conversely, non-obese and obese BC patients did not show any significant differences with general and clinicopathological characteristics (Supplementary Tables S1, S2). The extracted RNA from whole blood subjected for RNA-seq assay. The RNA-seq data detected a total of 31,698 RNA transcripts; among them, a total of 2,372 transcripts were found significantly dysregulated in obese BC patients compared with non-obese BC patients, of which 1,737 upregulated and 635 downregulated transcripts (Table 2). Moreover, 186 DEGs at the highly significance level (p ≤ 0.001) were identified, comprising 156 upregulated and 30 downregulated transcripts (Supplementary Table S3), furthermore, 31 DEGs were found as the most highly significantly (p ≤ 0.0001), of which 23 upregulated and 8 downregulated transcripts (Table 3). Among all identified DEGs; placenta growth factor (PGF) was the most downregulated gene with a logFC of −5.52, while adenylate cyclase type 1 (ADCY1) was the most upregulated gene with a logFC of 5.52 in obese BC patients compared with non-obese BC patients.

Table 1

Baseline characteristics of studied BC patients in RNA-seq analysis.

Parameters	Non-obese BC	Obese BC	p-value
N (%)	10 (47.62)	11 (52.38)
Age (years)	47.70 ± 2.16	49.09 ± 2.74	0.69
BMI (kg/m2)	22.10 ± 0.88	36.82 ± 1.87	<0.0001
Waist circumference (cm)	76.30 ±5.54	99.82 ± 3.61	0.0018
Hip circumference (cm)	90.20 ± 6.49	119.2 ± 3.92	0.0010

Data were presented as mean ± SEM. BMI, body mass index; N, number of samples.

Table 2

The total number of transcripts altered in obese BC compared with non-obese BC patients, the classification based on p-values.

	All transcripts		Upregulated transcripts		Downregulated transcripts		Non-change transcripts
p-value range	Transcripts number	LogFC range	Transcripts number	LogFC range	Transcripts number	LogFC range	Transcripts number
All	31,698	−5.53 to 6.98	19,105	0.01 to 6.98	12,341	−5.53 to −0.01	252
≤0.05	2,372	−5.53 to 6.98	1,737	0.30 to 6.98	635	−5.53 to −0.26	–
≤0.01	851	−5.53 to 6.98	664	0.41 to 6.98	187	−5.53 to −0.42	–
≤0.001	186	−5.53 to 6.98	156	0.62 to 6.98	30	−5.53 to −0.73	–
≤0.0001	31	−5.53 to 5.52	23	2.65 to 5.52	8	−5.53 to −0.95	–

LogFC rounded to two numbers after the decimal point. FC, fold change.

Table 3

The most highly significantly DEGs in obese BC as compared with non-obese BC patients, ordered depending on LogFC.

	Gene symbol	Gene name	Gene type	LogFC	p-value
1	ADCY1	Adenylate cyclase type 1	Coding	5.52	0.0001
2	ASPM	Abnormal spindle-like microcephaly-associated protein	Coding	4.33	<0.0001
3	E2F7	Transcription factor E2F7	Coding	4.18	0.0001
4	GALNT9	Polypeptide N-acetylgalactosaminyltransferase 9	Coding	4.10	0.0001
5	MYH10	Myosin-10	Coding	4.03	<0.0001
6	SEPT3	Neuronal-specific septin-3	Coding	3.90	0.0001
7	CDK1	Cyclin-dependent kinase 1	Coding	3.84	0.0001
8	AP001429.1	LncRNA-AP001429.1	Non-coding	3.74	<0.0001
9	IGFBP2	Insulin-like growth factor-binding protein 2	Coding	3.68	0.0001
10	BUB1B	Mitotic checkpoint serine/threonine-protein kinase BUB1 beta	Coding	3.50	<0.0001
11	CENPF	Centromere protein F	Coding	3.49	<0.0001
12	OLFM4	Olfactomedin-4	Coding	3.47	<0.0001
13	TOP2A	DNA topoisomerase 2-alpha	Coding	3.44	<0.0001
14	TICRR	Treslin	Coding	3.42	0.0001
15	CEP55	Centrosomal protein of 55 kDa	Coding	3.17	0.0001
16	UHRF1	ubiquitin like with PHD and ring finger domains 1	Coding	3.17	0.0001
17	SCN8A	Sodium channel protein type 8 subunit alpha	Coding	3.08	0.0001
18	SLCO4A1	Solute carrier organic anion transporter family member 4A1	Coding	3.02	<0.0001
19	CD109	CD109 antigen	Coding	2.97	0.0001
20	BRCA2	Breast cancer type 2 susceptibility protein	Coding	2.81	0.0001
21	MYB	Transcriptional activator Myb	Coding	2.74	0.0001
22	MKI67	Proliferation marker protein Ki-67	Coding	2.66	0.0001
23	ARHGEF10	Rho guanine nucleotide exchange factor 10	Coding	2.65	<0.0001
24	TIGD3	Tigger transposable element-derived protein 3	Coding	−0.95	0.0001
25	TPST1	Tyrosylprotein sulfotransferase 1,-like	Coding	−1.32	<0.0001
26	VSIG4	V-set and immunoglobulin domain-containing protein 4	Coding	−1.66	<0.0001
27	RNY1	RNA, Ro-Associated Y1	Non-coding	−2.24	0.0001
28	IGLV1-47	Immunoglobulin lambda variable 1-47	Coding	−2.75	<0.0001
29	IGKV1D-16	Immunoglobulin kappa variable 1D-16	Coding	−2.77	0.0001
30	IGHV6-1	Immunoglobulin heavy variable 6-1	Coding	−3.62	<0.0001
31	PGF	Placenta growth factor	Coding	−5.53	<0.0001

LogFC rounded to two numbers after the decimal point. FC, Fold change.

Gene Ontology Enrichment, Pathway Analysis and RNA Co-expression Network of Circulating Transcriptomic Data

GO enrichment and KEGG pathway analysis of the highly significantly DEGs were performed to identify the gene product enrichment in various GOs categories and determine the DEGs functions. Deregulated genes were enriched in 234, 1,058, 126 targets in the GO molecular function, GO biological process and GO cellular component, respectively. As shown in Figures 2A–C, the highest enriched GO's targeted were associated with the mitotic sister chromatid segregation (GO:0000070) in the GO biological process analysis. Meanwhile, the majority of the transcripts were stimulated patched binding (GO:0005113) in the GO molecular function analysis and related to the condensed nuclear chromosome kinetochore (GO:0000778) in the GO cellular component analysis. In the KEGG pathway analysis, there were 125 pathways for which the DEGs were enriched (Supplementary Table S4), the most predominant pathways being the cell cycle, one carbon pool by folate, progesterone-mediated oocyte maturation, vitamin B6 metabolism, homologous recombination, oocyte meiosis, p53 signaling pathway, Fanconi anemia pathway, cellular senescence, and notch signaling pathway (Figure 2D). The co-expression network was constructed to investigate the potential interaction among DEGs in obese BC compared with the non-obese BC patients. The DEGs interacted with 46 genes by 2,863 total links (Figure 2E). Therefore, each DEG correlates with a large number of mRNA targets, suggesting that the interconnection between DEGs and mRNAs may related to obesity.

Figure 2

GO enrichment, pathway analysis, and co-expression network of the transcriptomic data. (A–C) GO analysis of DEGs that associated with biological process, molecular function, and cellular component. (D) KEGG pathway analysis for DEGs. (E) Co-expression networks of DEGs. Dysregulated genes interacted with 46 total genes by 2,863 total links. The association sorted by combined score ranking. DEGs, differentially expressed RNAs; GO, gene ontology; KEGG, Kyoto encyclopedia of genes and genomes database.

Validation of the Transcriptomic Data via Investigated OLFM4 Expression Level in Large BC Cohort

Based on our RNA-seq data, the olfactomedin-4 (OLFM4) was chosen as it was among the most highly significant DEGs that altered. In addition, previous studies suggested a role of OLFM4 in immune cells and associate with increased risk of human cancers (39, 40). The expression level of OLFM4 was measured in a large validation cohort of BC patients (36 obese BC vs. 33 non-obese BC). From our RNA-seq data, the OLFM4 had an increased expression level in obese BC as compared with non-obese BC patients (3.47-folds; p < 0.0001). Moreover, the OLFM4 was still highly significantly upregulated (10.74-folds; p < 0.0001) in obese BC compared with non-obese BC patients in a large validation cohort (Figure 3). Therefore, the validation results showed concordance with the RNA-seq trend.

Figure 3

The OLFM4 expression level in the blood of obese BC compared with non-obese BC patients in the sequencing and validation cohort. The gene expression was detected by qRT-PCR and normalized by RPL11 expression; ***p < 0.0001.

Untargeted Metabolomics of Circulating Metabolites of Obese BC and Non-Obese BC Patients

We performed the metabolic study in obese BC patients and compared it with the non-obese BC, selected among the same RNA-seq patients using untargeted LC-MS/MS-based metabolomics. After analyzing the feature peaks, 173 features were detected in ESI+ mode. Two-dimensional principal component analysis (PCA) and three dimensional partial least squares-discriminant analysis (PLS-DA) models score plots of all samples showed no outliers in our study and revealed a significant difference in metabolomics between obese BC and non-obese BC samples (Figures 4A–D). We identified 173 metabolites, of which, 100 were downregulated while 73 metabolites were upregulated in obese BC compared with non-obese BC patients. Among these, 96 metabolites were significantly different (Supplementary Table S5), with 36 upregulated and 60 downregulated metabolites. The Spearman's correlation coefficient of the metabolomics data was evidenced by metabolites self-correlations (Figure 5A). However, in the hierarchical cluster analysis (HCA)-heatmap for the differential metabolites, the obese BC samples clustered and separated from non-obese BC (Figure 5B).

Figure 4

Untargeted metabolomics of obese and non-obese BC patients. (A) Total ion chromatogram of two groups in triplicates. (B) Two dimensional PCA score plot with experimental triplicate between samples. (C) Three-dimensional PLS-DA score plot between individual samples. (D) Significant metabolic features are marked in respective retention time and spot size indicates its abundance. PCA, principal component analysis; PLS-DA, partial least squares–discriminant analysis; OBS, obese BC.

Figure 5

Correlation of the metabolomics data. (A) Each metabolite in the square represent the Spearman's correlation coefficient between (r2) are calculated between each metabolite across all metabolites. Metabolite order is determined as in hierarchical clustering. Self-correlations are identified in red. (B) HCA-heatmap analysis showing positive (red) and negative (blue) comparison metabolites of obese BC (green) and non-obese BC (red) patients. P ≤ 0.05 for all the metabolites OBS, obese BC.

Pathway Enrichment Analysis of Deregulated Circulating Metabolites

The differential metabolites between obese BC and non-obese BC samples were used for pathway enrichment analysis. A total of 56 metabolic pathways were shown to be enriched in obese BC compared with non-obese BC patients (Supplementary Figure S1 and Supplementary Table S6), mainly involved in lipid, carbohydrate, and amino acid metabolism. As well as oxidative phosphorylation, and some other metabolic pathways, such as urea cycle, ammonia recycling, vitamins metabolism, etc. (Figure 6), which play important roles in ATP generation and cancer cell proliferation and metastasis, therefore, some can be utilized as novel therapeutic targets for cancer therapy (41).

Figure 6

The top 50 enriched pathway analysis of significant differential accumulated metabolites between obese and non-obese BC patients.

The most important upregulated functionally metabolites were related to epigenetic as well as metabolic pathways that are involved in energy metabolism and cell proliferation such as amino acid and citric acid cycle. Furthermore, the upregulated neurotransmitters metabolites as serotonin, histamine, and acetylcholine may play a different role in the immune system (Table 4).

Table 4

Metabolites identified in obese BC compared with non-obese BC patients.

Category	Compound	Log2FC	p-value
Metabolite involved in epigenetic	Ornithine	2.42	<0.0001
	7-Methyladenine	2.40	<0.0001
	2-Oxoarginine	1.74	<0.0001
Metabolite involved in citric acid cycle	L-Carnitine	1.83	<0.0001
	Glycerol 3-phosphate	2.38	<0.0001
	L-2-Hydroxyglutaric acid	0.71	0.0004
	Pentanoyl-CoA	−1.95	0.0001
Amino acid metabolism	Ornithine	2.42	<0.0001
	2-Oxoarginine	1.74	<0.0001
	Serotonin	0.73	0.0003
	Histamine	0.76	0.003
	Tryptophan	0.42	0.024
	L-Homoserine	0.44	0.038
	D-Leucine	2.34	<0.0001
	3-Hydroxyphenylacetic acid	−4.97	<0.0001
	Carbamoyl phosphate	−1.36	0.0001
	FAD	−2.28	0.0002
	Epinephrine	−0.33	0.0003
	Creatinine	−0.40	0.01
Cholesterol and fatty acid metabolites	25-Hydroxycholesterol	2.31	<0.0001
	Cholestenone	−7.85	<0.0001
	TG[16:0/14:1(9Z)/18:4(6Z,9Z,12Z,15Z)]	−1.52	0.0006
	Alpha-Linolenic acid	−0.77	0.002
	Hexanoyl-CoA	−1.71	0.0004
Neurotransmitters	Epinephrine	−0.33	0.0003
	Serotonin	0.73	0.0003
	SM[d19:1/24:1(15Z)]	6.55	0.0008
	Histamine	0.76	0.003
	Acetylcholine	0.55	0.007

FC, fold change.

Integration Analysis of Transcriptomic and Metabolomic Data

To provide more comprehensive understanding for the association between obesity and BC, the transcript–metabolite interaction network was generated for DEGs and the deregulated metabolites, in obese vs. non-obese BC samples. This provides a visualization of the interactions between functionally related metabolites and genes identified from transcriptomics and metabolomics. The gene–metabolite interaction consists of 65 nodes connected via 91 edges (Figure 7A). Furthermore, integration analysis at the pathway level was undertaken. Seven pathways were enriched during the integration of both transcriptomics and metabolomics data, that includes glutathione metabolism, glycine and serine metabolism, valine, leucine, and isoleucine degradation, purine metabolism, pyrimidine metabolism, thyroid hormone synthesis, and vitamin B6 metabolism (Figure 7B). The DEGs that relate to the integration pathways including RRM2, PSAT1, ADCY1, PAICS, TYMS, and BCAT1 were significantly upregulated, whereas the GPX3 gene was significantly downregulated. While, the differentially accumulated metabolites FAD, L-leucine, and carbamoyl phosphate, were significantly downregulated, whereas ornithine, dihydrouracil, and thymine metabolites, were significantly upregulated in obese compared to non-obese BC samples (Supplementary Table S7). Overall, the integration analysis successfully identified pathways and its related metabolites that can widely affect the functions of immune cells in obese BC patients.

Figure 7

Integration of transcriptomic and metabolomic data. (A) The transcript–metabolite interaction network of the integration network consists of 65 nodes connected through 91 edges. Nodes in orange indicated differential metabolites and nodes in gray indicated DEGs related to each metabolite. (B) Venn diagram of the transcriptomic and metabolomic enriched pathway.

Discussion

Transcriptomics and metabolomics reflect changes in genotype and phenotype, respectively and provide complementary information about genetic alteration, protein synthesis, metabolisms and cellular function (42, 43). Many studies focused on the differential transcripts and metabolites and their functional attributes to understanding the disease's biological interaction. It has been revealed that differentially expressed genes that lead to bio-fluid metabolome change may significantly contribute to the initiation and progression of many types of diseases including obesity and BC (44–47). Currently, the most widely techniques that are used to differentiation in the transcriptomic and metabolomic profiles are RNA-seq, LC-MS/MS, respectively (48, 49). Over the last few years, several transcriptomic and metabolomic studies identified the variation in transcripts and metabolites profiles related with diseases such as obesity and BC compared with non-obese and healthy control cases, respectively, to understand the mechanisms of disease initiation and progression as well as to biomarker identification to deduce the onset and progression of cancer (16, 20, 21, 50–55). Furthermore, other studies have performed the integration between transcriptomic and metabolomic data (14, 56–58). Interestingly, the molecular mechanisms underlying the association between obesity and BC risk are not well-understood and still unclear (18). In an effort to reduce the knowledge gap, we performed these experiments to clarify the potential relationship between obesity and BC by the integration of the peripheral blood differential transcriptomic and metabolomic profiles at the pathway level. Noteworthy, to the best of our knowledge, the approach to our study between obese and non-obese BC patients has not been previously applied. Furthermore, the outcomes of the above-referred transcriptomics and metabolomics studies deals with healthy non-obese cases vs. obese patients, as well as with healthy control cases vs. BC patients were different from our findings of transcriptomics and metabolomics investigation as well as their integration, in obese compared with non-obese BC patients.

In this study, during the comparison of gene expression levels among obese and non-obese BC patients, we identified 2,272 significant DEGs, of which 1,737 transcripts were upregulated and 635 transcripts were downregulated in obese BC. GO analysis and co-expression networks of DEGs were performed to delineates the molecular mechanism and identify interactions among the discovered genes. Unique deregulated transcripts were enriched in different cellular pathways such as: cell cycle, one carbon pathway, homologous recombination cellular senescence, P53, and notch signaling pathway, in obese BC patients when compared with non-obese BC patients. Therefore, DEGs might essentially contribute to the initiation and/or development of obesity that may lead to BC initiation and/or progresses. These findings were different from previously reported studies; such as the Merdad et al., transcriptomics study in the tissue of BC patients compared with normal controls that observed downregulated genes associated with lipid metabolism pathway (20). Additional to Sun et al., transcriptomics study that reveled three deregulated long non-coding RNA (lncRNA) (lncRNA-p5549, lncRNA-p21015, and lncRNA-p19461) in the circulation of obese vs. non-obese individuals (51).

The expression pattern of OLFM4 has been differentially reported among tissues type (59). According to the human protein atlas, the OLFM4 was mainly expressed in the gastrointestinal tract, bone marrow, and immune system. Albuquerque et al. (60) demonstrated that OLFM4 was upregulated in obese children. Moreover, OLFM4 was also highly expressed in colon, breast, and lung cancerous tissues (61), where it inhibits apoptosis and promotes cancer cell proliferation, suggesting it may serve as a diagnostic marker or a therapeutic target for human cancers (39). Our findings revealed that OLFM4 was upregulated in the blood of obese BC as compared with non-obese BC patients. Therefore, the OLFM4 was upregulated in blood immune cells in obesity with BC, suggesting that it may play a unique role in immune cells and associated with increased risk of BC (39, 40).

On the other hand, untargeted metabolome analysis was performed in serum samples selected among the same RNA-seq patient populations. We detected a total of 96 circulating metabolites deregulated in obese BC patients when compared to non-obese BC patients, enriched in 56 metabolic pathways, mainly involved in cellular functions regulation and playing an important roles in ATP generation, cancer cell proliferation, and metastasis, therefore, it can be used as novel therapeutic agents for cancer therapy (41). Moreover, the most important upregulated metabolites that act as epigenetic, play an immunoregulatory role and involved in energy metabolism and cell proliferation.

Finally, the integration analysis between transcriptomic and metabolomic data at the pathway level provided a visualization of the interactions between DEGs and differentially metabolites. Collectively, deregulated obesity-associated genes and metabolites involved in changes of pathways, effective cancer cells metabolic programs, and increase BC risk (62). Data integration revealed novel seven uniquely enriched pathways in obese BC patients when compared with non-obese BC patients. The glutathione metabolism as one of the integrated metabolism, utilized FAD in the glutathione biosynthesis, which promotes tumor progression, metastasis, and protects cancer cells (63). Interestingly, downregulation of FAD metabolite and glutathione peroxidase 3 (GPX3) was detected in obese BC patients. GPX3 protein protects the cells against oxidative damage, thereby the low expression of GPX3 was associated with breast carcinogenesis (64). Free amino acids as glycine, serine, and branched-chain amino acid (BCAA) were associated with obesity and various types of cancer including BC. Glycine and serine metabolism provides the essential precursors for proteins and nucleic acid biosynthesis (65–67). The phosphoserine aminotransferase 1 (PSAT1) and branched-chain amino acid transaminase 1 (BCAT1) were upregulated in many carcinoma tissues and associated with cell proliferation (67, 68). PSAT1 and BCAT1, which played an important enzymatic role in serine metabolism and BCAA degradation, respectively, were found to be upregulated in obese BC patients. Ribonucleotide reductase regulatory subunit M2 (RRM2), an enzyme involved in dNTP production, increased dNTP pools and related with purine and pyrimidine metabolism (69) is tumorigenic and upregulated in cancer cells (70), was also found to be upregulated. The dysregulation of the purine metabolism pathway was also demonstrated in the integration of metabolomic and transcriptomic data of BC patients compared with healthy subjects, which might affect the BC progression (57). All these findings, could enhance the integration between transcripts and metabolomics and provided resistance for BC cell to dodge the circulating immune cells of whole blood.

In conclusion, our results demonstrate alteration in pathways at transcriptome and metabolome level in obese BC patients. This may suggest that obesity-associated transcripts and metabolites reveals alteration in metabolic pathway networks and rewire metabolic programs in cancer cells. This information could provide an important tool in research and may add to the knowledge on the molecular interaction between obesity and BC.

There are some limitations in this work, including a small sample size in metabolomic analysis. In addition, the quantity of blood RNA in samples depended mainly on the content of white blood cells. Therefore, the volume of our blood samples may contain limited amount of circulating RNAs. Further control cross-sectional studies using healthy obese and non-obese patients as well as increase in metabolomics sample size, will be conducted in the near future. Finally, the identification of OLFM4 expressing immune cells and the functional role of OLFM4 in immune cells need further investigation.

Data Availability Statement

The original contributions presented in the study are publicly available. This data can be found here: the NCBI Gene Expression Omnibus (GSE148892).

Ethics Statement

The studies involving human participants were reviewed and approved by the Unit of Biomedical Ethics, Research Committee (Document number: HA-02-J-008), King Abdulaziz University Hospital, King Abdulaziz University, Jeddah, Saudi Arabia. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

MH, HC, KA-S, AA-M, MM, and MK designed and coordinated the experiments. SA diagnosed and recruit the patients. KA-S, JA-M, and SA obtained ethical approval, patient consent, and provided the samples. MH, MS, and AD performed the experiments. JA-M performed the pathological screening. KA-S, HC, and AA contributed in laboratory facilitates. MH and MS analyzed metabolomics data. JR and WW analyzed transcriptomics data. MH wrote the original manuscript draft. HC, AA, MK, and MM edited the manuscript. KA-S, AA, SA, HC, and JA-M provided project funding. All authors read and approved the final manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.