PMC9608777_01
Preserved over-expressed and under-expressed modules in sarcoidosis. (A) Heatmap demonstrating WGCNA module-trait relationships in the ACCESS cohort. Columns correspond to clinical feature and rows correspond to module eigengenes. Cells depict the corresponding Spearman correlation coefficient and p-value (in parenthesis). Cell color indicates the strength of the correlation from strong negative (darker blue) to strong positive (darker red). Significance was determined by p-value <0.05. Cyan, greenyellow and magenta modules were found to be significantly correlated with disease (sarcoidosis), whereas blue and brown modules were found to be significantly but negatively correlated with disease. There were no significant relationships between module expression and gender, age, or smoking history. (B) Scatterplot demonstrating module size and module preservation between cohorts based on Z-summary statistic. The magenta and blue modules were strongly preserved (Z-summary statistic > 10) between the ACCESS and UCSF cohorts. (C, D) Module eigengenes, a measure of module expression, demonstrated that the magenta module was similarly over-expressed in sarcoidosis cases (C) and the blue module was similarly under-expressed in sarcoidosis (D) between cohorts. NS, non-significant; **p-value < 0.01; ***p-value < 0.001; Mann-Whitney U Test.
fimmu-13-848759-g001.jpg
PMC9608777_01_fimmu-13-848759-g001.jpg
f1
Of the identified modules, 5 were significantly associated with disease status (Spearman's test p-value <0.05) and one trended towards significance (Spearman's test p-value =0.07) ( Figure 1A ). Module preservation between case-control datasets was determined using the Z-summary statistic with scores greater than 10 indicating strong preservation ( Figure 1B ). Compared to findings in the ACCESS cohort, the magenta module was similarly over-expressed and the blue module was similarly under-expressed in the UCSF cohort ( Figures 1C, D ).
Gene ontology biologic processes enriched in the preserved magenta and blue modules. (A) Highlighted processes in the magenta module are associated with apoptosis, cell death, differentiation, phagocytosis, and proliferation. They comprised approximately 9% (53/587) of all significantly enriched processes within the module (Benjamini-Hochberg adjusted [BH-adj] p-value < 0.1). (B) Highlighted processes in the blue module are all associated with the cell cycle. They comprised approximately 24.7% (18/73) of all significantly enriched processes within the module (BH-adj p-value < 0.1). Description of GO terms are available in 
Supplemental Tables S3A
 (magenta) and 
S3B
 (blue). Bubble plot y-axis represents biologic process significance (-Log10 adjusted p-values) and x-axis represents the GO plot Z-score indicative of the likelihood that a biologic process is up (positive value) or downregulated (negative value). The gold solid horizontal line on the bubble plot represents a BH-adj p-value = 0.1 and the gold dashed horizontal line represents a BH-adj p-value = 0.05.
fimmu-13-848759-g002.jpg
PMC9608777_01_fimmu-13-848759-g002.jpg
f2
UCSF cohort were filtered and normalized to identify 12,047 genes commonly expressed within both datasets. Utilizing WGCNA, 12 distinct modules of gene expression were identified within the 28 subjects sampled in the ACCESS cohort ( Supplemental Figure S3B  and  Supplemental Tables S2A, B ). However, most relevantly, functional analysis also revealed significant enrichment of cell death related pathways in the magenta module and decreased expression of cell survival pathways in the blue module ( Figures 2A, B ).
Absolute immune signal deconvolution of peripheral blood mononuclear cells based on the Human Protein Atlas. Violin plots demonstrate relative expression (log2-pTPM) of cell specific genes. (A) The Monocyte Activity Module (MAM, magenta) was characterized by high relative expression of monocyte and dendritic cell genes and low relative expression of lymphocyte genes. (B) The Lymphocyte Survival Module (LSM, blue) was characterized by high relative expression of lymphocyte genes and low relative expression of monocyte and dendritic cell genes. For (A, B) median log2-pTPM values per cell type (B-cells, Dendritic cells [DC], Monocytes, Natural Killer [NK] cells, and T-cells) are shown in the figure legends. Significant differences in relative expression of cell specific genes were identified between specific cell types within MAM and LSM (Kruskal-Wallis test p-value < 0.05). Post hoc analysis with Benjamini-Hochberg (BH) adjustment for multiple comparisons was performed with Dunn's test and BH.adj p-value <0.1 was considered significant. In MAM, all post hoc comparisons assessing differences in relative expression of cell type specific genes were found to be significant with the exception of B-cell and T-cell genes which were comparable (BH.adj p-value =0.33). Similarly, in LSM all post hoc comparisons were significant except for NK-cell versus T-cells (BH.adj p-value =0.37). pTPM = normalized to transcripts per million protein coding genes. The gold dashed horizontal line represents the median log2-pTPM values among all genes in each module.
fimmu-13-848759-g003.jpg
PMC9608777_01_fimmu-13-848759-g003.jpg
f3
UCSF cohort were filtered and normalized to identify 12,047 genes commonly expressed within both datasets. Utilizing WGCNA, 12 distinct modules of gene expression were identified within the 28 subjects sampled in the ACCESS cohort ( Supplemental Figure S3B  and  Supplemental Tables S2A, B ). Gene expression in these modules was relatively increased in comparison to controls with 59 genes in the cyan module, 12 genes in the greenyellow module, and 287 genes in the magenta module determined to be significantly over-expressed in sarcoidosis (BH.adj p-value <0.1) ( Supplemental Figure S4A  and  Supplemental Table S3A ). In contrast, the significant modules inversely correlated with sarcoidosis included 57 of 1,331 genes in the blue module and 17 of 823 genes in the brown module that were significantly under-expressed (BH.adj p-value <0.1) ( Supplemental Figure S4A  and  Supplemental Table S3A ). Even though there was similar modular expression in the UCSF cohort, the proportion of significantly over-expressed genes within the magenta module and under-expressed genes with the blue module was distinct (173/1,379 and 316/1,331 with BH.adj p-value <0.1; respectively) and likely attributable to clinical, demographic, and geographic heterogeneity between the cohorts ( Supplemental Figure S4B  and  Supplemental Table S3B ). The magenta module was therefore renamed the monocyte activity module (MAM) ( Figure 3A ). T-, B-, and NK lymphocyte protein coding genes, and was renamed the lymphocyte survival module (LSM) ( Figure 3B ). Taken together, these findings suggest that immune dysregulation leading to increased expression of MAM and decreased expression of LSM is likely to contribute to characteristic differential proportions of circulating cell subsets in sarcoidosis.
Significance of transcription factors based on differential transcription factor analysis. Dot plots demonstrate the significance of the top 50 transcription factors involved in the "transition" from health to sarcoidosis based on scaled differential transcription factor involvement (DTFI) within the MAM (A) and LSM (B) modules. DTFI was considered significant with a false discovery rate (FDR) p-value <0.1 utilizing the "z-score" method per the Monster algorithm. Gold solid vertical grid line represents a FDR p-value=0.1 and gold dashed vertical grid line represents FDR p-value=0.05. Numbers enclosed within parentheses next to transcription factor names on the y-axis represent scaled DTFI values standardized against the null distribution. Larger DTFI values indicate a greater magnitude of change in transcription factor gene targeting between states.
fimmu-13-848759-g004.jpg
PMC9608777_01_fimmu-13-848759-g004.jpg
f4
Gene expression in these modules was relatively increased in comparison to controls with 59 genes in the cyan module, 12 genes in the greenyellow module, and 287 genes in the magenta module determined to be significantly over-expressed in sarcoidosis (BH.adj p-value <0.1) ( Supplemental Figure S4A  and  Supplemental Table S3A ). In contrast, the significant modules inversely correlated with sarcoidosis included 57 of 1,331 genes in the blue module and 17 of 823 genes in the brown module that were significantly under-expressed (BH.adj p-value <0.1) ( Supplemental Figure S4A  and  Supplemental Table S3A ). Even though there was similar modular expression in the UCSF cohort, the proportion of significantly over-expressed genes within the magenta module and under-expressed genes with the blue module was distinct (173/1,379 and 316/1,331 with BH.adj p-value <0.1; respectively) and likely attributable to clinical, demographic, and geographic heterogeneity between the cohorts ( Supplemental Figure S4B  and  Supplemental Table S3B ). Despite dissimilarities, consensus gene expression was ascertained ( Supplemental Figure S4C ) and 20.4% of significantly over-expressed genes in the magenta module as well as 8.4% of significantly under-expressed genes in the blue module were common to both cohorts ( Supplemental Figures S4D, E ). Examination of DTFI in MAM identified 11 transcription factors (13.58%) with significantly different targeting patterns in sarcoidosis compared to controls (FDR < 0.1) ( Figure 4A ). Within 1,331 total genes, 36/104 (34.62%) transcription factors were determined to have significant changes in targeting patterns (FDR < 0.1) ( Figure 4B ).
Regulatory interactions between transcription factors and target genes in the "Cell Death" and "Cell Cycle" subnetworks. (A) The MAM cell death subnetwork was comprised of a total of 123 nodes and 557 edges. Nodes consisted of 11 transcription factors with significantly altered regulatory activity and 112 target genes. (B) The LSM cell cycle subnetwork was comprised of a total of 45 nodes and 210 edges. Nodes consisted of 36 transcription factors with significantly altered regulatory activity and 9 target genes. For (A, B), all transcription factors demonstrated significant differential transcription factor involvement (DTFI). Transcription factors are shown as hexagonal nodes if differentially over-expressed in MAM or differentially under-expressed in LSM and as rhomboid nodes if stably expressed (triangle shaped nodes in LSM represent 3 stably over-expressed transcription factors). Transcription factor color indicates direct or indirect relation with cell death or cell cycle GO terms. Target genes are depicted as magenta circles in (A) and blue circles in (B). For all nodes in (A, B), size is proportional to absolute log2-fold change based on differential expression analysis. Edges in (A, B) indicate an interaction between a transcription factor and target gene. Edge width represents the difference in weight of interaction between sarcoidosis cases and controls as determined by the Monster algorithm. MAM and LSM subnetwork structure is also depicted in 
Supplemental Tables S5A, B
.
fimmu-13-848759-g005.jpg
PMC9608777_01_fimmu-13-848759-g005.jpg
f5
The resultant MAM "cell death" subnetwork ( Figure 5A  and  Supplemental Table S5A ) was comprised of 11 DTFI and 112 target genes (nodes) with a total of 557 interactions (edges), whereas the LSM "cell cycle" subnetwork ( Figure 5B  and  Supplemental Table S5B ) was comprised of 36 DTFI and 9 target genes with a total of 210 interactions.
Transcription factor centrality in MAM Cell Death and LSM Cell Cycle subnetworks. Dot plots demonstrate significance of transcription factors regulating the MAM Cell Death (A) and LSM Cell Cycle (B) subnetworks based on Laplacian centrality in descending order.
fimmu-13-848759-g006.jpg
PMC9608777_01_fimmu-13-848759-g006.jpg
f6
SPI1 ( Figure 6A ). These transcription factors had the most impactful regulatory effects on differential gene expression within the MAM subnetwork and are likely key regulatory features driving cell death through heightened innate immune activity. Similar analyses performed on the smaller LSM cell cycle subnetwork identified 4 equally relevant transcription factors: ARID1B, CDK2, KMT2A, and CRY1 responsible for driving expression changes associated with an altered lymphocytic cell cycle ( Figure 6B ). Collectively, altered activity by the transcription factors driving gene expression in the MAM and LSM subnetworks suggests that disrupted T-cell proliferation and activity in a pro-inflammatory background underlies peripheral immunity in sarcoidosis.
Decreased lymphocyte abundance in sarcoidosis linked to MAM and LSM transcription factors. (A, B) Box plots demonstrating the distribution of CD14+ monocyte and CD4+ T-cell in ACCESS and UCSF cohorts based on case-control grouping. PBMC subpopulation proportions were derived from in-silico cellular deconvolution utilizing Cibersort and the Immunostates gene signature matrix. Sarcoidosis cases in both cohorts were found to have significantly lower percentages of CD4+ T-cells when compared to controls within the same cohort. Comparable cell proportions were found among cases and controls between cohorts. For each group, corresponding median, 1st quartile, and 3rd quartile values are shown. P-values are shown above brackets for each pairwise comparison: NS, non-significant; *p-value < 0.05, **p-value < 0.01; Mann-Whitney U Test. (C-F) Linear regression identified significant associations between module eigengenes and cell subset proportions (p-value <0.05). CD14+ monocytes were significantly associated with MAM expression in the ACCESS (C) and UCSF (D) cohorts. CD4+ T-cells were significantly associated with LSM expression in the ACCESS (E) and UCSF (F) cohorts.
fimmu-13-848759-g007.jpg
PMC9608777_01_fimmu-13-848759-g007.jpg
f7
From queried PBMC subsets, only significantly increased CD14+ monocyte and significantly decreased CD4+ T-cell proportions were noted in sarcoidosis when compared against controls within both cohorts (MWU test p-value <0.05) ( Figures 7A, B ). T-cells were directly associated with LSM expression in both ACCESS and UCSF cohorts ( Figures 7C-F ). Thereby, these findings suggest that MAM and LSM are not only differentially expressed in the ACCESS and UCSF cohorts, but that their member genes are likely to modulate activity and characteristic compositional differences of circulating immune cells in sarcoidosis. However, the predictive accuracy of the red random gene list was unexpectedly greater than what was observed in the UCSF cohort validation (>65% in the UIC STAR cohort models versus ~50% in the UCSF cohort models) ( Figures 9C-F  and  Supplemental Tables S7C, D ).
TLE3 and CBX8 target gene expression in MAM predicts imputed monocyte and lymphocyte proportions in the UCSF cohort. Receiver operating characteristic (ROC) curves demonstrate the area under the curve (AUC) as a measure of predictive performance for significant partial least squares discriminant analysis (PLS-DA) models in the validation cohort and negative controls derived from non-preserved modules (red and lightgreen) between ACCESS and UCSF cohorts. Predictive accuracy of CD14+ monocyte proportions (A) and CD4+ T-cell proportions (B) was significantly greater with TLE3 and CBX8 target genes in comparison to both negative control models (p-value < 0.05; Venkatraman test for AUC comparison with 10,000 bootstraps).
fimmu-13-848759-g008.jpg
PMC9608777_01_fimmu-13-848759-g008.jpg
f8
In contrast, target genes of TLE3 and CBX8 were identified as significant predictors of both CD14+ monocyte and CD4+ T-cell abundance and significantly outperformed red and lightgreen random gene lists ( Figures 8A, B ).
Imputed and clinical monocyte and lymphocyte proportions are strongly correlated in the UIC STAR cohort and predicted by TLE3 and CBX8 target gene expression in MAM. Significant correlations between imputed and clinical total monocyte percentage of PBMCs (A) and imputed and clinical CD4+ T-cell percentage of PBMCs (B) were identified (Pearson's and Spearman's test p-value <0.05). Concordance between imputed and clinical monocyte and lymphocyte "low" and "high" groups is depicted by colored circles. (C-F) Receiver operating characteristic (ROC) curves demonstrating the area under the curve (AUC) as a measure of predictive performance for significant partial least squares discriminant analysis (PLS-DA) models in the validation cohort and negative controls derived from non-preserved modules (red and lightgreen) between ACCESS and UIC STAR cohorts. Predictive accuracy of imputed CD14+ monocyte proportions (C) and total clinical monocyte proportions (D) was significantly greater with TLE3 and CBX8 target genes when compared to the lightgreen module random gene list (p-value < 0.05; Venkatraman test for AUC comparison with 10,000 bootstraps) but not the red module random gene list. Predictive accuracy of imputed CD4+ T-cell proportions (E) and clinical CD4+ T-cell proportions (F) was significantly greater with TLE3 and CBX8 target genes in comparison to the lightgreen and red module random gene lists (p-value < 0.05; Venkatraman test for AUC comparison with 10,000 bootstraps).
fimmu-13-848759-g009.jpg
PMC9608777_01_fimmu-13-848759-g009.jpg
f9
CD4+ T-cells and total monocytes, respectively) within the benchmarked ranges for Cibersort. Based on these findings, the corresponding median percentage of imputed CD14+ monocytes (16.84%) and total clinical monocytes (29.97%) as well as the corresponding median percentage of imputed (27.87%) and clinical (30.4%) CD4+ T-cells were utilized to designate subjects in the UIC STAR cohort into monocyte "high" and lymphocyte "low" groups. Concordance between imputed and clinical total monocyte (22/30) and imputed and clinical CD4+ T-cell (26/30) grouping was established ( Figures 9A, B ) and provides additional evidence to support the utilization of in-silico deconvolution to determine PBMC cell proportions in sarcoidosis. However, the predictive accuracy of the red random gene list was unexpectedly greater than what was observed in the UCSF cohort validation (>65% in the UIC STAR cohort models versus ~50% in the UCSF cohort models) ( Figures 9C-F  and  Supplemental Tables S7C, D ).