This study collected cervical cancer patients who received radiotherapy in Shanxi Bethune Hospital from May 2018 to October 2023. Inclusion criteria included: (a) patients with a newly diagnosed diagnosis of cervical squamous cell carcinoma (IB-IVA) without surgical treatment; (b) all patients were first-time receiving radiotherapy; (c) the patient has complete pathological, radiographic, and radiotherapy dose information; (d) patients who are clinically followed up and diagnosed with concurrent RP; (e) all patients were in good essential condition and had no apparent abnormalities in the rectum after MRI examination before radiotherapy. Some patients were excluded according to the following criteria: (a) a history of other malignancies; (b) lack of pathological, radiographic, and radiotherapy dose data of patients; (c) intolerance to radiotherapy or chemotherapy and failure to complete the treatment plan due to severe acute toxicity during treatment. Finally, 126 of the 356 patients met the criteria and were included in the study.

We followed up the above patients for 1 year and graded the rectal injury of 126 patients according to the European Organization for Research and Treatment of Cancer-Radiation Therapy Oncology Group(RTOG/EORTC) classification criteria (12). We divided the patients into two groups: proctitis symptoms (grade 2-5), and no obvious proctitis symptoms (grade 0-1) ( Figure 1 ).

All patients received external beam radiotherapy(EBRT)(1.8 ~ 2.0Gy/d, 5 times a week) and brachytherapy(BT)(5-7Gy/time, 4-6 times in total).MRI image data and clinical information of the above patients before radiotherapy and within 12 months after radiotherapy were collected, including age, the International Federation of Gynecology and Obstetrics (FIGO) staging of cervical cancer, minimum lymphocyte count value during treatment, radiation dosimetry parameters, the body mass index (BMI) and comorbidity. Among them, the radiation dosimetric parameters specifically include physical parameters related to external irradiation: the maximum dose (D_max), the minimum dose (D_min) received by the rectum, and the percentage of the total volume (V₃₀, V₄₀, V₄₅) irradiated by 30Gy, 40Gy and 45Gy. In order to solve the problem of calculating the equivalent dose of different fractionated radiotherapy, the equivalent dose in 2Gy/f (EQD2) is finally used as the standard equivalent dose (refer to the following formula):

Among them, n is the number of treatments, d is the fractional dose, 10 is the tumor in the a/β value, and 3 is the normal tissue (13). Finally, the minimum dose (D_1cc, D_2cc) accepted by the volume of 1cm³ and 2cm³ rectum after the superposition of internal and external irradiation is calculated.

All MRI images were obtained on a clinical whole-body 3.0T MRI scanner (GE Signa HDXT 3.0T MRI, GE Healthcare, USA) using a phased array 8-channel sensitivity coded abdominal coil. The patient took a supine position to maintain respiratory control. The scanning range covered the entire pelvis, from the upper edge of the iliac crest to the lower edge of the pubic symphysis. Sagittal T1 and T2 images were obtained from the picture archiving and communication system (PACS, Carestream, Canada) and exported in DICOM format. Tumor segmentation was performed manually using ITK-SNAP software, and the entire rectal region on the patient’s MRI image was defined as the region of interest (ROI). To reduce interfering information, a radiologist with 20 years of experience manually outlined the target region to include only structures within the rectum, excluding peripheral vasculature, peripheral tissues, and peripheral organs. Pelvic MRI images of 30 randomly selected patients were re-segmented by another radiologist with 20 years of experience. 2 weeks later, the first radiologist again re-outlined these randomly selected 30 images. Both radiologists were unaware of the patient’s clinical history and pathological information. The repeatability of segmentation was assessed using the intragroup correlation coefficient (ICC), and an ICC >0.8 was considered a good segmentation agreement (14).

Prior to feature extraction, the images were preprocessed to minimize variations due to different MRI scanners. All images were Z-value normalized to ensure that the images had a standard normal distribution, and then the images were resampled to a voxel size of 1 x 1 x 1 mm. Before the analysis, we first use the synthetic minority over-sampling technique (Smote) to preprocess the training set data so that the proportion of cases with different labels in the training set is 1:1. Then the variables with zero variance are excluded from the analysis, and the median is used to replace the missing values. Finally, the data is standardized. The pyradiomics in Python is used to extract imaging features, including first order features, shape features, texture features and transform features. After the feature extraction is completed, the patient label of ≥ 2 level is set to ‘ 1 ‘, and the patient label of < 2 level is set to ‘ 0 ‘. Then, the following formula is used to calculate the absolute change of features on the image after and before radiotherapy: delta-radiomics features represented by △RF. The calculation formula is as follows:

RF_MRI2 represents the radiomics features of MRI images within 12 months after radiotherapy, and RF_MRI1 represents the radiomics features of MRI images before radiotherapy.

Statistical analyses of the extracted raw features were performed using Python, starting with feature dimensionality reduction using logistic regression (LR) analysis to select features that were significantly different between the two groups. The p-value is usually set at p < 0.2 but can also be set at p < 0.05 or p < 0.1. This requires the researcher to adjust the p-value according to the sample size (15). Due to the limited amount of data in this study, we set p < 0.05 as the threshold. This was then filtered using pairwise Pearson correlation analysis, which sets the Pearson correlation coefficient to 0.7, i.e., any two features were correlated. When their correlation coefficient is greater than 0.7, a feature is removed (16). Finally, to further reduce overfitting and selection bias, the least absolute shrinkage and selection operator (LASSO) was chosen to reduce the coefficients of most of the uncorrelated Δ RFs to zero once the optimal λ had been determined by 5-fold cross-validation, where the maximum area under the curve is the final value of λ. The coefficients of the Δ RFs were then reduced to zero (17).The process of generating and selecting radiomic features was illustrated in Figure 2 .

Random forest (RF) was used to construct the traditional radiomics model and delta-radiomics model. A clinical model was constructed based on the selected clinical features. Finally, statistically significant clinical indicators were included in the radiomics and delta-radiomics models, and the imaging + clinical joint model was constructed, respectively.

SPSS 26.0 software was used for statistical analysis. Quantitative data conforming to normal distribution were expressed as x ¯ ± s. T validation was used to compare the differences between groups. Qualitative data were expressed as frequency (%), and the χ2 validation was used to compare the differences between groups. Python 3.10 was used for feature screening and RF model establishment. The receiver operating characteristic (ROC) curve was used to judge the performance of the machine learning model, and the area under the curve (AUC), sensitivity and specificity were calculated. The DeLong validation is used to evaluate the performance differences of different models in the classification task, and the prediction model’s calibration curve and decision curve analysis (DCA) are constructed. Finally, the Shapley Additive Explanations (SHAP) values are used to analyze the prediction results of the model. P < 0.05 was considered statistically significant.

All samples were randomly divided into training and validation sets in a 7:3 ratio. After processing the training set data using the SMOTE method, 16 samples labeled as ‘1’ were added. The final training set included 104 cases (52 cases in the ‘1’ group and 52 cases in the ‘0’ group), while the validation set contained 32 cases (16 cases in each group). There were no significant differences in baseline characteristics (e.g., age, BMI, comorbidity, etc.) between the two groups (all p-values > 0.05), indicating that the groups were well-balanced at the start of the study ( Table 1 ).

To identify characteristics associated with the severity of RP, we conducted t-tests and χ² tests. The results indicated significant differences in dosimetric parameters, including V₄₀, V₄₅, D_max, D_1cc, D_2cc, and the minimum lymphocyte count, across different severity groups (p < 0.05). Specifically, a lower minimum lymphocyte count correlated with increased severity of RP, while lower values of parameters such as V₄₀ and D_max were associated with milder RP. Notably, V₃₀, D_min, age, FIGO stage, BMI, and comorbidity did not demonstrate statistically significant differences ( Table 2 ).

A total of 1072 features were extracted from the T1 and T2 sequences before radiotherapy, with 12 features ultimately selected after a stepwise screening process. Similarly, 1072 features were extracted from the T1 and T2 sequences after radiotherapy, and 12 features were retained following gradual screening. Using the delta-radiomics calculation formula, a total of 1702 delta radiomic features (△RFs) were obtained, and 329 features remained after logistic regression analysis. Following Pearson correlation analysis, 64 features were retained, and finally, the top 10 radiomics features were selected using the LASSO.

Subsequently, we used RF to develop six prediction models: the radiomics models (model 1, combining T1 and T2 sequences before radiotherapy; model 2, combining T1 and T2 sequences after radiotherapy; model 3, the delta model), the clinical model (model 4), the post-radiotherapy T1 and T2 sequences combined with clinical features (model 5), and the combined clinical and delta-radiomics models (model 6). In the training set, the predictive efficiency of the combined model, consisting of clinical indicators and radiomics features, outperformed the radiomics-only model, with an AUC, sensitivity, and specificity of 0.99, 0.981, and 1.000, respectively ( Figures 3A, B ). Similarly, in the validation set, the combined model outperformed the other models, and the AUC, sensitivity, and specificity data are provided in Table 3 .

We subsequently applied DeLong’s test to compare the differences in predictive power between the models. In the validation set, there were statistically significant differences in AUC between model 1, model 2, and model 6, whereas no significant differences were observed among the remaining models ( Table 4 ).

In the calibration curve, the Brier score is used to evaluate the accuracy of the prediction model. It was observed that the Brier score of model 6 was closest to 0, indicating the highest calibration ( Figure 3C ). DCA demonstrated that models 3 and 6 had significant clinical utility in predicting the severity of RP ( Figure 3D ).

We calculated and visualized SHAP values for each feature in the radiomics model. Figure 4A visually illustrates the 20 key features of the RF model, including 5 dosimetric parameters and 12 imaging features. Unique color dots indicate the impact of each feature: red indicates higher feature values, and blue indicates lower wind feature values. The bar chart in Figure 4B illustrates the average absolute value of the SHAP values for each feature. The above figure shows that D_1cc is the most essential feature for the RF model in predicting the severity of RP occurrence. The higher the feature values of D_1cc, T1_wavelet-LLL_glcm_MCC, and D_2cc, the greater the positive impact on the model output, while the lower the feature values T2_original_firstorder_90Percentile and T1_wavelet-LHH_firstorder_Maximum, the greater the positive impact on the model output.