Diverse transcriptional states from combinatorial TF reprogramming.

A. Schematic of Reprogram-Seq lentiviral TF ORF barcoding system compatible with single-cell RNA-Seq detection. B. Schematic of Reprogram-Seq applied to Tabula Muris cell types. For 7 cell types, 10 to 29 candidate TFs were prioritized for combinatorial over-expression studies. C. For each of the 7 cell type batches, MEFs were infected with lentiviral libraries for the candidate TFs, followed by single-cell RNA sequencing and analysis. D. Distribution of the number of cells recovered for distinct TF cocktails. E. UMAP visualization of reprogrammed MEFs and reference cell types. Reference cells are colored by cell type. Reference cell types (excluding “Other primary cells”) are plotted as larger dots to improve visibility. (insert) Distribution of the number of TFs identified in sequenced cells. F. UMAP feature plot, with cells colored by reprogramming batch. G. UMAP feature plot, with reference cells colored by origin tissue. MEF-derived cells are grey. H. UMAP feature plot, with cells numbered and colored by cluster. I. (left) For each cell cluster indicated, shown is the origin of cells in the cluster (control MEF, reprogrammed MEF, primary cell). The “Primary” label on the bottom indicates the most frequent primary cell type in the cluster. (right) For each primary cell type indicated, the pie chart shows the distribution of the primary cells across the clusters. J. UMAP feature plot, where cells receiving the epicardial TF combination Atf3 + Gata6 + Hand2 are colored orange. K. UMAP feature plot, with cells colored by the number of distinct exogenous TFs detected.

TF combinations drive cell-type specific regulatory programs.

A. Schematic to visualize transcriptional phenotypes from combinatorial TF over-expression. Cells induced with identical TF combinations were combined. B. For the 1335 TF combinations analyzed, shown is the distribution of the number of TFs in each combination. C. The heatmap indicates the expression of 3244 variably expressed genes across the 1335 TF combinations. D. MDE embedding of transcriptomes for 1335 TF combinations. Each dot indicates a TF combination, colored by experimental batch. When a given TF combination is recovered by multiple batches, the size of each dot indicates the fraction of cells derived from the more abundant batch. E. Across the 1335 TF combinations, shown is the fraction of cells derived from the most abundant batch. F. MDE embedding, with TF combinations colored by perturbation cluster and annotated with enriched Gene Ontology terms. G. Heatmap illustrating the enrichment of Gene Ontology terms and pathways across the 51 perturbation clusters. H. The MA plot indicates the differentially expressed genes in perturbation cluster 1 (teal). Several genes with cardiac roles are indicated.

TF combinations drive modular gene regulatory networks.

A. Schematic to estimate gene regulatory networks driven by perturbation clusters. Cells in each perturbation cluster are grouped, followed by differential expression analysis. B. For each perturbation cluster, shown is the number of differentially expressed genes induced and repressed. (left) The number of unique TFs in the perturbation cluster. (right) The number of unique TF combinations. C. Across all differentially expressed genes identified in B, shown is the specificity across modules (1 - fraction of occurrence across modules). A specificity score of 0 indicates a gene is found in all modules. D. Expression of Mylpf in cells across the 51 perturbation clusters. E. Gene ontology enrichment for high specificity and low specificity genes. F. The heatmap indicates the overlap fraction of TFs across perturbation clusters. The histogram quantifies Jaccard similarity of TFs across perturbation clusters. G. Heatmap indicating gene expression patterns (bottom) across perturbation clusters. Each column indicates a distinct TF combination (middle), grouped by batch/assay and perturbation cluster (top). H. The heatmap illustrates the expression of 869 differentially expressed genes in cells induced with Gata6 alone or Gata6 + TFX compared to control MEFs. TFX spans 11 different TFs from different TF families (bold). For each Gata6 + TFX pair, the values on the right column indicate Euclidean distance relative to control cells.

Prediction of TF combinations for transcriptional reprogramming.

A. Pseudotime trajectory of epicardial reprogramming. (left) Indicated are control MEFs (blue), TF-induced MEFs (orange), and primary epicardial cells (yellow). (middle left) Feature plot colored by pseudotime. (middle right) Feature plot colored by the expression of epicardial marker gene Krt19. (right) Feature plot colored by cells induced with Atf3 + Gata6 + Hand2. B. Clusters of distinct gene expression patterns across epicardial transcriptional reprogramming. Pseudotime is indicated by the x-axis. C. The heatmap shows the expression of the gene clusters in C across all MEFs induced with epicardial TFs (ordered by pseudotime). D. Distinct TF combinations are ordered by transcriptional reprogramming performance according to pseudotime. Several TF combinations are indicated. Atf3 + Gata6 + Hand2 was previously identified as an epicardial reprogramming cocktail. E. Bulk qPCR validation of putative TFs for transcriptional reprogramming of epicardial cells. F. Immunostaining of tight junction protein ZO1 (a feature of epicardial cells), in MEFs after over-expression of putative TFs for epicardial reprogramming. G. Quantification of ZO1 positive cells in panel F. H. Across all reprogramming batches and targeted primary cells, shown is the clustering of control MEFs (blue), MEF-derived cells (orange), and primary cells (yellow). I. Pseudotime analysis to rank putative TFs for neuroendocrine and stromal cell reprogramming.

Transcription factor cooperativity in transcriptional reprogramming.

A. A linear model to decompose combinatorial TF perturbations from single TF data. The c1δa and c2δb terms model independent contributions of individual TFs A and B, while the c3 x δa x δa term models cooperativity. B. Distribution of R2 across linear models across 801 TF pairs. C. Across all 801 TF pairs, shown is the distribution of c1 and c2 values in the linear model. D. Distribution of the number of cells used in each application of the linear model. E. Visual embedding of TF pairs based on parameters of the linear model. F. The bar chart indicates the number of TF pairs with the following relationship between constituent TFs: redundant, dominant, cooperative, neomorphic, and none. G. Feature plots indicating TF pairs with redundant, dominant, cooperative, and neomorphic relationships between constituent TFs. Example TF pairs are indicated. H. Example TF pairs with redundant (Hey2 + Nkx2.5), dominant (Hand2 + Wt1), and cooperative (Atf3 + Gata6) relationships. For each of the three plots, shown is the transcriptome expression pattern for: cells induced with each individual TF (top 2 rows), cells induced with both TFs (third row), and the prediction from the linear model (bottom row). The model parameters are indicated above each plot. I. TF interaction relationship in the context of epicardial reprogramming. J. For the epicardial reprogramming batch, cells were ranked by pseudotime and colored by whether they contain the indicated TF pairs. K. qPCR analysis of Gpm6a and Upk1b expression for pairwise and individual TFs.

A. Schematic of approach to model transcriptional responses from single-TF experiments, and to use this model to predict TF perturbations from multi-TF experiments. B. (bottom) Heatmap illustrating the performance of a linear model to predict TF perturbations from transcriptional responses in cells with exactly 1 TF. Columns indicate predictions and rows indicate observations. Heatmap units are the fraction of cells in each row (which sums to 1). (top) TF predictions for cells with Tbx2. C. UMAP embedding of cells with exactly 1 TF. Feature plots indicate TF perturbations. D. Receiver operating characteristic (ROC)-like curves illustrating the performance of predicting single TF perturbations from single-cell transcriptomes, after training on cells perturbed with 1 TF. Dotted line indicates performance where TF labels are randomly shuffled. E. (bottom) Heatmap illustrating the performance of a linear model to predict TF perturbations from transcriptional responses in cells with exactly 2 TFs. (top) F. As in D, but illustrating the performance of predicting one of TFs of the double TF perturbations from single-cell transcriptomes, after training on cells perturbed with 1 TF. G. As in D, but illustrating the performance of predicting two of the TFs of double TF perturbations from single-cell transcriptomes, after training on cells perturbed with 1 TF. H. When predicting the perturbations of cells with 2 TFs (a and b), shown is the frequency that the top prediction is either a or b. I. As in (E), but predicting perturbations in cells that have exactly 3 TFs. J. As in (E), but predicting perturbations in cells that cluster near primary mouse cell types.

A. Schematic of iterative strategy to normalize lentiviral TF titers. As the lentiviral titer of each TF may be different, normalization of relative titers will be important for downstream analysis. To assess the relative titers of TFs in each reprogramming batch, we will create an equimolar pool of plasmids, package pooled lentivirus, infect MEFs for 7 days, and measure the titer for each TF by next-generation sequencing of TF barcodes from mRNA. We then adjust the concentration of each plasmid and iterate this analysis to confirm normalization of relative titers. B. Controlling perturbation complexity by controlling titer. We infected MEFs with a pool of 30 normalized TFs, quantified GFP+ cells by flow cytometry as a measure of absolute titer, and confirmed these measurements by scRNA-Seq. C. The median number of TFs in each batch of perturbation assay. D. Visualization of MEFs demultiplexing. Four batches of MEFs were labeled and pooled together prior to sequencing (cell hashing). The clustering of visualization was performed using the abundance of hashtag oligos. Colors indicate cell groups after demultiplexing. Grey dots indicate multiplets. E. Expression heatmap of infected fluorescent proteins of ‘D5; 3F’ MEFs (infected with 3 fluorescent proteins and cultured for 5 days). F. Cell type compositions of each cluster. Only cell type annotations from Tabula Muris were shown (non-grey colored). G. Cell type compositions of each cluster. Only cell type annotations from reprogrammed cells were shown (non-grey colored).

A. MDE visualization of perturbation clusters. From top left to bottom right are perturbation clusters, number of cells for each perturbation cluster, and number of TFs in each perturbation cluster. B. Differentially expressed genes (DEGs) for selected perturbation clusters. Teel indicates DEGs.

A. Across all differentially expressed genes identified in Fig. 3B, shown is the occurrence across modules). 1 indicates a gene is found only in one module. B. The heatmap indicates the correlation of TFs across perturbation clusters. C. Heatmap indicating gene expression patterns (bottom) across perturbation clusters. Each column indicates a distinct TF combination (middle), grouped by batch/assay and perturbation cluster (top). D. Distribution of the TF family annotations for the 105 TFs 88. E. The MA plot indicates the differentially expressed genes between TF pairs Gata6 + Bnc1 and Gata6 + Wt1 (teal). Several genes are indicated.

A. Pseudotime trajectory of epicardial reprogramming. Colors indicate different cellular states. B. The distribution of the numbers of TFs in perturbation cocktails and the average pseudotimes. C. The distribution of pseudotime ranking and the numbers of cells of these perturbation cocktails. D. Across all reprogramming batches and targeted primary cells, shown is the clustering of control MEFs (blue), MEF-derived cells (orange), and primary cells (yellow). E. Pseudotime analysis to rank putative TFs for cardiomyocytes, lung endothelial cells and mesenchymal stem cells reprogramming.

Visual embedding of TF pairs based on parameters of the linear model, the colors indicating R2 values (A), dcor correlation (B), relationship of c1 and c2 (C), equality of contribution (D), TF A dominant (E), TF B dominant (F), mean absolute error (H), magnitude (I) and variance inflation factor (G).

J. The correlation of average pseudotimes and the number of cooperative TF paris in the perturbation cocktails.

A. Prediction accuracy for 66 TFs using XGBoost algorithm with different cutoffs (0.2266). B. Heatmap visualization of the distribution of predictions (the sum of row is 1). C. The distributions of predictions for Erg-203, Erg-205, and Gata6. D. MA-plots of Wt1 and Mafb perturbed cells. Colors indicate differentially expressed genes. E. Number of cells for each single TF perturbation. Color indicates the prediction accuracy of each TFs.