The abundance of tumor-infiltrating lymphocytes has been associated with a favorable prognosis in estrogen receptor-negative breast cancer. However, a high degree of spatial heterogeneity in lymphocytic infiltration is often observed and its clinical implication remains unclear. Here we combine automated histological image processing with methods of spatial statistics used in ecological data analysis to quantify spatial heterogeneity in the distribution patterns of tumor-infiltrating lymphocytes. Hematoxylin and eosin-stained sections from two cohorts of estrogen receptor-negative breast cancer patients (discovery: n=120; validation: n=125) were processed with our automated cell classification algorithm to identify the location of lymphocytes and cancer cells. Subsequently, hotspot analysis (Getis-Ord Gi*) was applied to identify statistically significant hotspots of cancer and immune cells, defined as tumor regions with a significantly high number of cancer cells or immune cells, respectively. We found that the amount of co-localized cancer and immune hotspots weighted by tumor area, rather than number of cancer or immune hotspots, correlates with a better prognosis in estrogen receptor-negative breast cancer in univariate and multivariate analysis. Moreover, co-localization of cancer and immune hotspots further stratified patients with immune cell-rich tumors. Our study demonstrates the importance of quantifying not only the abundance of lymphocytes but also their spatial variation in the tumor specimen for which methods from other disciplines such as spatial statistics can be successfully applied.