Single-cell transcriptomic analyses now frequently involve elaborate study designs including samples from multiple individuals, experimental conditions, perturbations, and batches from complex tissues. Dimensionality reduction is required to facilitate integration, interpretation, and statistical analysis. However, these datasets often include subtly different cellular subpopulations or state transitions, which are poorly described by clustering. We previously reported a Bayesian matrix factorization algorithm called single-cell hierarchical Poisson factorization (scHPF) that identifies gene co-expression patterns directly from single-cell RNA-seq (scRNA-seq) count matrices while accounting for transcript drop-out and noise. Here, we describe consensus scHPF, which analyzes scHPF models from multiple random initializations to identify the most robust gene signatures and automatically determine the number of factors for a given dataset. Consensus scHPF facilitates integration of complex datasets with highly multi-modal posterior distributions, resulting in factors that can be uniformly analyzed across individuals and conditions. To demonstrate the utility of consensus scHPF, we performed a meta-analysis of a large-scale scRNA-seq dataset from drug-treated, human glioma slice cultures generated from surgical specimens across three major cell types, 19 patients, 10 drug treatment conditions, and 52 samples. In addition to recapitulating previously reported cell type-specific drug responses from smaller studies, consensus scHPF identified disparate effects of the topoisomerase poisons etoposide and topotecan that are highly consistent with the distinct roles and expression patterns of their respective protein targets.