Bournonite (CuPbSbS₃) is an earth-abundant mineral with potential thermoelectric applications. This material has a complex crystal structure (space group Pmn2₁ #31) and has previously been measured to exhibit a very low thermal conductivity (κ < 1 W m⁻¹ K⁻¹ at T ≥ 300 K). In this study, we employ high-throughput density functional theory calculations to investigate how the properties of the bournonite crystal structure change with elemental substitutions. Specifically, we compute the stability and electronic properties of 320 structures generated via substitutions {Na–K–Cu–Ag}{Si–Ge–Sn–Pb}{N–P–As–Sb–Bi}{O–S–Se–Te} in the ABCD₃ formula. We perform two types of transport calculations: the BoltzTraP model, which has been extensively tested, and a newer AMSET model that we have developed and which incorporates scattering effects. We discuss the differences in the model results, finding qualitative agreement except in the case of degenerate bands. Based on our calculations, we identify p-type CuPbSbSe₃, CuSnSbSe₃ and CuPbAsSe₃ as potentially promising materials for further investigation. We additionally calculate the defect properties, finding that n-type behavior in bournonite and the selected materials is highly unlikely, and p-type behavior might be enhanced by employing Sb-poor synthesis conditions to prevent the formation of Sb_Pb defects. Finally, we discuss the origins of various trends with chemical substitution, including the possible role of stereochemically active lone pair effects in stabilizing the bournonite structure and the effect of cation and anion selection on the calculated band gap.