Maintenance of genome integrity is paramount to molecular programs in multicellular organisms. Throughout the lifespan, various endogenous and environmental factors pose persistent threats to the genome, which can result in DNA damage. Understanding the functional consequences of DNA damage requires investigating their preferred genomic distributions and influences on gene regulatory programs. However, such analysis is hindered by both the complex cell-type compositions within organs and the high background levels due to the stochasticity of damage formation. To address these challenges, we developed Paired-Damage-seq for joint analysis of oxidative and single-stranded DNA damage with gene expression in single cells. We applied this approach to cultured HeLa cells and the mouse brain as a proof of concept. Our results indicated the associations between damage formation and epigenetic changes. The distribution of oxidative DNA damage hotspots exhibits cell-type-specific patterns; this selective genome vulnerability, in turn, can predict cell types and dysregulated molecular programs that contribute to disease risks.