4F. Construct a Gene Regulatory Network (GRN)

• The focus is on what is regulating the transformation process in smooth muscle cells and cells of the renin lineage
• Related to Figure 2.

renin_obj <- trj_object@data[, rownames(trj_object@data@meta.data[trj_object@data$cell_type_manual %in% c("SMC", "tSMC_1", "tSMC_2", "JG"), ])]
regions <- Links(trj_object)
renin_obj <- initiate_grn(renin_obj,
                          regions = regions,
                          peak_assay = "ATAC",
                          rna_assay = "RNA",
                          exclude_exons=TRUE)
DefaultAssay(renin_obj@data) <- "RNA"
Idents(renin_obj@data) <- renin_obj@data$condition
DEGs <- FindMarkers(renin_obj@data, features=gene_names,
                    ident.1 = "trt", ident.2 = "ctrl", test.use="MAST")
sig_gene_names <- rownames(DEGs[DEGs$p_val_adj < 0.01,])
# Must also remove any genes beginning with a number for the infer_grn() fcn 
sig_gene_names <- sort(sig_gene_names[!grepl("^[0-9](.*)", sig_gene_names)])
motif_dt <- readRDS("motif_datatable.rds")
seqlevelsStyle(BSgenome.Mmusculus.UCSC.mm39) <- "UCSC"
# This calculates the motif presence
renin_obj <- find_motifs(
    renin_obj,
    pfm = motifs,
    genome = BSgenome.Mmusculus.UCSC.mm39,
    motif_tfs = motif_dt
)
# This calculates the motif enrichment
DefaultAssay(renin_obj@data) <- "ATAC"
seqlevelsStyle(BSgenome.Mmusculus.UCSC.mm39) <- "NCBI"
renin_obj@data <- RunChromVAR(renin_obj@data,
    genome = BSgenome.Mmusculus.UCSC.mm39)
DefaultAssay(renin_obj) <- 'chromvar'
seqlevelsStyle(BSgenome.Mmusculus.UCSC.mm39) <- "UCSC"
renin_grn <- infer_grn_mouse(
        renin_obj,
        genes = sig_gene_names,
        peak_to_gene_method = 'GREAT',
        method = "glm",
        parallel = T,
        org = "mouse"
    )
saveRDS(renin_grn, paste0("trj-object_renin-lineage_GRN-GREAT-glm_SeuratLinkPeaks-regions_", Sys.Date(), ".rds"))
renin_grn <- find_modules(renin_grn,
                    p_thresh = 0.25,
                    rsq_thresh = 0.1,
                    min_genes_per_module = 5)
grn_coef <- as.data.table(coef(renin_grn))
grn_sig <- grn_coef[padj < 0.1,]
# Highest correlations
grn_sig[order(abs(corr), decreasing=T),]
renin_grn <- get_network_graph(renin_grn,
                               graph_name='umap_graph')
modules <- NetworkModules(renin_grn)
module_metadata <- modules@meta
DEGs_GRN_plt <- plot_network_graph(renin_grn, graph='umap_graph', color_nodes=TRUE)
umap_graph <- NetworkGraph(renin_grn, graph='umap_graph')
# Get NODE data
gene_centrality <- umap_graph %N>% as_tibble()
gene_centrality <- gene_centrality[order(abs(gene_centrality$centrality), decreasing=T),]
# Get EDGE data
gene_estimate <- umap_graph %E>% as_tibble()
gene_estimate <- gene_estimate[order(abs(gene_estimate$estimate), decreasing=T),]
fwrite(gene_centrality, paste0("trj-object_renin-lineage_GRN-GREAT-glm_SeuratLinkPeaks-regions_UMAP-graph_gene-centrality-table_", Sys.Date(), ".csv"), sep=",")
genes_show <- gene_centrality$name[1:100]
current_labels <- DEGs_GRN_plt$data$name
updated_labels <- copy(current_labels)
for (gene in current_labels) {
    if (gene %ni% genes_show) {
        position <- grep(gene, updated_labels)
        updated_labels[position] <- ""
    }
}
DEGs_GRN_plt$data$name <- updated_labels
ggsave(paste0("GRN/trj-object_renin-lineage_GeneModulesQ05-gene-list_SeuratLinkPeaks-regions_GRN-GREAT-glm_", Sys.Date(), ".png"),
           plot = DEGs_GRN_plt, width = 360, height = 360, units = "mm", dpi = 600)
ggsave(paste0("GRN/trj-object_renin-lineage_GeneModulesQ05-gene-list_SeuratLinkPeaks-regions_GRN-GREAT-glm_", Sys.Date(), ".svg"),
           plot = DEGs_GRN_plt, width = 360, height = 360, units = "mm", dpi = 300)