Gene regulatory network inference from CRISPR perturbations in primary CD4+ T cells elucidates the genomic basis of immune disease

Abstract

The effects of genetic variation on complex traits act mainly through changes in gene regulation. Although many genetic variants have been linked to target genes in cis, the trans-regulatory cascade mediating their effects remains largely uncharacterized. Mapping trans-regulators based on natural genetic variation has been challenging due to small effects, but experimental perturbations offer a complementary approach. Using CRISPR, we knocked out 84 genes in primary CD4⁺ T cells, targeting inborn error of immunity (IEI) disease transcription factors (TFs) and TFs without immune disease association. We developed a novel gene network inference method called linear latent causal Bayes (LLCB) to estimate the network from perturbation data and observed 211 regulatory connections between genes. We characterized programs affected by the TFs, which we associated with immune genome-wide association study (GWAS) genes, finding that JAK-STAT family members are regulated by KMT2A, an epigenetic regulator. These analyses reveal the trans-regulatory cascades linking GWAS genes to signaling pathways.

Keywords: CRISPR; GWAS; RNA-seq; T cells; gene regulatory networks; immunology; inbornb errors of immunity; network inference.

Graphical abstract

Figure 1

Study overview Schematic describing the three gene sets that were perturbed with CRISPR knockouts and modeling of the gene network, network inference analyses, and gene module identification and integration with immune GWAS data.

Figure 2

The gene network of the 84 perturbed genes (A) Estimate of the directed network that describes how the 84 perturbed genes interact. The radius of each point is proportional to the degree of that gene. Arrows are used to indicate directionality of the edges, such that an arrow pointing into a gene indicates that it is being regulated by another gene. For emphasis, the opacity of the edges from or to inborn error of immunity transcription factors is increased, and all other edges are displayed with greater transparency. Positive values in the color scale indicate that the parent gene is a positive regulator of the child gene. (B) A sub-network centered around STAT1. (C) A scatterplot of the indegree and outdegree of each of the 84 genes. (D) Association analyses between gene properties and their indegree, outdegree, and total degree.

Figure 3

The landscape of downstream effects (A) The statistical model used to relate the 84 perturbed genes to the expressed genes. (B) The distribution of the number of downstream effects for each of the 84 genes, stratified by gene group. Genes that are outliers with respect to their gene group distribution are labeled. (C) The distribution of indegree for each of the non-perturbed genes. Outlier genes are labeled. (D) Association between the properties of downstream genes and the gene set of the upstream regulators. Coefficients are estimated with negative binomial regressions of the gene-set-specific indegree. Downstream gene annotations are indicated on the y axis, and the facets are used to indicate the gene set of the upstream regulator. ∗The 306 denotes the total number of RNA-seq observations, which includes 84 genes perturbed in three donors and 54 samples from control guides.

Figure 4

The discovery of gene modules Hierarchical clustering is used to identify clusters of shared downstream effects. The upstream gene members within each module are labeled in the left-handed margin of the plot, and the gene group of each gene is indicated by the text color. The total number of genes in the module, including both upstream and downstream effects, is included under the list of genes.

Figure 5

Gene module characterization (A) Enrichment analyses of KEGG genetic, immune, and signaling pathways for each of the 84 perturbed genes, stratified by gene module. The JAK-STAT pathway is highlighted with a dashed red box. The color bar maximum is set to 4. (B) The JAK-STAT sub-network, which is organized such that cytokine genes are at the bottom and upstream regulators are at the top. (C) Effects of KOs in the gene modules on a proliferation assay. Each point represents an individual gene perturbation sample plotted as the log₂ fold change sample count as compared to AAVS1 KO control samples from the same donor (∗p < 0.05 and ∗∗∗∗p < 0.001; n = 3 donors per KO, the number of KOs per cluster is reflected in Figure 4).

Figure 6

Autoimmune and allergy association of SNPs linked with the gene modules (A) Estimated $\tau$ coefficients from LD score regression are plotted for each gene module and phenotype. Module 0 is defined as genes that were not included in any module but are still expressed in CD4⁺ T cells. (B) Exemplar analysis annotating the fine-mapped genes from a Finngen dermatitis GWAS based on their presence in module 2A.

Figure 7

The transcriptional logic linking module 4 to GWAS loci (A) The sub-network of module 4 and Th17 cytokines. (B) Locus plot including tracks describing the functional characteristics of the region. Each track is constructed from publicly available ChIP-seq data (STAR Methods) or ATAC-seq data from Freimer et al. Gray boxes indicate significantly different regions between the respective KO and AAVS1 control KO ATAC data (adjusted p value [p_adj] < 0.05, n = 3 donors per KO). The y axis displays normalized counts.