We aimed to investigate the association between gene co-expression modules and responses to neoadjuvant chemotherapy in breast cancer by using a systematic biological approach. The gene expression profiles and clinico-pathological data of 508 (discovery set) and 740 (validation set) patients with breast cancer who received neoadjuvant chemotherapy were analyzed. Weighted gene co-expression network analysis was performed and identified seven co-regulated gene modules. Each module and gene signature were evaluated with logistic regression models for pathological complete response (pCR). The association between modules and pCR in each intrinsic molecular subtype was also investigated. Two transcriptional modules were correlated with tumor grade, estrogen receptor status, progesterone receptor status, and chemotherapy response in breast cancer. One module that constitutes upregulated cell proliferation genes was associated with a high probability for pCR in the whole (odds ratio (OR) = 5.20 and 3.45 in the discovery and validation datasets, respectively), luminal B, and basal-like subtypes. The prognostic potentials of novel genes, such as MELK, and pCR-related genes, such as ESR1 and TOP2A, were identified. The upregulation of another gene co-expression module was associated with weak chemotherapy responses (OR = 0.19 and 0.33 in the discovery and validation datasets, respectively). The novel gene CA12 was identified as a potential prognostic indicator in this module. A systems biology network-based approach may facilitate the discovery of biomarkers for predicting chemotherapy responses in breast cancer and contribute in developing personalized medicines.