Before filtering:

After filtering:

Differentially expressed genes subset shows the best clusterization of PAM50 subtypes

First 2 components of PCA explain only ~16% of variance

t-SNE better distinguishes Luminal A and Luminal B subtypes

Adding PCA (`n_components` explaining 80% of variance) before t-SNE led to a more distinguished representation of subtypes

UMAP captures global tendencies better than t-SNE

            classifiers_degs = {
                'Random Forest': RandomForestClassifier(random_state=42),
                'SVM': make_pipeline(StandardScaler(), SVC(probability=True, random_state=42)),
                'Gradient Boosting': HistGradientBoostingClassifier(random_state=42),
            }
            

            param_space_gb = {
                'learning_rate': [0.01, 0.1, 0.2],
                'max_iter': [100, 200, 300],
                'max_leaf_nodes': [21, 31, 41]
            }

            grid_search_gb = GridSearchCV(
                classifiers_degs['Gradient Boosting'],
                param_space_gb,
                n_jobs=-1,
                scoring='f1_weighted',
                verbose=1,
            )

            grid_search_gb.fit(X_train, y_train)

            print(grid_search_gb.best_params_)
            print(grid_search_gb.best_score_)
          

            param_space_rf = {
                'n_estimators': [100, 200, 300],
                'max_features': ['sqrt', 'log2', None],
                'max_depth': [None, 10, 20, 30]
            }

            grid_search_rf = GridSearchCV(
                classifiers_degs['Random Forest'],
                param_space_rf,
                n_jobs=-1,
                scoring='f1_weighted',
                verbose=1,
            )

            grid_search_rf.fit(X_train, y_train)

            print(grid_search_rf.best_params_)
            print(grid_search_rf.best_score_)
          

            param_space_svm = {
                'svc__C': [0.1, 1, 10],
                'svc__gamma': ['scale', 'auto'],
                'svc__kernel': ['linear', 'poly',
                                          'rbf', 'sigmoid']}

            grid_search_svm = GridSearchCV(
                classifiers_degs['SVM'],
                param_space_svm,
                n_jobs=-1,
                scoring='f1_weighted',
                verbose=1,
            )

            grid_search_svm.fit(X_train, y_train)

            print(grid_search_svm.best_params_)
            print(grid_search_svm.best_score_)
          

              classifiers = {
                  'rf1': RandomForestClassifier(
                      **{'max_depth': None, 'max_features': None,
                          'n_estimators': 100}),
                  'rf2': RandomForestClassifier(
                      **{'max_depth': None, 'max_features': None,
                          'n_estimators': 300}),
                  'gb': HistGradientBoostingClassifier(
                      **{'learning_rate': 0.1, 'max_iter': 200,
                          'max_leaf_nodes': 31}),
                  'svm': make_pipeline(StandardScaler(), SVC(
                      **{'C': 0.1, 'gamma': 'scale', 'kernel': 'linear'}))
              }

              rkf = RepeatedKFold(n_splits=5, n_repeats=5,
                                  random_state=42)

              scores = {}
              for name, clf in (pbar := tqdm(classifiers.items())):
                  pbar.set_description(name)
                  scores[name] = cross_val_score(
                      clf, X_train, y_train, cv=rkf,
                      n_jobs=-1, verbose=2,
                      scoring='f1_weighted')
            

Choosing Random Forest with 100 estimators, as its performance is not significantly worse
than RF with 300 estimators or Gradient Boosting, but it is less computationally expensive

• Based on the mean decrease in impurity (Gini importance)
• Computationally inexpensive
• Easy to interpret
• Biased toward high cardinality features
• Unequally distributes importance between highly correlated features
• Based on feature permutations
• Measures the direct impact of each feature
• Easy to interpret
• Reduces importance of highly correlated features*
• Computationally expensive
• Based on SHAP (SHapley Additive exPlanations) values
• Provides a consistent and accurate measure of feature importance
• Accounts for feature interactions and provides local interpretability
• Computationally expensive
• More complex

* The impact of correlation was reduced by clustering features using Spearman correlation and silhouette scores,
retaining only a core feature from each cluster.

TPX2, CENPA, ESR1, KRT14, SGOL1, KRT5 (highlighted in bold)
were ranked among the top 10 features in at least 2 feature importance evaluation techniques

Luminal vs other:

• TPX2 is associated with poor prognosis in luminal breast carcinoma; knockdown of TPX2 impacted the cell cycle and colony formation in MDA-MB-231 cells (Kahl et al., 2022).
• Mutations in ESR1, particularly in the ligand-binding domain, lead to endocrine therapy resistance by causing constitutive receptor activity, contributing to poorer patient outcomes in luminal breast cancer (De Santo et al., 2019).
• In many breast cancer models, luminal cells that express basal cell markers, including KRT14, become the leaders of invasion (Khalil et al., 2024).
• KRT5 and KRT14 expression in luminal breast cancer are associated with a basal-like phenotype, tumor heterogeneity, and poor prognosis due to their roles in cancer stem cell properties (Li et al., 2022).

Luminal A vs B:

• CENPA was identified in an mRNA signature correlated with a lower survival ratio in luminal A breast cancer and not in luminal B (Xiao et al., 2018).
• CENPA and FOXM1 are key transcription factors, driving cell cycle, proliferation, and poor outcomes in luminal A breast cancer (García-Cortés et al., 2021).

KRT5, ESR1, TPX2, KRT14, MLPH, SGOL1, NEIL3 (highlighted in bold)
were ranked among the top 10 features in at least 2 feature importance evaluation techniques

• Mostly the same key genes as in luminal A subtype except for CENPA.
• KRT5, SFRP1, and KRT14 were downregulated in all TNBC samples, and these changes are more likely to be characteristic of the luminal B subtype, according to (Ossovskaya et al., 2012).
• EGFR expression was higher in luminal B than in luminal A subtype and was associated with poor prognosis in ER+ breast cancers (Jeong et al., 2019).

• Couldn't find any information regarding biological role of SGOL1 (both luminal subtypes) and NEIL3 (luminal B subtype).
• ESR1, MLPH, FOXA1 distinguish both luminal subtypes from basal-like and Her2-enriched subtypes and will be adressed later.

Only ESR1 (highlighted in bold)
was ranked among the top 10 features in at least 2 feature importance evaluation techniques

MLPH, FOXA1, ESR1, TTC6, HJURP (highlighted in bold)
were ranked among the top 10 features in at least 2 feature importance evaluation techniques

ESR1, TPX2, KRT14, KIF20A (highlighted in bold)
were ranked among the top 10 features in at least 2 feature importance evaluation techniques

• Normal-like breast cancer subtype is still not well understood.
• Given the classifier's poorest performance on this subtype, the most important features identified are unreliable.
• From the top 10 genes identified via all three feature importance techniques only KIF20A was found in the context of normal-like breast cancer subtype:
  KIF20A mRNA levels were a negative prognostic indicator in normal-like breast cancer patients. KIF20A expression correlates with distant metastasis-free survival in patients with luminal B, HER2+ and normal-like breast cancer, but not in patients with basal or luminal A subtype cancer (Mamoor, 2021).

Genes already being used in clinic:

• ERBB2 testing is standard in breast cancer diagnosis and treatment planning (Hilal and Romond, 2016). Targetable ERBB2 mutations are also being explored as markers for personalized therapy(Kurozumi et al., 2020). ERBB2 is part of both BluePrint and Prosigna/PAM50 molecular classification systems for breast cancer subtyping (Krijgsman et al., 2012, Parker et al., 2009).
• CENPA is included in some RNA-based prognostic assays like MammaPrint^® and the Molecular Grade Index (MGI^SM) (Rajput et al., 2011).
• MLPH is part of both BluePrint and Prosigna/PAM50 (Krijgsman et al., 2012, Parker et al., 2009).
• FOXA1 is used in clinical tests to predict prognosis and treatment response (Metovic et al., 2022). FOXA1 is part of both BluePrint and Prosigna/PAM50 (Krijgsman et al., 2012, Parker et al., 2009).
• GRB7 expression is included in the 21-gene recurrence score assay used to predict recurrence in breast cancer (Sparano et al., 2015). GRB7 is part of both BluePrint and Prosigna/PAM50 (Krijgsman et al., 2012, Parker et al., 2009).
• ESR1 mutations are common in metastatic estrogen receptor-positive, HER2-negative breast cancer and are used in clinical tests to predict resistance to estrogen deprivation therapies (link). ESR1 is part of both BluePrint and Prosigna/PAM50 (Krijgsman et al., 2012, Parker et al., 2009).

Genes being investigated as potential biomarkers:

SFRP1 (Wu et al., 2020, Baharudin et al., 2020), EGFR (Kjær et al., 2020), LINC00504 (Hou et al., 2021),
STARD3 (Lodi et al., 2023, Korucu and Inandiklioglu, 2024).

1. Differential Expression Analysis:
• Significant differentially expressed genes (DEGs) were identified when comparing each breast cancer subtype to healthy samples, with Luminal B showing the highest number of DEGs (24.19%) and Luminal A the lowest (17.25%).
• Significant differentially expressed genes (DEGs) were identified when comparing each breast cancer subtype to other subtypes combined, with Basal-like showing the highest number of DEGs (14.46%) and Normal-like the lowest (0.58%).
2. Dimension Reduction and Visualization:
• PCA combined with t-SNE provided the best clustering of PAM50 subtypes, while UMAP captured global tendencies better than t-SNE alone.
3. Classifier Performance:
• Random Forest with 100 estimators was chosen for its balance of performance and computational efficiency, achieving an accuracy of 86% and high precision and recall for most subtypes.
• Most of key genes identified through feature importance analysis were found to be biologically relevant and associated with breast cancer prognosis and treatment response.