Abstract: We present a realistic numerical model of sunspot and active region formation through the emergence of flux tubes generated in a solar convective dynamo. The magnetic and velocity fields in a horizontal layer near the top boundary of the solar convective dynamo simulation are used as a time-dependent bottom boundary to drive the radiation magnetohy drodynamic simulations of the emergence of the flux tubes through the upper most layer of the convection zone to the photosphere. The emerging flux tubes rise in the mean convective speed. They interact with the convection and break into small scale magnetic elements that further rise to the photosphere. At the photosphere, several bipolar pairs of sunspots are formed through the coalescence of the small scale magnetic elements. The sunspot pairs in the simulation successfully reproduce the fundamental observed properties of solar active regions, including the more coherent leading spots with a stronger field strength, and the correct tilts of the bipolar pairs. These asymmetries originate from the intrinsic asymmetries in the emerging fields imposed at the bottom boundary. The leading sides of the emerging flux tubes are up against the downdraft lanes of the giant cells and strongly sheared downward. This leads to the stronger field strength of the leading polarity fields. The azimuthal average of the flows around a (leading) sunspot reveals a predominant down flow inside the sunspots and a large-scale horizontal inflow at the depth of about 15 Mm. The inflow pattern becomes an outflow in upper most convection zone in the vicinity of the sunspot.