Simultaneous Pathway Activity Inference and Gene Expression Analysis Using RNA Sequencing

Abstract

Reporter gene assays are a venerable tool for studying signaling pathways, but they lack the throughput and complexity necessary to contribute to a systems-level understanding of endogenous signaling networks. We present a parallel reporter assay, transcription factor activity sequencing (TF-seq), built on synthetic DNA enhancer elements, which enables parallel measurements in primary cells of the transcriptome and transcription factor activity from more than 40 signaling pathways. Using TF-seq in Myd88(-/-) macrophages, we captured dynamic pathway activity changes underpinning the global transcriptional changes of the innate immune response. We also applied TF-seq to investigate small molecule mechanisms of action and find a role for NF-κB activation and coordination of the STAT1 response in the macrophage reaction to the anti-inflammatory natural product halofuginone. Simultaneous TF-seq and global gene expression profiling represent an integrative approach for gaining mechanistic insight into pathway activity and transcriptional changes that result from genetic and small molecule perturbations.

Figure 1. Overview of transcription factor activity sequencing (TF-Seq)

(A) TF-Seq's 58 lentiviral reporter vectors are distinguished by the DNA response element (RE) driving Luc2P transcription, the RE sequence tags (RE-tag) and their associated unique molecular identifiers (UMIs). Emerald-GFP (EmGFP) is downstream of the reporter gene with an invariant SV40 promoter and used to monitor the rate of lentiviral transduction. (B) 58 unique REs were created to represent more than 40 of the most widely studied signaling pathways, with a total degenerate vector complexity of 191,724. n = 3 index-tagged PCR amplicons of the pooled TF-seq vector library, data points represent mean ± s.d. (C) To perform TF-seq, an equimolar ratio of all 58 TF-seq vectors are first transfected into HEK293T cells. After 48 h, the conditioned media containing a pool of lentiviral particles representing the TF-seq reporter library is collected and transduced into a cell type of interest. After 72 h, the heterogeneous population of cells transduced with TF-seq is collected and replated in 96-well microtiter plates for stimulation and assay collection. (D) TF-seq libraries are prepared with gene-specific reverse transcriptase priming of the Luc2P reporter gene. Additional 5′ sequence is included in the TF-seq reverse transcriptase primers to add well-tags and a second UMI to approximate single molecule mRNA counting and pool samples from a 96-well microtiter plate for Illumina adapter PCR and multiplexed sequencing. (E) Multiplexing pathway activity inferences using TF-seq correlates with luciferase pathway activity inferences r = 0.68 (P < 4.2 × 10⁻⁷ based on Spearman correlation test). n = 4 wells for luciferase measurements and n = 4 for the TF-seq measurements. See also Figure S1 and Tables S1 and S2.

Figure 2. TF-seq rediscovers Myd88-dependent innate immune signaling

(A) A graphical abstract of the innate immune receptors and their dependence on different signaling adaptor proteins. (B) The temporal TF-seq activity patterns in bone marrow-derived macrophages after stimulation with LPS. Five manually selected activity profiles highlight the major patterns of responses. Each line represents the smoothened activity profile of distinct transcription factors inferred by TF-seq. Solid lines represent wildtype and dashed lines represent Myd88^−/−. (C) A temporal heatmap that summarizes the TF-seq activity patterns for all PAMP stimulations. Each colored bar represents a condition where TF-seq was significantly different from non-treated samples (FWER < 0.05). The color corresponds to log₂ fold change and the number shows ANOVA P value for the full time series. The bar graph immediately below the heatmap shows the number of differentially expressed genes detected for the entire time series using the same ANOVA test. Bone marrow from two Myd88^−/− and two wildtype littermates were pooled during BMDM differentiation. n = 2 wells for each stimulus and time point. See also Figure S2 and Table S3.

Figure 3. TF-seq pathway activity measurements cannot be reproduced by global RNA-seq or ChIP-seq data alone

(A-B) Comparison of the concordance in the statistical testing results of inferred pathway activity based on TF-seq and gene expression derived pathway activity inferences. Two approaches to infer the pathway activity data were used: using transcription factor mRNA level as a proxy of its activity (A) and averaging the expression levels of reported target genes from published ChIP-seq data (B). n = 2 wells for all TF-seq and gene expression data. See Experimental Procedures for detailed descriptions of expression based pathway activity inferences.

Figure 4. Functional association of signaling pathway activity and global gene expression changes using TF-seq

(A) A heatmap of the gene expression changes associated with 765 genes (75%) of the global gene expression changes after PAMP stimulation and clustered by their correlation (cor > 0.4) with TF-seq pathway activity. The DNA binding motif table on the right represents the significantly enriched DNA binding motifs identified within the open chromatin regions of the genes in the pathway activity clusters. (B) Changes in STAT1 binding occupancy in the genes associated with STAT1 activity by TF-seq exhibits the highest enrichment for bona fide direct targets of STAT1. n = 2 wells for all TF-seq and gene expression data. See also Figure S3 and Table S4.

Figure 5. Small molecule mechanism of action discovery using TF-seq

(A) Multidimensional scaling of samples from small molecule screen based on gene expression data. n = 8 wells. (B) Samples plotted using t-SNE coordinates and colored according to TF-seq pathway activity inferences. n = 8 wells. (C) Halofuginone dose-dependent patterns of selected reporters and expressed as fold change against DMSO control samples. n = 2 wells. (D) Genes that were differentially regulated by halofuginone were associated with signaling pathways based on their expression profiles across six dose-response curves (halofuginone and halofuginone + proline in three time points). The panel depicts the distribution of the genes between six pathway reporters. n = 2 wells. See also Figures S4 and S5 and Tables S5-S7.