Antitumor immunotherapy can enable promising and durable responses following their clinical application. However, heterogeneity in the tumor immune microenvironment leads to differences in the individual response rates. In this study, we identified novel immune-related molecular subclasses of breast cancer using a non-negative matrix factorization analysis. We enrolled 4184 patients with breast cancer, including 1104 patients from The Cancer Genome Atlas as a training cohort and 3080 patients from another four independent datasets as validation cohorts. In the training cohort, 36.9% of patients who exhibited significantly higher immunocyte infiltration and enrichment of immune response-associated signatures were categorized into an immune class, which was confirmed by probing the expression of immunocyte markers (CD3, CD19, and CD163). Within the immune class, 53.3% of patients belonged to an immune-suppressed subclass, characterized by the activation of stroma-related signatures and immune-suppressive cells. The remaining patients in the immune class were allocated to an immune-activated subclass. The interferon-γ and granzyme B levels were higher in the immune-activated subclass, whereas the transforming growth factor-β1 and programmed cell death-1 (PD-1) levels were higher in the immune-suppressed subclass. The established molecular classification system was recapitulated in validation cohorts. The immune-activated subclass was predicted to have a better response to anti-PD-1 immunotherapy. The immune-related subclasses were associated with differences in copy number alterations, tumor mutation burden, neoantigens, tumor-infiltrating lymphocyte enrichment, PD-1/programmed death-ligand 1 expression, mutation landscape, and various infiltration immunocytes. Overall, we established a novel immune-related molecular classification of breast cancer, which may be used to select candidate patients for immunotherapy.