Materials and methods

Study design and data collection

This retrospective cohort study included pregnant women who received antenatal care and delivered at Suzhou Municipal Hospital and the First Affiliated Hospital of Soochow University between 2015 and 2024. The Reproductive Medicine Ethics Committee of Suzhou Municipal Hospital reviewed and approved this research (K-2025-200-K01). All data were collected from the hospitals' electronic medical record systems. Laboratory markers test results collected before 14 weeks of gestation were considered as potential features for the development of the model for each participant. After excluding 140 cases of fetal demise, termination, or spontaneous abortion, 12,715 pregnancies remained for analysis. Given the retrospective nature of the study, informed consent was waived by the institutional review board. All patient information was anonymized in accordance with the Declaration of Helsinki.

Laboratory parameter measurement and statistical analysis

Blood samples for laboratory markers testing were collected during the routine first-trimester screening visit, scheduled between 11 and 13+6 weeks of gestation. The median gestational age at blood sampling for the entire cohort was 12.43 weeks (interquartile range: 11.81–13.10 weeks). In the control group, the median was 12.43 weeks (interquartile range: 11.86–13.10 weeks), while in the PE group, it was 12.29 weeks (interquartile range: 11.71–13.00 weeks); the difference between groups was not statistically significant (p = 0.283) (Table S1).

Data collected included demographic characteristics, laboratory markers test results, and ultrasound screening parameters, all obtained within the same early pregnancy clinical assessment protocol (11–13+6 weeks). Clinical variables included MAP, maternal age, weight, height, gravidity, parity, plurality (e.g., twin vs. singleton gestation), and mode of conception (natural vs. assisted reproduction–in vitro fertilization/IVF). Additional maternal history variables included smoking status, history of hypertension, diabetes, fetal growth restriction, PE, obstetric complications, and medication use during pregnancy. Laboratory markers included the angiogenic marker PLGF; liver function indices: total bile acids (TBA), albumin-to-globulin ratio (A/G), total bilirubin (TBIL), direct bilirubin (DBIL), indirect bilirubin (IBIL), prealbumin (PALB), globulin (GLB), ALT, alkaline phosphatase (ALP), gamma-glutamyltransferase (GGT), lactate dehydrogenase (LDH), and cholinesterase (CHE); metabolic and renal markers: fasting blood glucose (FBG), UA, blood urea, and creatinine (Cr); electrolytes and minerals: total carbon dioxide (CO2), sodium (Na), potassium (K), calcium (Ca), magnesium (Mg), copper (Cu), iron (Fe), chloride (Cl), phosphorus (P), and zinc (Zn); lipid profile: triglycerides (TG), total cholesterol (TC), high-density lipoprotein cholesterol (HDL-C), low-density lipoprotein cholesterol (LDL-C), and very low-density lipoprotein cholesterol (VLDL-C); inflammation marker C-reactive protein (CRP); derived indices: triglyceride-glucose index multiple of the median (TyG MoM) and uric acid-to-albumin ratio multiple of the median (UAR MoM); other markers: total protein (TP), albumin (ALB), AST, hydroxybutyrate dehydrogenase (HBDH), creatine kinase (CK); ultrasound parameters: ultrasound data included the UtA-PI measured during first-trimester screening. Pregnancy outcomes and complications: Information on pregnancy outcomes and obstetric complications was collected from hospital medical records and discharge summaries. Furthermore, missing laboratory or clinical data were imputed using the median value calculated from the available data for each respective variable. This approach ensures that imputed values align with the overall data distribution and helps mitigate potential biases due to missing data. Variables with a missing rate exceeding 20% were excluded from the analysis.

The primary outcome of this study was the development of PE, defined as new-onset hypertension accompanied by proteinuria. Hypertension was defined as a systolic blood pressure ≥ 140 mmHg and/or a diastolic blood pressure ≥ 90 mmHg. Proteinuria was defined by any of the following criteria: at least two readings of ≥ 1+ on a urine dipstick test, ≥ 300 mg of protein in a 24-h urine collection, or a protein-to-creatinine ratio ≥ 30 mg/mmol. In the absence of proteinuria, a diagnosis of PE could also be made if any of the following severe features were present: thrombocytopenia (platelet count < 100 × 109/L), impaired liver function (serum transaminases > 2 times the upper limit of normal), renal insufficiency (serum creatinine > 1.1 mg/dL or > 2 times baseline), pulmonary edema, or new-onset neurological symptoms including visual disturbances. PE was further categorized based on gestational age at delivery into preterm PE (< 37 weeks of gestation) and term PE (≥ 37 weeks of gestation).

The baseline characteristics of the study participants were presented as the median (first quartile, third quartile) for continuous variables and the number (percentage) for categorical variables. The Mann–Whitney U test and Chi-square test were used to assess differences between groups for continuous and categorical variables, respectively. Multivariate logistic regression models (Model 1, Model 2, and Model 3) were employed to evaluate the association between maternal early laboratory markers and PE. Model 1 (Crude model) did not adjust for covariates, Model 2 adjusted for maternal age, body mass index (BMI), and plurality, and Model 3 adjusted for maternal age, BMI, plurality, parity, and IVF conception. Additional trend analyses were also conducted. Laboratory markers were first divided into quartiles, and the medians of each group were then introduced as continuous variables into the logistic regression model to identify potential linear trends in the association. Results were expressed as odds ratios (OR) and 95% confidence intervals (CI). Additionally, we conducted restricted cubic spline (RCS) analysis to assess the nonlinear relationship between laboratory markers and the risk of PE, using likelihood ratio tests to detect nonlinearity. The threshold effects of laboratory markers on PE risk were further analyzed using a two-stage linear regression model.

Machine learning model development and validation

A total of 12,715 pregnant women were included, comprising 12,159 without PE, 262 with preterm PE, and 294 with term PE. Two analyses were performed: one for preterm PE (12,421 women: 262 cases and 12,159 controls) and one for term PE (12,453 women: 294 cases and 12,159 controls). Each dataset was randomly split into training and testing sets (8:2). Feature selection was performed on the training set using the Boruta algorithm based on the random forest classifier. This algorithm identifies relevant features by comparing the importance scores of actual variables with those of shadow features, which are randomly permuted versions of the original variables. As a result, 22 out of 26 differentially expressed features were deemed important and selected for preterm PE, and 19 out of 26 for term PE. These features were incorporated into seven different machine learning models: logistic regression, random forest, k-nearest neighbors (KNN), adaptive boosting (AdaBoost), categorical boosting (CatBoost), extreme gradient boosting (XGBoost), and light gradient boosting machine (LightGBM).

For each model, hyperparameters were optimized using Bayesian optimization to maximize the area under the receiver operating characteristic (ROC) curve (AUC) on the training set. Ten-fold cross-validation was employed to ensure robustness and reliability. ROC curves and sensitivity values were calculated on the testing set to evaluate model performance. The model with the best combination of AUC and sensitivity was selected as the optimal model. To enhance interpretability, SHapley Additive exPlanations (SHAP) were computed to quantify the contribution of each feature to the model's predictions, and SHAP summary plots were generated to visualize feature impacts.

Among the cohort, 5901 pregnant women underwent PE risk assessment using the FMF clinical model, which incorporated maternal risk factors, MAP, PLGF, and UtA-PI. The study population included 80 cases of preterm PE, 153 cases of term PE, and 5668 controls without PE. For model development, two analyses were conducted: one for preterm PE (80 cases and 5668 controls) and one for term PE (153 cases and 5668 controls). The datasets were randomly divided into training and testing sets in an 8:2 ratio. Input variables included maternal risk factors, MAP, PLGF, UtA-PI, and prediction scores from the preliminary models (AdaBoost for preterm PE and LightGBM for term PE). Seven machine learning algorithms were evaluated, and the optimal model was selected based on the AUC and sensitivity.

Cell-free DNA sequencing and gene function analysis

Maternal plasma samples were collected during mid-pregnancy, before any clinical diagnosis of PE, from 51 pregnant individuals undergoing routine first-trimester screening. Of these, 18 individuals were later diagnosed with preterm PE, defined as PE requiring delivery before 37 weeks of gestation, based on criteria from the American College of Obstetricians and Gynecologists (ACOG).

Approximately 5 mL of maternal peripheral blood was drawn into EDTA-containing tubes. Plasma was isolated by sequential centrifugation at 1600 g for 10 min and 16,000 g for 10 min at 4 °C. Plasma aliquots were carefully transferred to fresh Eppendorf tubes and stored at −80 °C until processing.

Cell-free DNA (cfDNA) was extracted from 200 μL of plasma using a commercial Circulating DNA from Plasma Kit (GenMag Biotech, Beijing, China), following the manufacturer's instructions. Sequencing libraries were prepared according to the manufacturer's protocols. Library concentrations were measured using a Qubit 2.0 fluorometer (Invitrogen, Carlsbad, California, USA), and fragment size distribution was verified using the Agilent High Sensitivity DNA Kit and a 2100 Bioanalyzer (Agilent Technologies, Palo Alto, California, USA). Indexed libraries were pooled and sequenced on the MGISEQ-2000 platform using paired-end 100 bp reads. A minimum of 128 million reads was generated per sample.

Sequencing reads were aligned to the human reference genome (hg19) using BWA-MEM.9,10 PCR duplicates were removed using the MarkDuplicates function from Picard tools (version 1.119). Sequencing depth and coverage were assessed using the mpileup function of samtools. Genomic annotations were retrieved from RefSeq (University of California, Santa Cruz).11

To quantify promoter coverage, the ±1 kb region flanking the transcription start site (TSS) for each gene was defined as the promoter region. Read counts were computed in 10 bp bins using bamCoverage to generate bigWig files. The computeMatrix function was used to extract read coverage across each promoter region. To correct for variations in sequencing depth across samples, raw coverage values were normalized using a modified reads-per-million method according to the formula: promoter proofing of TSSi = coverage of TSSi proofing ÷ coverage of total TSSi proofing × 106.

Normalized TSS coverage profiles were compared between preterm PE cases and controls. Genes with significant differential coverage were identified using the Mann–Whitney U test, and Benjamini-Hochberg correction was applied to control the false discovery rate. Gene set enrichment analysis (GSEA) was then performed using the fgsea R package with Hallmark, Kyoto Encyclopedia of Genes and Genomes (KEGG), Reactome, and WikiPathways gene sets sourced from the Molecular Signatures database (MSigDB). Enrichment scores were calculated for each pathway, and significance was determined based on normalized enrichment score (NES), nominal p value, and adjusted false discovery rate < 0.25.

Acquisition of the GEO datasets and differential expression analysis

Three placental tissue-derived gene expression microarray datasets (GSE10588,12 GSE25906,13 and GSE7501014) were retrieved from the GEO database (http://www.ncbi.nlm.nih.gov/geo/) according to the following criteria: sample size ≥ 5 per group (PE vs. controls); clinically confirmed diagnosis of PE based on ACOG criteria; exclusive use of placental tissue samples (excluding cell lines or other tissue). GSE10588, based on the GPL2986 platform, comprised 26 control and 17 PE samples; GSE25906, based on GPL6102, included 37 controls and 23 PE samples; and GSE75010, based on GPL6244, contained 77 controls and 80 PE samples. Data preprocessing, including standardization, probe annotation, and normalization, was performed using the limma R package. The three datasets were combined to form a unified expression matrix. Given that the chip data was measured by different platforms, the sva R package was applied to remove batch effects, and principal component analysis was performed to evaluate sample distribution before and after correction. Differentially expressed genes (DEGs) between PE and control groups were identified using the limma package, applying thresholds of |log2 fold change| > 0.2 and adjusted p value (P.adj) < 0.05. Volcano plots were generated with the ggplot2 package to visualize expression patterns. Subsequently, the identified DEGs were intersected with mid-pregnancy cfDNA-derived DEGs (P.adj < 0.05), with results visualized through a Venn diagram.