This single-center retrospective study received approval from our institutional review board and complied with ethical regulations. The requirement for informed consent was waived for this retrospective study of anonymized data.

With the approval of institution ethics committee, we collected two independent patient cohorts. The detailed inclusion and exclusion criteria are shown in Figure 1 . The lung metastasis cohort of 157 RC patients with lung metastasis risk diagnosed between Jan 2018 and Jun 2023. Inclusion criteria were as follows: (a) pathological confirmation of rectal adenocarcinoma; (b) availability of pre-treatment pelvis MRI scans before the initiation of therapy; and (c) utilization of high-resolution contrast enhanced chest CT scans for lung metastasis diagnosis. The exclusion criteria included: (a) concurrent others primary malignant neoplasms or previous anticancer treatment; (b) missing MRI data or insufficient image quality; (c) missing or incomplete electronic medical records; and (d) simultaneous occurrence of other DM.

Besides, we collected a non-overlapping cohort of 662 patients diagnosed with RC at the same institution between Jan 2018 and Jun 2023, forming a distinct non-metastasis cohort. Inclusion criteria (a) and (b) mirrored those of the lung metastasis cohort, while criteria (c) were omitted in the non-metastasis cohort due to the absence of metastasis in these patients. The exclusion criteria remained consistent with those applied to the lung metastasis cohort.

Clinical information, including age, gender, pre-treatment histologic grade (24), and CEA levels, was extracted from medical records. It is important to note that CEA levels were determined through routine blood tests conducted within one week before treatment.

MRI was conducted using a 1.5-T MR system (Optimal 360, GE Healthcare, Waukesha, Wis) or a 3.0-T MR system (MR 750w, GE Medical System) with a phased-array body coil (eight-channel and sixteen-channel phased-array body coil). The standard procedure included axial oblique T2WI sequences and transverse DWI sequences. For the T2WI sequences, the in-plane pixel spacings ranged from 0.366 mm to 0.703 mm, with an average of 0.490 mm, and slice thicknesses ranged from 3.500 mm to 7.000 mm, with an average of 4.692 mm. For the DWI sequences, the in-plane pixel spacings ranged from 0.976 mm to 1.953 mm, with an average of 1.526 mm, and slice thicknesses ranged from 3.500 mm to 7.000 mm, with an average of 4.694 mm. Among the 819 patients in the lung metastasis cohort and non-metastasis cohort, images obtained from different scanners were randomly distributed in the training, validation, and independent test sets.

Radiologists with over 10 years of experience in MRI manually delineated the entire RC tumor using ITK-SNAP 3.9 on pre-treatment T2WI and pre-treatment DWI at b=1000 s/mm². The resulting tumor masks were cropped into image patches as the inputs of 3D networks. To mitigate the impact of variability in acquisition and sequence parameters, image pre-processing was implemented before analysis. All MR images were used the Simple-ITK toolkit to correct the bias field artifacts (25). Gray-level normalization was applied to harmonize the gray values of MR images, compensating for intensity variations across different MRI scanners.

The architecture of our proposed transformer-based DL model for predicting RCLM is illustrated in Figure 2 . We used two pre-trained UniMiSS models, which were built upon the ViT framework, as our primary feature extraction network. T2WI and DWI scans underwent individual processing through each branch network, and the resulting features were fused through concatenation using a post-fusion approach. This fusion strategy engendered a comprehensive representation, combining complementary information from both T2WI and DWI. To distill spatial information, an adaptive max pooling operation was employed to reduce the feature length to 1000. These features were then concatenated with clinical information, including CEA level, age, and gender. The combined features were subsequently fed into a fully connected layers to effectively learn the imaging and clinical information for accurate predicting RCLM. This architecture enables the model to leverage both imaging and clinical data synergistically to enhance prediction performance.

The performance of the models was evaluated on an independent test set, which was not included in the model training. The remaining data were randomly divided into training and validation sets for parameter tuning. Once optimal model parameters were determined using the training and validation sets, the models were evaluated on the independent test set. The models were implemented using PyTorch v1.7. Our models process 3D volumes, with a data batch size of 32 and input image dimensions of 128 (height) × 128 (width) × 16 (depth). In our preprocessing pipeline, we use manually annotated RC tumor masks to crop image patches from the 3D volume data to a specified size, while ensuring that the tumor region is fully retained. During the cropping process, the cropping area is dynamically adjusted based on the mask to ensure that the entire tumor is included within the cropped region. All MR images were linearly transformed to the range of [-1, 1] via gray-level normalization. Data augmentation techniques, including random rotation in the range of [-30, 30] and random scaling in the scale of [0.95, 1.05], were applied for model training to enhance generalization. For optimization, the Adam optimizer with an initial learning rate of 0.0002 was employed, and the models were trained for 100 epochs. Cross-entropy loss was used as the loss function for model training. To ensure a fair comparison, all compared models were trained using the same image size, learning rate, and number of iterations.

Among these, the model using only clinical characteristics employed the Logistic Regression model. The remaining five models of seven combinations of MRI sequences and clinical features used the pre-trained UniMiSS model (built on the ViT framework) to extract image features. Specifically, the “Images-only” model concatenated the DWI and T2WI images along the channel dimension before feeding them into the feature extraction network. When clinical features are included, a post-fusion method was applied to combine the image and clinical features. In addition, the two demographic characteristics (age and gender) and one CEA used in the comparison experiment were pre-treatment. Gender was labeled with 0 and 1, age was standardized with z-score, and CEA was scaled using minimum-maximum normalization. The latter two characteristics of information were subsequently transformed into one-dimensional representations via Word Embedding.

For a thorough evaluation, our method is compared with state-of-the-art deep learning models, including ResNet18, EfficientNetb0, MobileNet, ShuffleNet, and DenseNet. ResNet18 employs residual blocks with skip connections, enabling the network to learn residual functions. EfficientNetb0 uses a compound scaling method to balance depth, width, and resolution, optimizing performance and computational efficiency. MobileNet leverages depthwise separable convolutions to reduce complexity. ShuffleNet introduces a channel shuffle operation to enhance the efficiency of group convolutions while minimizing computational costs. DenseNet incorporates dense connections among all layers, fostering feature reuse and improving gradient flow, resulting in reduced parameters and enhanced training efficiency.

All well-trained models were evaluated on both validation and independent test sets. The performance of the prediction model was assessed using various metrics, including the area under the receiver operating characteristic curve (AUC), accuracy (ACC), sensitivity (SEN), specificity (SPE) and the wilcoxon rank-sum test for statistical comparison. The optimal threshold for the AUC value was determined by maximizing the sum of the sensitivity and specificity values. The 95% CIs were obtained using bootstrapping to assess variability, and p < 0.05 considered indicative of statistically significant difference.

Applying specific inclusion and exclusion criteria, our study enrolled a total of 819 patients, comprising 524 males and 295 females, with an average age of 57.25 ± 12.14 years (ranging from 17 to 86 years). The average CEA level among the patients was 12.58 ng/mL (ranging from 0.50 to 99.62 ng/mL). Among these patients, 657 were randomly assigned to the training set, 81 to the validation set, and 81 to the independent test set. Of the total patients, 157 patients (19.17%) were RC with lung metastasis, while 662 patients (80.83%) were RC without metastasis. In terms of tumor staging, there were 19 patients at the T1 stage (with 1 lung metastasis), 122 patients at the T2 stage (with 9 lung metastasis), 485 patients at T3 stage (with 82 lung metastasis), and 193 patients at T4 stage (with 65 lung metastasis). Regarding the N stage, there were 271 patients at the N0 stage (with 31 lung metastasis), 282 patients at the N1 stage (with 54 lung metastasis), and 266 patients at N2 stage (with 72 lung metastasis). The details of the demographic and clinical characteristics of patients were presented in Table 1 .

We developed a transformer-based model for the prediction of RCLM. As shown in Table 2 , the transformer-based model achieved an Area Under Curve (AUC) of 84.24% (95% CI, 73.87%-92.68%) on the validation set and 83.74% (95% CI, 72.60%-92.83%) on the independent test set. Our model outperformed the performance of ResNet18, EfficientNetb0, MobileNet, ShuffleNet, and DenseNet. Figures 3A, B illustrates the Receiver Operating Characteristic (ROC) curves for six deep learning models. These models yielded AUC values of 70.30% (95% CI, 56.18%-82.92%, p = 0.0196), 72.02% (95% CI, 58.73%-83.08%, p = 0.0186), 72.32% (95% CI, 56.36%-86.31%, p = 0.0158), 75.96% (95% CI, 61.84%-88.36%, p = 0.0495), 77.17% (95% CI, 64.26%-87.92%, p = 0.0412), and 83.74% (95% CI, 72.60%-92.83%) for ResNet18, EfficientNetb0, MobileNet, ShuffleNet, DenseNet and our method on independent test set, respectively.

Table 3 presented the performance of our transformer-based model trained with seven combinations of MRI sequences and clinical features. These combinations include single image models (DWI-only or T2WI-only), models using a single image and clinical features (DWI-Clinical features or T2WI-Clinical features), a model using only clinical features (Clinical features-only), an image-only model with DWI and T2WI (Images-only), and a model using both DWI, T2WI, and clinical features (Images-Clinical features (Ours)). Incorporating two demographic characteristics (age and gender) and one clinicopathologic factors (CEA) into our model achieved the best performance. As outlined in Table 3 , with increases of 2.86% in AUC, 23.07% in accuracy, and 30.00% in specificity compared to the original Image-only model on the independent test set, the combination of T2WI, DWI, and clinical features achieved the AUC of 83.74% (95% CI, 72.60%-92.83%) on the independent test set. Figures 3C, D illustrate the ROC curves for seven different combinations of image and clinical features. These combinations yielded AUC values of 56.57% (95% CI, 41.69%-71.22%, p = 0.0233), 61.62% (95% CI, 43.67%-78.57%, p = 0.0046), 72.73% (95% CI, 60.40%-83.75%, p = 0.0124), 69.70% (95% CI, 52.96%-84.55%, p < 0.001), 51.52% (95% CI, 37.73%-64.58%, p < 0.001), 81.41% (95% CI, 69.93%-90.72%, p = 0.0143), and 83.74% (95% CI, 72.60%-92.83%) for DWI-only, T2WI-only, DWI-clinical characteristics, DWI-clinical characteristics, clinical characteristics-only, images-only, and images-clinical characteristic on independent test set, respectively.

Table 4 and Figure 4 presented comprehensive results on both the validation set and the independent test set, focusing on different T and N stages. In Table 4 , we provided a directly overview of various models’ performance at N stages. Furthermore, Figure 4 visually illustrated the AUC, accuracy (ACC), sensitivity (SEN), and specificity (SPE) for stages T3 and T4 in different models. Our method outperformed the aforementioned state-of-the-art DL methods in terms of AUC in general. Specifically, we achieved an AUC of 76.32% (95% CI, 59.80%-89.55%) in the T3 stage, an AUC of 96.67% (95% CI, 87.14%-100%) in the T4 stage, an AUC of 85.00% (95% CI, 66.67%-100%) in the N0 stage, an AUC of 54.46% (95% CI, 29.89%-82.76%) in the N1 stage, and an AUC of 96.83% (95% CI, 88.67%-100%) in the N2 stage on the independent test set.

Three experts (Li X, Li B, Wan Y, with 32 years, 10 years and 4 years of experience, respectively) dedicated to imaging diagnosis based on MRI data in the independent test set. They consecutively and independently evaluated the MRI data from the independent test set. The diagnostic performance of subjective evaluation by three radiologists was presented in Table 5 . The results of our model had been binarized for fair comparison. Table 5 showed performance with AUCs of 56.52% (95% CI, 45.30%-69.65%, p < 0.001), 56.82% (95% CI,44.16%-70.89%, p < 0.001), 62.27% (95% CI, 47.88%-76.32%, p < 0.001) and 83.74% (95% CI,72.60%-92.83%) for three experts and our model. The diagnostic time for each case was 22.32s, 38.62s, 11.22s, and 0.67s for three experts and our model, respectively. As shown in Table 5 , the experts’ results showed largely individual differences. In contrast, our model achieved the best performance in predicting RCLM.