We have studied four complex organic molecules (COMs), methyl formate ($CH_3OCHO$), dimethyl ether ($CH_3OCH_3$), formamide ($NH_2CHO$), and ethyl cyanide ($C_2H_5CN$), towards a large sample of 39 high-mass star-forming regions representing different evolutionary stages, from early to evolved phases. We aim to identify potential correlations between the molecules and to trace their evolutionary sequence through the star formation process. We analysed spectra obtained at 3, 2, and 0.9 mm with the IRAM-30m telescope. We derived the main physical parameters for each species by fitting the molecular lines. We compared them and evaluated their evolution, also taking several other interstellar environments into account. We report detections in 20 sources, revealing a clear dust absorption effect on column densities. Derived abundances are ~$10^{-10}-10^{-7}$ for $CH_3OCHO$ and $CH_3OCH_3$, ~$10^{-12}-10^{-10}$ for $NH_2CHO$, and ~$10^{-11}-10^{-9}$ for $C_2H_5CN$. The abundances of $CH_3OCHO$, $CH_3OCH_3$, and $C_2H_5CN$ are very strongly correlated (r>0.92) across ~4 orders of magnitude. $CH_3OCHO$ and $CH_3OCH_3$ show the strongest correlations in most parameters, and a nearly constant ratio (~1) over a remarkable ~9 orders of magnitude in luminosity for a wide variety of sources: pre-stellar to evolved cores, low- to high-mass objects, shocks, Galactic clouds, and comets. This indicates that COMs chemistry is likely early developed and then preserved through evolved phases. Moreover, the molecular abundances clearly increase with evolution. We consider $CH_3OCHO$ and $CH_3OCH_3$ to be most likely chemically linked: they could e.g. share a common precursor, or be formed one from the other. We propose a general scenario for all COMs, involving a formation in the cold, earliest phases of star formation and a following increasing desorption with the progressive heating of the evolving core.