We present a systematic study of the cosmic variance that existed in the formation of first stars and galaxies. We focus on the cosmic variance induced by the large-scale density and velocity environment engraved at the epoch of recombination. The density environment is predominantly determined by the dark-matter overdensity, and the velocity environment by the dark matter-baryon streaming velocity. Toward this end, we introduce a new cosmological initial condition generator BCCOMICS, which solves the quasi-linear evolution of small-scale perturbations under the large-scale density and streaming-velocity environment and generates the initial condition for dark matter and baryons, as either particles or grid data at a specific redshift. We also describe a scheme to simulate the formation of first galaxies inside density peaks and voids, where a local environment is treated as a separate universe. The resulting cosmic variance in the minihalo number density and the amount of cooling mass are presented as an application. Density peaks become a site for enhanced formation of first galaxies, which compete with the negative effect from the dark matter-baryon streaming velocity on structure formation.