Perturbation effects on factors (stimulated cells)

Fisrt of all, we look at the estimated effects of gene perturbations on factors inferred by GSFA.

We found that targeting of 9 genes, ARID1A, CBLB, CD5, CDKN1B, DGKA, LCP2, RASA2, SOCS1, and TCEB2, has significant effects (PIP > 0.95) on at least 1 of the 20 inferred factors.

Estimated effects of perturbations on factors (Figure S4A):

dotplot_beta_PIP(t(gibbs_PM$Gamma1_pm), t(gibbs_PM$beta1_pm),
                 marker_names = KO_names,
                 reorder_markers = c(KO_names[KO_names!="NonTarget"], "NonTarget"),
                 inverse_factors = F) +
  coord_flip()

Past versions of Fig_S4A-1.png

Here is a closer look at the estimated effects of selected perturbations on selected factors (Figure 3A):

targets <- c("ARID1A", "LCP2", "CD5", "CBLB", "RASA2", 
             "DGKA", "TCEB2", "SOCS1", "CDKN1B")
complexplot_perturbation_factor(gibbs_PM$Gamma1_pm[-nrow(gibbs_PM$Gamma1_pm), ],
                                gibbs_PM$beta1_pm[-nrow(gibbs_PM$beta1_pm), ],
                                marker_names = KO_names, reorder_markers = targets,
                                reorder_factors = c(2, 4, 9, 12))

Past versions of Fig_3A-1.png

We can also assess the correlations between each pair of perturbation and inferred factor.
The distribution of correlation p values show significant signals in stimulated cells.

## Indices of stimulated cells:
stim_cells <-
  (1:nrow(G_mat))[startsWith(rownames(G_mat), "D1S") | 
                    startsWith(rownames(G_mat), "D2S")]
gibbs_res_tb <- make_gibbs_res_tb(gibbs_PM, G_mat, compute_pve = F,
                                  cell_indx = stim_cells)
heatmap_matrix <- gibbs_res_tb %>% select(starts_with("pval"))
rownames(heatmap_matrix) <- 1:nrow(heatmap_matrix)
colnames(heatmap_matrix) <- colnames(G_mat)

summ_pvalues(unlist(heatmap_matrix),
             title_text = "GSFA Stimulated\n(21 Targets x 20 Factors)")

Past versions of unnamed-chunk-5-1.png

Perturbation effects on factors (unstimulated cells)

In unstimulated cells, only three pairs of associations were detected at PIP > 0.95, which is unsurprising given the role of these targeted genes in regulating T cell responses (Figure S4B):

dotplot_beta_PIP(t(gibbs_PM$Gamma0_pm), t(gibbs_PM$beta0_pm),
                 marker_names = KO_names,
                 reorder_markers = c(KO_names[KO_names!="NonTarget"], "NonTarget"),
                 inverse_factors = F) +
  coord_flip()

Past versions of Fig_S4B-1.png

Correlation within factors

Since the GSFA model does not enforce orthogonality among factors, we first inspect the pairwise correlation within them to see if there is any redundancy. As we can see below, the inferred factors are mostly independent of each other.

plot_pairwise.corr_heatmap(input_mat_1 = gibbs_PM$Z_pm,
                           corr_type = "pearson",
                           name_1 = "Pairwise correlation within factors (Z)",
                           label_size = 10)

Past versions of unnamed-chunk-6-1.png

plot_pairwise.corr_heatmap(input_mat_1 = (gibbs_PM$F_pm > 0.95) * 1,
                           corr_type = "jaccard",
                           name_1 = "Pairwise correlation within \nbinarized gene loadings (F_pm > 0.95)",
                           label_size = 10)

Past versions of unnamed-chunk-7-1.png

Gene loading in factors

To understand these latent factors, we inspect the loadings (weights) of several marker genes for T cell activation or proliferation states in them.

We visualize both the gene PIPs (dot size) and gene weights (dot color) in all factors (Figure S4C):

complexplot_gene_factor(genes_df, interest_df, gibbs_PM$F_pm, gibbs_PM$W_pm)

Past versions of Fig_S4C-1.png

A closer look at some factors that are associated with perturbations (Figure 3C):

complexplot_gene_factor(genes_df, interest_df, gibbs_PM$F_pm, gibbs_PM$W_pm,
                        reorder_factors = c(2, 4, 9, 12))

Past versions of Fig_3C-1.png

GO enrichment analysis in factors

To further characterize these latent factors, we perform GO (gene ontology) enrichment analysis of genes loaded on the factors using WebGestalt.

Foreground genes: Genes w/ non-zero loadings in each factor (gene PIP > 0.95);
Background genes: all 6000 genes used in GSFA;
Statistical test: hypergeometric test (over-representation test);
Gene sets: GO Slim "Biological Process" (non-redundant).

## The "WebGestaltR" tool needs Internet connection.
enrich_db <- "geneontology_Biological_Process_noRedundant"
PIP_mat <- gibbs_PM$F_pm
enrich_res_by_factor <- list()
for (i in 1:ncol(PIP_mat)){
  enrich_res_by_factor[[i]] <- 
    WebGestaltR::WebGestaltR(enrichMethod = "ORA",
                             organism = "hsapiens",
                             enrichDatabase = enrich_db,
                             interestGene = genes_df[PIP_mat[, i] > 0.95, ]$ID,
                             interestGeneType = "ensembl_gene_id",
                             referenceGene = genes_df$ID,
                             referenceGeneType = "ensembl_gene_id",
                             isOutput = F)
}

Several GO “biological process” terms related to immune responses or cell cycle are enriched in factors 2, 4, 9, and 12 (Figure 3D):

factor_indx <- 2
terms_of_interest <- c("kinetochore organization", "chromosome segregation",
                       "cell cycle G2/M phase transition", "cytokinesis")
barplot_top_enrich_terms(enrich_res_by_factor[[factor_indx]],
                         terms_of_interest = terms_of_interest,
                         str_wrap_length = 50) +
  labs(title = paste0("Factor ", factor_indx))

Past versions of Fig_3D_factor2-1.png

factor_indx <- 9
terms_of_interest <- c("microtubule cytoskeleton organization involved in mitosis",
                       "chromosome segregation", "cytokinesis", "cell cycle checkpoint")
barplot_top_enrich_terms(enrich_res_by_factor[[factor_indx]],
                         terms_of_interest = terms_of_interest,
                         str_wrap_length = 35) +
  labs(title = paste0("Factor ", factor_indx))

Past versions of Fig_3D_factor9-1.png

factor_indx <- 4
terms_of_interest <- c("response to chemokine", "cell killing", "leukocyte migration",
                       "response to interferon-gamma", "cytokine secretion")
barplot_top_enrich_terms(enrich_res_by_factor[[factor_indx]],
                         terms_of_interest = terms_of_interest,
                         str_wrap_length = 35) +
  labs(title = paste0("Factor ", factor_indx))

Past versions of Fig_3D_factor4-1.png

factor_indx <- 12
terms_of_interest <- c("leukocyte cell-cell adhesion", "extrinsic apoptotic signaling pathway",
                       "cell killing", "T cell activation", "NIK/NF-kappaB signaling")
barplot_top_enrich_terms(enrich_res_by_factor[[factor_indx]],
                         terms_of_interest = terms_of_interest,
                         str_wrap_length = 35) +
  labs(title = paste0("Factor ", factor_indx))

Past versions of Fig_3D_factor12-1.png