We study the self-assembly of amphiphilic cubic colloids using molecular dynamics as well as rejection-free kinetic Monte Carlo simulations. We vary both the number and location of the solvophobic faces (patches) on the cubes at several colloid volume fractions and determine the resulting size and shape distributions of the self-assembled aggregates. When the binding energy is comparable to the thermal energy of the system, aggregates typically consist of only few spontaneously associating/dissociating colloids. Increasing the binding energy (or lowering the temperature) leads to the emergence of highly stable aggregates, e.g., small dimers in pure suspensions of one-patch cubes or large (system-spanning) aggregates in suspensions of multipatch colloids. In mixtures of one- and multipatch cubes, the average aggregation number increases with increasing number of solvophobic faces on the multipatch cubes as well with increasing fraction of multipatch cubes. The resulting aggregate shapes range from elongated rods over fractal objects to compact spheres, depending on the number and arrangement of solvophobic patches on the cubic colloids. Our findings establish the complex self-assembly pathways for a class of building blocks that combine both interaction and shape anisotropy, with the potential of forming hierarchically ordered superstructures.