Data were extracted from the Surveillance, Epidemiology and End Results (SEER) research database. The SEER program is the most authoritative and premier source of cancer statistics in the US, collecting demographics, tumor characteristics and survival data. The SEER data of the version November 2019 (https://seer.cancer.gov/data-software/documentation/seerstat/nov2019/) were downloaded. The downloaded data contained four compressed files: file1Bvek1.rar, file9c7Clt.rar, fileb0bhOt.rar and fileYhNIxM.rar. Each of the above files contained two files, titled xxx.dic and xxx.txt. The dic file is the column id of the txt file is the content. An in-house python script was used to read and process these files. In brief, a function named extract_colname was used to extract the column names from the dic file, which were assigned to the corresponding txt file. The patients with BC were then selected and the histology was rated according to criteria in the 3rd edition of the International Classification of Diseases for Oncology (10).

Patients with BC diagnosed after 2010 were included because molecular subtypes were available from then on. Only cases with complete data, without missing values on important covariates (age, ethnicity, tumor site, grade, size) were eligible. Cases that were reported from a death certificate or autopsy were excluded and a 2-month latency exclusion was set to further distinguish SPMs from simultaneous cancers. The identified patients with BC were then categorized into two cohorts: The OOPM cohort and the SPM cohort. The study design and workflow are presented in Fig. 1.

According to the SEER rules for classifying multiple primary cancers, the definition was dependent on the cancer site of origin, date of diagnosis, histology, tumor behavior (i.e., in situ vs. invasive) and laterality of any paired organs. In general, all SPMs occurring 2 or more months after the initial diagnosis were considered separate primary cancers unless the medical record stated that the tumor was recurrent or metastatic (11,12). There were also two key variables indicating multiple primary malignancies in SEER, the ‘total number of in situ/malignant tumors for patient’ and the ‘sequence number’ of the multiple primary malignancies. The former was used to identify patients with an SPM and the latter to index the sequence of multiple malignancies. A random variable that was named ‘indicator of SPM’ was defined to indicate whether the patient had developed one.

An unsupervised machine learning method called factor analysis of mixed data (FAMD), which is generally used to analyze datasets containing both quantitative and qualitative variables, was used to transform data. In brief, FAMD may be regarded as a mix between principal component analysis (PCA) and multiple correspondence analysis. It acts as PCA for quantitative variables and as multiple correspondence analysis for qualitative variables. This was achieved by an R package FactoMineR (13).

The total dataset was randomly split 75/25% into the training set and testing set, stratified by the existence of SPMs. A popular supervised machine learning method called random forest was applied to the training set to predict the likelihood of developing SPMs. The performance of the random forest classifier was evaluated in the testing set with 50 repetitions to reduce the influence of randomization. The outcome was visualized by the receiver operating characteristic (ROC) curve. Feature importance was calculated by their contribution to the prediction ability of the model. The above steps were implemented in python using the sci-kit-learn package.

To compare distributions between variables, the χ² test was generally applied for discrete variables, Student's t-test for continuous variables satisfying a normal distribution and the Mann-Whitney U-test for continuous variables otherwise. For survival analysis, both the non-parametric Kaplan-Meier model and the semi-parametric Cox proportional hazard model were used to evaluate the influence of variables on overall survival (OS); when both methods produced a significant P-value, the result was regarded to be significant. P<0.05 was considered to indicate statistical significance.

A total of 286,047 patients with BC were identified from the SEER database. Of them, 26,657 (9.32%) developed SPMs within a maximum follow-up of ~7 years. The characteristics of the patients with OOPM and SPM are compared in Table I. In general, ~99% of patients with BC were females and ~50% of BCs were in the left breast. Specifically, the SPM cohort had significantly more patients with well (23.91%) or moderately (44.22%) differentiated pathology and with more widowed patients (14.87%), while the OOPM cohort had a significantly higher percentage of patients poorly differentiated tumors (31.86%) and of married patients (55.04%). The SPM cohort was significantly older than the OOPM cohort (median age, 63 vs. 60 years; P<0.001). Of note, the SPM frequency was significantly (P=0.003) higher in stage M0 (25,316/245,406, 10.32%) than in stage M1 (1,089/11,590, 9.39%; Table I).

The cancer types of the SPMs are profiled in Fig. 2A. Half of the SPMs occurred in the breast and half were found evenly in the other systems, including the respiratory, female germline, lymphatic/leukocytic, colorectal, other digestive and urinary systems. of note, 0.2% of SPMs were detected in the male germline system.

The histological type distribution in the SPM cohort was compared with that in the OOPM cohort (Fig. 2B). Most SPMs were ductal carcinoma in both cohorts, but the frequency of ductal carcinoma was significantly lower in the SPM cohort (69.1%) than in the OOPM cohort (75.2%, Fisher's exact P<0.001). On the other hand, the frequency of lobular carcinoma, which was the second most common carcinoma in both cohorts, was significantly more common in the SPM cohort (23.3%) than in the OOPM cohort (17.1%, Fisher's exact P<0.001).

Molecular status had been determined by a combination of immunohistochemistry, fluorescence in situ hybridization, chromogenic in situ hybridization and other methods by the SEER group (https://seer.cancer.gov/seerstat/databases/ssf/her2-derived.html). The HER2^-/HR⁺ subtype was significantly enriched in the SPM cohort compared to the OOPM cohort (78.1% vs. 72.8%, P<0.001; Fig. 2C). Of note, HER2^-/HR⁺ was the most common subtype in both the OOPM cohort and SPM cohort, and the HER2^-/HR⁺ subtype was more likely to have an SPM (10%) than other subtypes. The proportion of HER2⁺/HR⁺, HER2⁺/HR^- and triple-negative subtypes was generally lower in the SPM cohort than in the OOPM cohort.

Since it is at times difficult to distinguish metastasis and SPM, the SPM frequency was compared between stage M0 and stage M1, stratified by histological and molecular subtype. Although the SPM frequency was slightly higher in stage M0 than in stage M1 across numerous histological types, no significant difference in SPM frequency was detected in any histological type (Fig. 2D). As for molecular subtypes, the SPM frequency was significantly higher in stage M0 in the HER2^-/HR^- subtype (Fig. 2E).

The OS of the SPM cohort was significantly worse than that in the OOPM cohort (hazard ratio: 1.49; 95% CI: 1.44-1.53; log-rank P<0.001), indicating the role of SPMs in accelerating patient death (Fig. 3A). As more than half of SPMs occurred in the breast again, while the rest seemed to occur in other organs at random, their OS rates were compared to explore the influence of SPM location on survival. The OS of breast SPMs was significantly better than that of other single-organ SPMs (hazard ratio: 0.46; 95% CI: 0.45-0.49; log-rank P<0.001) and that of multiple-organ SPMs (hazard ratio: 0.44; 95% CI: 0.39-0.50; log-rank P<0.001; Fig. 3B), while there was no significant OS difference between the latter two groups (hazard ratio: 0.98; 95% CI: 0.87-1.10; log-rank P=0.69). Median OS was not achieved (i.e. the survival rate remained >50%) in patients with SPMs of the breast, while it was 81 and 76 months in patients with SPMs of other organs and in patients with multiple SPMs, respectively.

Since distant metastasis is a key factor influencing OS, this was validated in the present dataset (Fig. 3C). To better understand the role of distant metastasis and SPM on OS, multivariate Cox survival analysis was performed for distant metastasis, SPM and their interacting effect as covariates. The results indicated that both distant metastasis and SPM were significantly associated with OS after adjusting for each other (for SPM, hazard ratio: 1.71; 95% CI: 1.70-1.82; for M1, hazard ratio: 12.64, 95% CI: 12.38-13.05; Table II). The significant interaction P-value implies that the influence of SPM on OS was more significant in patients with stage M0.

To detect whether the SPM and OOPM cohorts may be distinguished by certain features, unsupervised transformation was performed using FAMD, which was an extension of PCA. The general purpose of PCA is to find transformed features that may cluster the patients into two or more clusters and the transformed features are a combination of the original variables. In the present study, there were 23 original variables. After transformation, the top 5 features were extracted. Only slightly >10% of the variance of the data was able to be explained by the top five transformed features (Fig. S1A), while the top two transformed features were only able to explain 6.5% of the variance. The tumor stage contributed the most to the variance detected by the top two transformed features (Fig. S1B). Patients with SPMs were not able to be clustered together using the top two transformed features (Fig. S1C). Therefore, using supervised learning, a random forest model was created to predict the probability of SPM in patients with BC.

The patient population was randomly split into a training set (75%) and a testing set (25%), each stratified by the presence of SPMs. Parameters including maximum depth and class weight were learned from the training set and the parameters that generated the highest positive predictive value (PPV) in the out-of-bag mode were adopted to create the random forest model. The model generated an overall area under the curve (AUC) of 0.95 in the training set (Fig. 4A). To reduce the influence of randomization, the model was tested 50 times using different random seeds. The mean consistency of the testing set was 0.91, but in the imbalance dataset, the PPV and negative predictive value (NPV) were more important features. The mean AUC, NPV and PPV in the testing set were 0.57, 0.91 and 0.32, respectively (Fig. 4B). The median value of the above parameters was the same as their mean value.

A representative ROC curve with its AUC in the testing set is illustrated in Fig. 4C. The top 10 features contributing to the estimation in the testing set were displayed in Fig. 4D. Age at diagnosis and tumor size had the highest weight, at ~25 and 20%, respectively.