We combined multi-wavelength data for blazars from the Roma-BZCAT catalog and analyzed hundreds of X-ray spectra. We present the fluxes and Spectral Energy Distributions (SEDs), in 12 frequency bands from radio to $\gamma$-rays, for a final sample of 2214 blazars. Using a model-independent statistical approach, we looked for systematic trends in the SEDs; the most significant trends involved the radio luminosities and X-ray spectral indices of the blazars. We used a Principal Component Analysis (PCA), to determine the basis vectors of the blazar SEDs and, in order to maximize the size of the sample, imputed missing fluxes using the K-nearest neighbors method. Using more than an order of magnitude more data than was available when Fossati et al. (1997, 1998) first reported trends of SED shape with blazar luminosity, we confirmed the anti-correlation between radio luminosity and synchrotron peak frequency, although with greater scatter than was seen in the smaller sample. The same trend can be seen between bolometric luminosity and synchrotron peak frequency. Finally, we used all available blazar data to determine an empirical SED description that depends only on the radio luminosity at 1.4~GHz and the redshift. We verified that this statistically significant relation was not a result of the luminosity-luminosity correlations that are natural in flux-limited samples (i.e., where the correlation is actually caused by the redshift rather than the luminosity).