omicsDNA: A multilayer network toolkit for grouped omics and Dynamic Network Analysis (and other high‑dimensional) data
Mehran Piran
2025-09-22
\(~\)
1 Introduction
Gene–gene relationships rarely remain constant across development, disease, or cell types. omicsDNA is an R package for building, stabilising, and comparing layered co‑expression networks so that such shifts are straightforward to see and test. Starting from either bulk expression grouped by condition or single‑cell RNA‑seq organised by cell type, the package estimates one network per group (a “layer”), repeats the estimation under balanced resampling, and aggregates the repeats into a consensus edge set. From there, omicsDNA assembles a multilayer object, detects and visualises communities, quantifies layer‑wise structure, and compares overlaps. Functional interpretation is integrated via g:Profiler, including term–gene “concept” networks and term‑only overlap networks. All major outputs—CSVs, RDS files, and plots—are timestamped for reproducibility. In a vignette based on human foetal kidney scRNA‑seq, we contrast progenitor populations, differentiating lineages, and endothelial cells, highlighting shared and layer‑specific modules.
2 Package description
Most transcriptomic analyses ask whether two genes are co‑expressed within a single cohort. In practice, biology is stratified: developmental time points, disease stages, tissues, and cell types often carry distinct regulatory programmes. The more informative question is how co‑expression patterns change across these contexts. Single‑cell RNA‑seq sharpens this question by allowing network estimation within narrow populations, but sparsity and uneven sampling can make edges unstable.
omicsDNA is designed for these realities. It centres on layered networks: you define groups—age bins, case/control, cell types—and the package constructs one gene–gene network per group. Two settings are supported:
- Bulk:
buildAdjacency()takes a gene × sample matrix plus sample metadata and returns one adjacency per group. - Single‑cell:
sc_buildAdjacency()reads an assay/slot from a Seurat object and treats a metadata column (e.g., cell type) as the grouping variable. Correlations are computed across cells within each type (no pseudobulking), so each layer reflects cell‑level variation.
Because small per‑layer sample sizes and dropouts can yield brittle
edges, stability is a recurring theme. omicsDNA implements
balanced resampling, repeatedly re‑estimating each
layer using the same number of samples (or cells) per draw.
consensusEdges() then aggregates repeats, retaining edges
that recur consistently (e.g., in ≥70% of resamples) and summarising
weights by a chosen statistic (mean or median). This filters noise while
preserving structure that replicates.
With stable per‑layer networks in hand, omicsDNA supports multilayer analysis. Adjacency matrices can be converted to edge lists, combined into a multilayer object, and explored with community detection and layer‑wise metrics (degree, density, centrality). Comparisons become natural: which modules recur across layers, which are layer‑specific, and which genes change community membership as you move from progenitors to differentiating cells?
Interpretation matters as much as topology. For each layer, omicsDNA runs g:Profiler and writes explicit term–gene intersections. Two compact visual summaries are provided: a concept network linking enriched terms to their member genes, and a term‑only overlap network in which edges encode shared genes and node labels show intersection sizes. These views make it easy to spot shared themes and layer‑specific pathways.
Reproducibility is built in. The package writes adjacency lists and
consensus results (RDS), metrics and tables (CSV), and plots (PNG/PDF)
with timestamped filenames and a simple run manifest.
Correlations and p‑values are computed with WGCNA when
available (Pearson, Spearman, biweight midcorrelation) or fall back to
Hmisc, so the workflow runs across a wide range of
environments. Gene‑set choice is explicit: pass feature_ids
to fix the universe across layers, use Seurat’s
VariableFeatures(), or select top‑variance genes; in all
cases, matrices are padded to a common gene order to
streamline downstream comparisons.
The vignette accompanying this release analyses human foetal kidney scRNA‑seq. We build networks for nephron progenitors, differentiating lineages, and endothelial cells; stabilise edges with repeated, balanced draws (e.g., 100 cells per type per repeat); and assemble a multilayer object. Community detection and enrichment highlight programmes shared across compartments (e.g., core developmental signalling) as well as lineage‑specific modules. In short, omicsDNA brings a layer‑centric, statistically careful workflow into one place: build per‑group networks, stabilise them, combine and compare layers, and interpret the results—end to end and reproducibly.
\(~\)
3 Run the basic functions on the seurat object
library(omicsDNA)
library(tidyverse)
library(textshape)
library(igraph)
library(multinet)
library(ndtv)
library(org.Hs.eg.db)
library(msigdbr)
library(fgsea)
library(ggplot2)
library(Seurat)
library(WGCNA)
library(clusterProfiler)
library(gprofiler2)\(~\)
so <- readRDS("~/Niche_analysis/Hsa_fetal_kidney/Hsa_fetal_30perc_mit.rds")
so## An object of class Seurat
## 58393 features across 8413 samples within 3 assays
## Active assay: SCT (22655 features, 0 variable features)
## 3 layers present: counts, data, scale.data
## 2 other assays present: RNA, integrated
## 3 dimensional reductions calculated: pca, umap, umap3dtable(so$Clusters)##
## Endothelium Immune N.CS N.DCT
## 710 255 56 39
## N.Distal_EN N.DistTubule_Dev N.LoH N.Medial_EN
## 296 116 133 280
## N.NP_Str N.NPC N.NPC_CC N.NPC_Primed
## 70 931 110 592
## N.NPC_PTA N.NPC_PTA_CC N.PEC N.Pod_Dev
## 566 135 181 292
## N.Pod_Mat N.ProxTubule_Dev N.ProxTubule_Mat N.RV
## 516 260 222 212
## S.IC S.Med S.Mesangial S.NP_Str
## 718 765 148 50
## S.OC_NZ U.Cortical U.Med_Inner U.Med_Outer
## 477 214 6 33
## U.Tip
## 30DimPlot(so , label = T)
\(~\)
The main call, sc_buildAdjacency(), builds gene–gene
adjacency matrices per cell type (“layer”) from a Seurat object
so. The key inputs are the assay/slot to read expression
from (assay = "SCT", slot = "data"), and the
grouping column (group_col = "Clusters") that defines
layers. Pairwise Pearson correlations are computed
within each layer, and edges are kept only where both the
absolute correlation exceeds
corr_threshold = 0.35. and the FDR‑adjusted
p‑value (Benjamini–Hochberg) is ≤
pval_cutoff = 0.05. It would be better the increase that
correlation threshold when analysing bulk data. With
resample = TRUE, the function performs balanced
resampling: in each repeat it draws up to
samples_per_group = 100 cells per layer (using all
available if fewer), computes correlations, and thresholds as above.
You’ve requested n_repeats = 10, so this process runs ten
times. group_order = unique(so$Clusters) fixes the order of
layers in the output; verbose = TRUE prints progress/QC
notes. Results are saved to disk (save_rds = TRUE) in
out_dir = "omicsDNA_results" with the base name
file_prefix = "sc_adjacency" and xz
compression.
The expected output is a list of length 10 (one
element per repeat). Each element is a named sub‑list
with one square adjacency matrix per layer (cell type). Matrices are
genes × genes over a common feature set determined
internally (variable features if present, otherwise top‑variance genes);
the diagonal is zero and non‑zero entries are the
signed correlation coefficients that passed the thresholds. Each matrix
carries an attribute "n_cells" giving the number of cells
used for that layer in that repeat. Layers with too few cells (below the
function’s min_cells_per_layer, default 20 unless you
override it) will be reported in the console and returned as all‑zero
matrices of the correct shape. The top‑level return object includes an
"rds_file" attribute with the full path to the saved RDS,
so you can reload the whole structure later without recomputing.
Pearson method is suggested if scale data is the input.
adj_sc <- sc_buildAdjacency(
seurat_obj = so,
assay = "SCT",
slot = "data",
group_col = "Clusters",
cor_method = "pearson",
pval_adjust = "fdr",
n_var_features = 3000,
pval_cutoff = 0.05,
corr_threshold = 0.35,
resample = TRUE,
samples_per_group = 100,
n_repeats = 10,
min_cells_per_layer = 100,
save_rds = TRUE,
out_dir = "omicsDNA_results",
file_prefix = "sc_adjacency",
compress = "xz",
verbose = TRUE
)\(~\)
# number of iterations
adj_sc %>% length()
adj_sc[[1]] %>% names # layers in the first iteration
adj_sc[[1]]$N.DCT %>% dim()
adj_sc[[1]]$N.DCT[0:10 , 0:10]\(~\)
edgesFromAdjacency() takes adj_sc—your
nested list of per‑cell‑type adjacency matrices produced by
sc_buildAdjacency() across bootstrap repeats—and converts
each matrix into an edge table. With flatten = FALSE, it
preserves the original structure, returning a list of
length n_repeats where each element is a named sub‑list
by layer; each layer element is a data frame with the standard
columns from, to, and weight (one
row per non‑zero matrix entry, i.e., an observed edge and its weight).
The id_cols = c("repeat","group") argument defines how
hierarchy levels would be labeled if you chose a single
flattened data frame (e.g., columns named “repeat” and “group”
for repeat index and layer/cell‑type); since you set
flatten = FALSE, the output remains nested and these names
mainly document the intended identifiers.
# Will save automatically into ./omicsDNA_results/ with a timestamped name
edges_list <- edgesFromAdjacency(
adj_sc,
flatten = FALSE,
id_cols = c("repeat","group")
)length(edges_list)## [1] 10edges_list[[1]]$N.DCT %>% head()## [1] from to weight
## <0 rows> (or 0-length row.names)\(~\)
The function call to consensusEdges() takes the
resampled edge lists from edges_list and distills them into
robust consensus networks. By requiring each edge to be present in at
least 70% of the repeats (prop_present = 0.7) and
summarizing its weight by the mean (summary = "mean"). With
as_list = TRUE, the output cons_list is a list
of per-layer consensus edge tables and statistics on how consistently
the edge appeared across repeats.
cons_list <- consensusEdges(
edges_list,
prop_present = 0.7,
summary = "mean",
as_list = TRUE
)\(~\)
build_multiNet() function builds a
multilayer network from your per‑layer consensus edge lists and a table
of node metadata. The inputs are cons_list
(a named list where each element is a data frame of edges for one layer,
typically with columns from, to, and
weight) and genes_info (a data frame of node
attributes with at least GeneName and
GeneType). The arguments tell the function
how to assemble and annotate the graph:
featureID_col = "GeneName" says “match graph actors to rows
in genes_info using the GeneName column”;
nodeAttrCols = "GeneType" copies that attribute onto each
vertex; directed = FALSE treats edges as undirected; and
aggregate_duplicates = "mean" averages weights when the
same actor pair appears more than once in a layer (recommended for
stable edge sets). save_to_rds = TRUE writes an RDS
snapshot to your results folder, and verbose = TRUE prints
per‑layer summaries (e.g., vertex counts and attribute‑match rates). The
output is a multilayer network object (often shown like
ml-net[...]) whose layers mirror
names(cons_list), whose vertices are the union of actors
referenced in each layer, and whose edges carry your consensus weights;
node attributes (e.g., GeneType) are attached where names
match.
genes_info <- read.csv("~/Abir/Pipeline_Data-n-Codes/countdataGenes_Name.Type.Length_lncRNAs.mRNAs.miRNAs_v4.csv", skip=1)
head(genes_info)## GeneName GeneType Length
## 1 MIR6859-1 miRNA 68
## 2 MIR1302-2HG lncRNA 1556
## 3 MIR1302-2 miRNA 138
## 4 FAM138A lncRNA 1528
## 5 OR4F5 protein_coding 6167
## 6 AL627309.1 lncRNA 44429dim(genes_info)## [1] 38730 3net <- build_multiNet(
edgeListPerLayer = cons_list, # output of consensusEdges(..., as_list=TRUE)
nodesMetadata = genes_info,
featureID_col = "GeneName",
nodeAttrCols = "GeneType",
directed = FALSE,
aggregate_duplicates = "mean", # optional but recommended
save_to_rds = TRUE,
verbose = TRUE
)## Layer summary (unique actors per layer):
## Layer ActorCount
## 1 Endothelium 226
## 2 Immune 1035
## 3 N.Distal_EN 487
## 4 N.DistTubule_Dev 684
## 5 N.LoH 594
## 6 N.Medial_EN 182
## 7 N.NPC 45
## 8 N.NPC_CC 494
## 9 N.NPC_Primed 67
## 10 N.NPC_PTA 343
## 11 N.NPC_PTA_CC 455
## 12 N.PEC 617
## 13 N.Pod_Dev 197
## 14 N.Pod_Mat 81
## 15 N.ProxTubule_Dev 178
## 16 N.ProxTubule_Mat 843
## 17 N.RV 342
## 18 S.IC 136
## 19 S.Med 66
## 20 S.Mesangial 418
## 21 S.OC_NZ 80
## 22 U.Cortical 317net## ml-net[2076, 22, 7887, 51838 (51838,0)]\(~\)
add_network_attributes() enriches the
existing multilayer network net with additional
node and edge attributes by joining
metadata. The inputs are the same
genes_info node table and the per‑layer edge tables
cons_list. The arguments specify how to
join: for nodes,
featureID_col = "GeneName" and
nodeAttrCols = "GeneType" add GeneType to
every vertex matched by gene name; for edges,
edgesMetadata = cons_list with
edge_from = "from" and edge_to = "to"
identifies the endpoint columns, and
edgeAttrCols = c("n_present","n_repeats","prop_present")
copies those consensus metrics onto each edge. To make matching robust,
actor_key_normalize = c("strip_version","trim","tolower")
aligns identifiers (e.g., turns “TPM1.1” into “tpm1”),
map_edges_by_actor_key = TRUE applies that normalization to
edge endpoints too, and auto_detect_keys = FALSE avoids
guessing any other key columns. directed_network = FALSE
keeps edges undirected; verbose = TRUE prints a join
summary; and save_report = TRUE writes a small audit file
(and, if needed for debugging, export_edge_join_csv = FALSE
can be flipped to TRUE). The output is the
same net object returned with enriched vertex attributes
(GeneType) and enriched edge attributes
(n_present, n_repeats,
prop_present) across all layers.
net <- add_network_attributes(
net,
nodesMetadata = genes_info,
featureID_col = "GeneName",
nodeAttrCols = "GeneType",
edgesMetadata = cons_list,
edge_from = "from",
edge_to = "to",
edgeAttrCols = c("n_present","n_repeats","prop_present"),
actor_key_normalize = c("strip_version","trim","tolower"),
auto_detect_keys = FALSE,
map_edges_by_actor_key = TRUE,
directed_network = FALSE,
verbose = TRUE,
save_report = TRUE, # writes a small report
export_edge_join_csv = FALSE # flip to TRUE if debugging joins
)## Attached 3 edge attribute(s) to 51838 edges.
## Attributes present in network:
## name type
## 1 GeneType string\(~\)
Print a report of attributes attached to a multilayer network.
# After you enriched the network with add_network_attributes(...):
rep <- network_attributes(
net,
show_examples = TRUE, # print small tables of actors/edges with attrs
show_max = 12, # up to 12 names/rows per section
verbose = TRUE
)## Actor-attribute attachment report:
## Coverage: 93.4% Matched: 1939/2076 Normalizer: strip_version -> trim -> tolower
## Unmatched examples: VEGFA, WASH4P, SEPT10, KIRREL, LINC00969, FAM63B, HRSP12, UFD1L, ATP5O, H2AFY
## Edge-attribute attachment report:
## Edges in net: 51838 | Metadata rows in: 51838 | Dropped (unmapped): 0
## Edge attrs added (3): n_present, n_repeats, prop_present
##
## === Attribute Report =====================================
##
## [Added (from stored reports)]
## Actor attrs added: <not recorded; showing 'present' below>
## Edge attrs added (3): n_present, n_repeats, prop_present
##
## [Present now]
## Actor attributes (1): GeneType
## name type
## GeneType string
## Edge attributes (4): weight, n_present, prop_present, n_repeats
## layer name type
## N.NPC_Primed weight double
## N.NPC_Primed n_present double
## N.NPC_Primed prop_present double
## N.NPC_Primed n_repeats double
## N.NPC_PTA n_present double
## N.NPC_PTA prop_present double
## N.NPC_PTA weight double
## N.NPC_PTA n_repeats double
## N.ProxTubule_Dev weight double
## N.ProxTubule_Dev n_repeats double
## N.ProxTubule_Dev prop_present double
## N.ProxTubule_Dev n_present double
## Network-level attributes (0): <none>rep$added$edge_attrs # what edge attrs were attached (per report)## [1] "n_present" "n_repeats" "prop_present"rep$present$actor # current actor-attribute definitions## name type
## 1 GeneType stringrep$samples$edges # first rows of edges with attribute columns## NULL\(~\)
detectCom() takes your multilayer network
net and runs a community‐detection algorithm on it. In your
call you choose method = "louvain" (a fast modularity‑based
method) and tell the routine which edge weights to use via
edgeWeight = "count" (i.e., use the edge attribute named
"count" if it exists; otherwise most workflows treat edges
as unweighted). The size/robustness filters keep only communities with
at least min.actors = 15 unique actors and
present in at least min.layers = 2 layers.
With save_to_rds = TRUE and write_csv = FALSE,
the function saves an RDS snapshot (and skips CSVs). The
output (comm) is a tidy data frame
assigning each (actor, layer) to a community (typically
with columns like actor, layer, and either
cid and/or com), filtered by your minimum
size/layer criteria.
comm <- detectCom(
net,
method = "louvain",
edgeWeight = "count",
min.actors = 15,
min.layers = 2,
seed = 1, # optional: reproducibility
save_to_rds = TRUE,
write_csv = FALSE, # set TRUE if you want a CSV too
verbose = TRUE
)\(~\)
plotCom() visualizes those communities on the underlying
network. You pass the same net and the assignment table
comm. The layout is set to
layout = "multiforce" (a force‑directed layout adapted for
multilayer graphs), with node sized via vertex.size = 5 and
labels hidden (show.labels = FALSE);
gravity = 0.3 control the force layout.
plotCom(
net, comm,
layout = "multiforce",
vertex.size = 1,
show.labels = FALSE,
gravity = 0.3,
seed = 1,
save_plot = TRUE,
format = "png",
width = 10, height = 8, units = "in", dpi = 300,
# key bits for your requests:
show_in_rstudio = TRUE, # show in Plots pane so you can Export
save_df = TRUE, # save the long actor-layer-cid table
df_format = "rds" # or c("rds","csv") if you want both
)
\(~\)
com_vs_samples() relates how many
communities each layer contains to how many
samples/cells you had in that layer.
Idents(so) <- "Clusters"
n <- multinet::layers_ml(net)
n## [1] "N.NPC_Primed" "N.NPC_PTA" "N.ProxTubule_Dev" "N.ProxTubule_Mat"
## [5] "S.IC" "S.Med" "S.Mesangial" "N.DistTubule_Dev"
## [9] "N.NPC_CC" "S.OC_NZ" "U.Cortical" "N.RV"
## [13] "N.Medial_EN" "Immune" "N.LoH" "Endothelium"
## [17] "N.Distal_EN" "N.NPC_PTA_CC" "N.PEC" "N.Pod_Dev"
## [21] "N.Pod_Mat" "N.NPC"num_samples <- subset(so , idents = n)$Clusters %>% table()
num_samples## .
## Endothelium Immune N.Distal_EN N.DistTubule_Dev
## 710 255 296 116
## N.LoH N.Medial_EN N.NPC N.NPC_CC
## 133 280 931 110
## N.NPC_Primed N.NPC_PTA N.NPC_PTA_CC N.PEC
## 592 566 135 181
## N.Pod_Dev N.Pod_Mat N.ProxTubule_Dev N.ProxTubule_Mat
## 292 516 260 222
## N.RV S.IC S.Med S.Mesangial
## 212 718 765 148
## S.OC_NZ U.Cortical
## 477 214comm$layer %>% table()## .
## Endothelium Immune N.Distal_EN N.DistTubule_Dev
## 226 1034 487 657
## N.LoH N.Medial_EN N.NPC N.NPC_CC
## 584 182 45 470
## N.NPC_Primed N.NPC_PTA N.NPC_PTA_CC N.PEC
## 67 343 453 615
## N.Pod_Dev N.Pod_Mat N.ProxTubule_Dev N.ProxTubule_Mat
## 197 81 178 835
## N.RV S.IC S.Med S.Mesangial
## 342 136 66 414
## S.OC_NZ U.Cortical
## 80 317comm$layer <- trimws(as.character(comm$layer))
num_samples <- setNames(as.integer(num_samples), trimws(names(num_samples)))
keep <- intersect(unique(comm$layer), names(num_samples))
num_samples_comm <- num_samples[keep]
layer_cid_counts <- com_vs_samples(comm, num_samples_comm)## Endothelium Immune N.DistTubule_Dev N.Distal_EN
## 6 7 7 7
## N.LoH N.Medial_EN N.NPC N.NPC_CC
## 7 6 5 7
## N.NPC_PTA N.NPC_PTA_CC N.NPC_Primed N.PEC
## 7 7 6 7
## N.Pod_Dev N.Pod_Mat N.ProxTubule_Dev N.ProxTubule_Mat
## 6 6 7 7
## N.RV S.IC S.Med S.Mesangial
## 7 6 4 7
## S.OC_NZ U.Cortical
## 6 7
\(~\)
analyze_actor_overlap(net, reorder = TRUE) and
analyze_edge_overlap(net, reorder = TRUE) are complementary
functions for exploring how layers in a multilayer network relate to one
another. The actor overlap function measures the similarity of node sets
(e.g., genes) across layers, producing a Jaccard similarity matrix
visualized as an “Actors overlap between layers” heatmap, while
the edge overlap function compares the connectivity patterns of those
nodes, generating a Jaccard matrix shown as an “Edge overlap between
layers” heatmap. In both cases, the argument
reorder = TRUE reorders layers by similarity.
actor_jaccard_matrix <- analyze_actor_overlap(net, reorder = TRUE)
\(~\)
edge_jaccard_matrix <- analyze_edge_overlap(net, reorder = TRUE)
\(~\)
\(~\)
the next function takes a multilayer network as input and computes a comprehensive suite of graph statistics at both the global and node levels. The function returns: (i) summary, a per-layer table reporting global metrics such as the number of nodes and edges, density, number of connected components, average path length, diameter, and degree distribution statistics; (ii) nodes, a named list containing node-level centrality measures for each layer, including degree, betweenness, eigenvector centrality, closeness, PageRank, clustering coefficient, and coreness.
res <- layer_metrics(net)
# You should now see non-zero nodes/edges for the layers that have edges
res$summary[, c("layer","n_nodes","n_edges","density")]## layer n_nodes n_edges density
## Endothelium Endothelium 226 570 0.02241888
## Immune Immune 1035 17225 0.03219055
## N.Distal_EN N.Distal_EN 487 5219 0.04410137
## N.DistTubule_Dev N.DistTubule_Dev 684 1864 0.00797993
## N.LoH N.LoH 594 2593 0.01472283
## N.Medial_EN N.Medial_EN 182 585 0.03551697
## N.NPC N.NPC 45 100 0.10101010
## N.NPC_CC N.NPC_CC 494 1443 0.01185011
## N.NPC_Primed N.NPC_Primed 67 148 0.06693804
## N.NPC_PTA N.NPC_PTA 343 1829 0.03118340
## N.NPC_PTA_CC N.NPC_PTA_CC 455 2395 0.02318827
## N.PEC N.PEC 617 6138 0.03229914
## N.Pod_Dev N.Pod_Dev 197 819 0.04242204
## N.Pod_Mat N.Pod_Mat 81 106 0.03271605
## N.ProxTubule_Dev N.ProxTubule_Dev 178 651 0.04132546
## N.ProxTubule_Mat N.ProxTubule_Mat 843 5716 0.01610581
## N.RV N.RV 342 1524 0.02613572
## S.IC S.IC 136 243 0.02647059
## S.Med S.Med 66 123 0.05734266
## S.Mesangial S.Mesangial 418 1242 0.01425080
## S.OC_NZ S.OC_NZ 80 151 0.04778481
## U.Cortical U.Cortical 317 1154 0.02304037# Inspect one layer’s node metrics
head(res$nodes[["N.LoH"]])## actor degree betweenness eigenvector closeness pagerank clustering
## 1 EDNRB 45 0.0203628784 0.0015568126 0.0005005005 0.005908063 0.3141414
## 2 ITGA8 13 0.0007771128 0.0003786623 0.0004325260 0.001921354 0.6153846
## 3 KIF26B 24 0.0020400543 0.0013878918 0.0004849661 0.003128291 0.5579710
## 4 MEIS1 25 0.0024930142 0.0013963593 0.0004849661 0.003240313 0.5133333
## 5 MECOM 9 0.0101068595 0.0131338903 0.0005252101 0.002127543 0.3333333
## 6 CKB 9 0.0068917275 0.0056445326 0.0005035247 0.001857488 0.2500000
## coreness
## 1 14
## 2 9
## 3 12
## 4 12
## 5 5
## 6 6# inspect
res$summary # per-layer summary stats## layer n_nodes n_edges density n_components
## Endothelium Endothelium 226 570 0.02241888 14
## Immune Immune 1035 17225 0.03219055 6
## N.Distal_EN N.Distal_EN 487 5219 0.04410137 8
## N.DistTubule_Dev N.DistTubule_Dev 684 1864 0.00797993 71
## N.LoH N.LoH 594 2593 0.01472283 40
## N.Medial_EN N.Medial_EN 182 585 0.03551697 5
## N.NPC N.NPC 45 100 0.10101010 5
## N.NPC_CC N.NPC_CC 494 1443 0.01185011 60
## N.NPC_Primed N.NPC_Primed 67 148 0.06693804 6
## N.NPC_PTA N.NPC_PTA 343 1829 0.03118340 5
## N.NPC_PTA_CC N.NPC_PTA_CC 455 2395 0.02318827 31
## N.PEC N.PEC 617 6138 0.03229914 20
## N.Pod_Dev N.Pod_Dev 197 819 0.04242204 10
## N.Pod_Mat N.Pod_Mat 81 106 0.03271605 15
## N.ProxTubule_Dev N.ProxTubule_Dev 178 651 0.04132546 9
## N.ProxTubule_Mat N.ProxTubule_Mat 843 5716 0.01610581 30
## N.RV N.RV 342 1524 0.02613572 7
## S.IC S.IC 136 243 0.02647059 6
## S.Med S.Med 66 123 0.05734266 7
## S.Mesangial S.Mesangial 418 1242 0.01425080 20
## S.OC_NZ S.OC_NZ 80 151 0.04778481 8
## U.Cortical U.Cortical 317 1154 0.02304037 16
## avg_path_len diameter degree_mean degree_sd degree_median
## Endothelium 3.915489 12 5.044248 4.614738 3
## Immune 3.452593 9 33.285024 38.875924 16
## N.Distal_EN 2.734260 9 21.433265 31.666133 9
## N.DistTubule_Dev 4.900307 13 5.450292 8.513989 2
## N.LoH 4.325151 12 8.730640 12.622346 3
## N.Medial_EN 3.703986 9 6.428571 6.541230 4
## N.NPC 2.136139 5 4.444444 3.834980 4
## N.NPC_CC 3.632328 10 5.842105 8.939699 2
## N.NPC_Primed 2.336222 7 4.417910 3.080482 4
## N.NPC_PTA 3.342803 8 10.664723 13.714981 4
## N.NPC_PTA_CC 3.690696 11 10.527473 17.992987 3
## N.PEC 3.026828 9 19.896272 31.749562 5
## N.Pod_Dev 2.966959 7 8.314721 11.627397 3
## N.Pod_Mat 2.244318 6 2.617284 2.266980 2
## N.ProxTubule_Dev 3.229630 8 7.314607 8.405131 4
## N.ProxTubule_Mat 3.148749 11 13.561091 26.106206 4
## N.RV 3.588752 10 8.912281 11.356057 4
## S.IC 5.076053 12 3.573529 2.987957 3
## S.Med 2.587145 5 3.727273 3.145426 3
## S.Mesangial 4.424155 11 5.942584 7.336813 3
## S.OC_NZ 2.729645 6 3.775000 3.035549 3
## U.Cortical 4.919281 14 7.280757 9.128827 3
## degree_min degree_max degree_cv degree_skew bet_mean
## Endothelium 1 22 0.9148517 1.4304559 0.003367924
## Immune 1 169 1.1679704 1.4046870 0.002328586
## N.Distal_EN 1 181 1.4774293 2.4493447 0.003330651
## N.DistTubule_Dev 1 75 1.5621160 3.4208072 0.003357646
## N.LoH 1 71 1.4457527 2.4926735 0.004124976
## N.Medial_EN 1 33 1.0175246 1.6341091 0.012814104
## N.NPC 1 18 0.8628705 1.7477164 0.010782241
## N.NPC_CC 1 59 1.5302188 2.6169254 0.002050630
## N.NPC_Primed 1 13 0.6972713 0.7712116 0.005364784
## N.NPC_PTA 1 87 1.2860138 2.1987257 0.006399035
## N.NPC_PTA_CC 1 111 1.7091460 3.1282141 0.004298272
## N.PEC 1 195 1.5957543 2.4661451 0.002862393
## N.Pod_Dev 1 61 1.3984110 2.4643629 0.006847081
## N.Pod_Mat 1 9 0.8661574 1.4887672 0.001711205
## N.ProxTubule_Dev 1 46 1.1490886 1.9919376 0.009987997
## N.ProxTubule_Mat 1 215 1.9250815 3.7439257 0.002131922
## N.RV 1 65 1.2742032 2.2716234 0.006835356
## S.IC 1 15 0.8361361 1.5196208 0.011719995
## S.Med 1 15 0.8438948 1.4739527 0.009353147
## S.Mesangial 1 41 1.2346167 2.2294122 0.004946071
## S.OC_NZ 1 16 0.8041189 1.4039184 0.006722655
## U.Cortical 1 70 1.2538294 2.5594221 0.009574199
## eig_mean close_mean pr_mean clust_mean core_mean
## Endothelium 0.08769089 0.076751930 0.0044247788 0.4307828 3.274336
## Immune 0.11417292 0.009949303 0.0009661836 0.5686163 20.034783
## N.Distal_EN 0.16628860 0.024394510 0.0020533881 0.5921031 12.209446
## N.DistTubule_Dev 0.05361494 0.184172561 0.0014619883 0.2235913 3.146199
## N.LoH 0.07366429 0.124954379 0.0016835017 0.3517280 5.272727
## N.Medial_EN 0.13659061 0.033464385 0.0054945055 0.4441646 4.120879
## N.NPC 0.21121140 0.139177465 0.0222222222 0.5411640 3.222222
## N.NPC_CC 0.07574089 0.213721636 0.0020242915 0.2393817 3.548583
## N.NPC_Primed 0.12592765 0.078387774 0.0149253731 0.5123113 3.268657
## N.NPC_PTA 0.11879757 0.017717310 0.0029154519 0.4787526 6.087464
## N.NPC_PTA_CC 0.11298276 0.122961194 0.0021978022 0.3720808 6.134066
## N.PEC 0.12911455 0.058094849 0.0016207455 0.4272671 11.525122
## N.Pod_Dev 0.15750129 0.059568399 0.0050761421 0.4437689 5.040609
## N.Pod_Mat 0.10096463 0.246141397 0.0123456790 0.2589947 2.012346
## N.ProxTubule_Dev 0.17259605 0.076342010 0.0056179775 0.4469548 4.337079
## N.ProxTubule_Mat 0.09116109 0.058922519 0.0011862396 0.4303334 7.540925
## N.RV 0.12203295 0.026150774 0.0029239766 0.4722308 5.122807
## S.IC 0.06203789 0.043852421 0.0073529412 0.4246327 2.389706
## S.Med 0.16589813 0.121801532 0.0151515152 0.5115539 2.727273
## S.Mesangial 0.10222743 0.078366429 0.0023923445 0.3452017 3.523923
## S.OC_NZ 0.13981896 0.123530911 0.0125000000 0.4682616 2.612500
## U.Cortical 0.08451692 0.079404824 0.0031545741 0.4259823 4.451104res$paths # files/folders created## $run_dir
## [1] "/mnt/market/anclab-rstudio-server/home/mpir0002/Abir/sc_omicsDNA/omicsDNA_results/layer_metrics_2025-10-07_003202"
##
## $summary_csv
## [1] "/mnt/market/anclab-rstudio-server/home/mpir0002/Abir/sc_omicsDNA/omicsDNA_results/layer_metrics_2025-10-07_003202/layer_metrics_summary.csv"
##
## $nodes_xlsx
## [1] "/mnt/market/anclab-rstudio-server/home/mpir0002/Abir/sc_omicsDNA/omicsDNA_results/layer_metrics_2025-10-07_003202/layer_metrics_nodes.xlsx"
##
## $per_layer_dir
## NULL\(~\)
annotateCom() takes a community assignment table (e.g.,
from detectCom) and enriches it with metadata about each
actor. In your case, the input was comm, a
long-format table listing actors, their assigned community
(com/cid), and layers, plus
genes_info, a metadata table mapping each
gene (GeneName) to its type (GeneType). The
function outputs a community table with a new GeneType
column, so each actor is annotated with its category (e.g., lncRNA,
protein_coding, or Unknown). Key arguments include
featureID_col = "GeneName" to tell the function which
column in genes_info contains IDs, and
attribute = "GeneType" to specify which metadata attribute
to append. To improve matching, actor_normalize
standardizes IDs (removing version suffixes, trimming whitespace,
converting to lowercase), while fill_missing = "Unknown"
ensures that unmatched actors are still labeled rather than dropped. In
your run, table(comm_annot$GeneType) confirmed that after
annotation, genes were classified into three categories:
lncRNA, protein_coding, and Unknown.
comm_annot <- annotateCom(
communities = comm,
nodesMetadata = genes_info,
featureID_col = "GeneName",
attribute = "GeneType",
actor_normalize = c("strip_version","trim","tolower"),
fill_missing = "Unknown",
write_csv = TRUE,
output_prefix = "communities_with_GeneType", # alias kept for convenience
verbose = TRUE
)
# Quick sanity checks:
head(comm_annot)## actor com layer method GeneType
## 1 A2M C5 Endothelium louvain protein_coding
## 2 ACTB C3 Endothelium louvain protein_coding
## 17 ACTG1 C6 Endothelium louvain protein_coding
## 22 AFF3 C2 Endothelium louvain protein_coding
## 50 ATP1B1 C1 Endothelium louvain protein_coding
## 70 AGTR2 C2 Endothelium louvain protein_codingtable(comm_annot$GeneType, useNA = "ifany")##
## lncRNA protein_coding Unknown
## 120 6987 702\(~\)
The function sumComFeat() is designed to summarise how
particular feature types (e.g., protein-coding genes, lncRNAs, or
transcription factors) are distributed across communities in a
multilayer network. It takes as input a community assignment table such
as comm_annot, which contains actor IDs, layer membership,
community identifiers, and annotations like GeneType. By
specifying feature_col = "GeneType" and choosing one or
more feature_type values (for example
"protein_coding", "lncRNA", or
c("lncRNA","TF")), the function filters the communities to
only those actors of interest. It then produces three concise outputs: a
per-actor summary showing the communities and layers where each selected
feature appears, a per-community composition summary listing which
features belong to each community, and a count table of unique features
per community × layer.
# protein-coding
pc_summary <- sumComFeat(comm_annot, feature_col = "GeneType", feature_type = "protein_coding" , write_csv = TRUE)
# lncRNA
lnc_summary <- sumComFeat(
communities = comm_annot,
feature_col = "GeneType",
feature_type = "lncRNA"
)
attr(lnc_summary, "files") # where the three CSVs were written## NULL# multiple feature types at once
mix_summary <- sumComFeat(comm_annot, feature_col = "GeneType", feature_type = c("protein_coding","lncRNA"), write_csv = TRUE)\(~\)
get_FeatureDeg() quantifies network
centrality—specifically degree—for a user-supplied set of actors
(featureList, here all protein-coding genes extracted from
comm_annot) across selected layers
(layers = all_layers) of the multilayer network
net. Core arguments include directed = FALSE
and mode = "all" (treat edges as undirected and count all
incident edges). The function returns a table of per-layer degrees for
the requested features.
# Use any feature set you like (lncRNAs, TFs, enzymes, …)
feat_list <- unique(subset(comm_annot, GeneType == "protein_coding")$actor)
head(feat_list)## [1] "A2M" "ACTB" "ACTG1" "AFF3" "ATP1B1" "AGTR2"all_layers <- multinet::layers_ml(net)
deg_genes <- get_FeatureDeg(
net, featureList = feat_list, layers = all_layers,
directed = FALSE, mode = "all",
write_csv = TRUE,
output_file = "protein_coding_degree_byLayer.csv"
)
deg_genes[ , 1:10]## Layer A2M ACTB ACTG1 AFF3 ATP1B1 AGTR2 BEX2 CADM1 B2M
## 1 N.NPC_Primed 0 0 0 0 0 0 0 0 0
## 2 N.NPC_PTA 0 1 28 0 0 0 0 3 0
## 3 N.ProxTubule_Dev 0 0 0 0 0 0 0 0 0
## 4 N.ProxTubule_Mat 0 1 33 1 13 0 0 0 0
## 5 S.IC 0 5 2 0 1 1 2 0 0
## 6 S.Med 0 0 0 0 0 0 0 5 0
## 7 S.Mesangial 3 8 17 0 2 0 0 0 1
## 8 N.DistTubule_Dev 6 0 13 8 4 0 2 0 0
## 9 N.NPC_CC 0 1 2 0 0 0 0 0 0
## 10 S.OC_NZ 0 0 0 0 1 0 0 0 0
## 11 U.Cortical 0 3 7 0 21 1 0 0 15
## 12 N.RV 0 5 0 0 16 0 0 0 0
## 13 N.Medial_EN 0 0 0 0 3 0 0 8 0
## 14 Immune 162 91 70 16 87 0 94 158 7
## 15 N.LoH 0 1 24 12 44 14 0 1 0
## 16 Endothelium 1 5 4 6 1 3 16 3 3
## 17 N.Distal_EN 0 8 2 1 13 0 0 157 0
## 18 N.NPC_PTA_CC 1 1 1 8 1 1 0 0 1
## 19 N.PEC 0 23 30 0 99 0 0 0 0
## 20 N.Pod_Dev 1 1 1 0 0 0 0 1 0
## 21 N.Pod_Mat 0 2 1 0 0 0 0 0 2
## 22 N.NPC 0 0 0 0 0 0 0 0 0\(~\)
\(~\)
4 Comparative visualisation of multi-layers
animate_multiNet() function renders an
animated traversal through the layers of a multilayer network
(net). In its minimal form it returns an animation object
that shows the same node–edge layout while the active layer changes over
time / layer-to-layer. The output is an animation object suitable for
preview in RStudio and, if a file/format is specified, it writes a
timestamped movie under the results directory.
# Minimal run (no community colouring)
# anim <- animate_multiNet(net)
# With community colouring (using your detectCom() output)
anim <- animate_multiNet(
net,
layer_order = names(num_samples_comm),
communities = comm, # must have actor, layer, com (or cid)
layout = "kamadakawai", # typo-safe: "force", "fr" also map here
slice.par = list(start=0, interval=1, aggregate.dur=1),
seed = 1
)\(~\)
\(~\)
filmstrip_multiNet() produces a grid
(“filmstrip”) of static per‑layer network panels for a
multilayer net, optionally coloured by communities
(communities).
fs <- filmstrip_multiNet(
net,
layer_order = names(cons_list),
communities = comm,
layout = "kamadakawai",
format = "png",
width = 12, height = 8, dpi = 300,
ncol = 4,
vertex.cex = 0.8,
label.cex = 0.12,
label.col = "black",
displaylabels = TRUE,
edge.col = "#55555555",
seed = 1,
show_in_rstudio = TRUE, # <<— now it also appears in the Plots pane
verbose = TRUE
)\(~\)
plotActivityTimeline() function
summarises how actors or edges persist across layers by
drawing a longitudinal “activity” chart. In your call,
type = "edge" asks for an edge‑centric timeline and
top_n = 300 limits plotting to the 300 edges with the
longest span (or highest activity) to keep the figure legible;
edge.col sets the colour of the ribbons/segments.
layer_order aligns the x‑axis with the rest of your
pipeline. The output is a single figure (and often a companion CSV if
requested) that makes it easy to see which interactions recur across
developmental stages or conditions.
# Longest 300 edges only (keeps plot manageable)
res_edges <- plotActivityTimeline(
layer_order = names(cons_list),
net,
type = "edge",
top_n = 300,
edge.col= "#2c7fb8"
)
\(~\)
\(~\)
plot_layer_diff() function compares two
layers of the multilayer network head‑to‑head, classifying edges as
only in L1, only in L2, or
common.
num_samples
# Minimal: just draw to the Plots pane (no file), undirected, default styling
out1 <- plot_layer_diff(net, "N.ProxTubule_Dev", "N.ProxTubule_Mat", layout = "fr", vertex_label = FALSE)
# Inspect the CSV-ready data frame (weights now populated correctly)
#head(out1$edges)
# Save figure + CSV with custom aesthetics
out2 <- plot_layer_diff(
net, "N.ProxTubule_Dev", "N.ProxTubule_Mat",
layout = "kk",
col_only_L1 = "#d62728",
col_only_L2 = "#1f9d55",
col_common = "#7f8c8d",
edge_alpha = 0.55,
edge_width_only = 2.0,
edge_width_common = 1.0,
vertex_size = 4.5,
vertex_col = "grey80",
vertex_frame_col = NA,
vertex_label = TRUE,
label_cex = 0.5,
legend_pos = "bottomright",
seed = 1,
file = "diff_M3_vs_M4.png", # saved under results_dir if relative
format = "png",
width = 10, height = 7, dpi = 300,
save_edge_csv = TRUE, # CSV written alongside
csv_file = "edges_diff_M3_vs_M4.csv",
weight_attr = c("weight","prop_present","w"), # try 'prop_present' if needed
weight_aggregate = "sum"
)
\(~\)
\(~\)
The following function,
grid_layer_diffs(), higher‑level utility
automates pairwise, consecutive layer comparisons
(e.g., E1→E2→…): it loops over adjacent layer pairs in
layer_order, renders a per‑pair difference
figure for each (using the same aesthetics you would pass to
plot_layer_diff()), assembles a timestamped grid
image of all pairs (grid_format), and writes a
single combined CSV stacking all per‑pair edge
tables.
# A) Minimal, PNG everywhere (all outputs timestamped)
#gmin <- grid_layer_diffs(
# net,
# layer_order = names(cons_list),
# ncol = 4,
# vertex_label = FALSE,
# pair_format = "png",
# grid_format = "png"
#)
# B) Vector per‑pair (PDF) + JPEG grid, custom aesthetics
gfull <- grid_layer_diffs(
net,
layer_order = names(cons_list),
ncol = 4,
layout = "kk",
col_only_L1 = "#d62728",
col_only_L2 = "#1f9d55",
col_common = "#7f8c8d",
edge_alpha = 0.55,
edge_width_only = 2.0,
edge_width_common = 1.2,
vertex_size = 4.0,
vertex_col = "grey85",
vertex_frame_col = NA,
vertex_label = T,
legend_pos = "bottomright",
seed = 1,
grid_format = "jpg",
pair_format = "pdf",
pair_file_prefix = "diff_pair",
panel_width = 6, panel_height = 5,
combined_csv_file = "diff_edges_all_pairs.csv" # will get timestamp appended
)
\(~\)
animate_multiNet_mp4(): This renderer
creates a video (MP4) or animated GIF
of the multilayer network. You pass the net, optionally
communities (with actor, layer,
com/cid) to colour nodes consistently, and
choose the layout (e.g., "kamadakawai").
Temporal resolution is controlled by fps and
frames_per_layer; canvas size by width and
height; and show.labels toggles node labels.
The format argument selects "gif" or
"mp4".
# GIF
res_gif <- animate_multiNet_mp4(
net,
communities = comm, # optional; actor, layer, com/cid
layout = "kamadakawai",
format = "gif",
fps = 12,
frames_per_layer = 12,
width = 1000, height = 800,
show.labels = FALSE
)
res_gif$file # path to the GIF
# MP4
res_mp4 <- animate_multiNet_mp4(
net, communities = comm, format = "mp4", fps = 15, frames_per_layer = 10
)\(~\)
\(~\)
5 Per cell-type enrichment analysis
go_enrichment_report() — per‑layer (multinet)
mode, based on clusterProfiler This call performs
over‑representation testing separately for each layer of the
multilayer network using clusterProfiler’s
enrichGO. The function first harmonizes gene identifiers to
ENTREZ (via OrgDb = org.Hs.eg.db,
id_from = "SYMBOL"), then runs GO enrichment for the
Biological Process ontology (ont = "BP"). With
universe_mode = "union", the statistical background is the
union of all genes mapped across layers, improving comparability of
p‑values between layers. Plotting options generate per‑layer
dotplots and cnet plots (barplots suppressed by
include_barplot_multinet = FALSE), with label spacing tuned
by cnet_left_margin_pt and the number of terms shown
controlled by cnet_showCategory. Results are previewed in
RStudio (show_in_rstudio = TRUE) and saved as one
combined, timestamped PDF (save_pdf = TRUE,
one_pdf = TRUE, file_prefix = "GO_BP_byLayer",
out_dir = getOption("mlnet.results_dir", …)). The output
object contains, for each layer, the enrichResult table,
the plots used for export, and file paths to the saved artifacts.
go_enrichment_report() — single set mode, based
on clusterProfiler Here the same function is used in
“single‑set” mode by supplying a vector of genes
(genes = g) rather than a network, again mapping
SYMBOL to ENTREZ and running enrichGO for
GO:BP. The call yields
an enrichResult object plus publication‑ready plots
(dotplot and cnet; a barplot can be included by default in single‑set
mode), displayed in RStudio and saved using the provided prefix
(file_prefix = "g5l7_1_GO_BP"). The output includes the
enrichment table (adjusted p‑values per term), the plots, and the paths
to timestamped files written under the results directory.
# Per‑layer enrichment; dot + cnet per layer; NO barplot in multinet
out <- go_enrichment_report(
net = net,
OrgDb = org.Hs.eg.db,
id_from = "SYMBOL",
ont = "BP",
universe_mode = "union", # background = union of layer genes (mapped)
include_barplot_multinet = FALSE, # now recognized by the function
cnet_showCategory = 8,
cnet_left_margin_pt = 100,
show_in_rstudio = TRUE,
save_pdf = TRUE,
one_pdf = TRUE, # one combined PDF (no duplicates)
file_prefix = "GO_BP_byLayer",
out_dir = getOption("mlnet.results_dir","omicsDNA_results")
)
#g <- c("RHCG","SPRR1A","SPRR2A","SPRR3","KRT13","KRT4","CRNN")
#res <- go_enrichment_report(genes = g, OrgDb = org.Hs.eg.db, ont = "BP", file_prefix = "g5l7_1_GO_BP", show_in_rstudio = TRUE)You can see ClusterProfiler enrichment results for loop of Henle cells as follow:
\(~\)
\(~\)
gp_enrich_multinet() — per‑layer ORA based on
g:Profiler (gprofiler2) This function
performs per‑layer over‑representation analysis via
g:Profiler (through gprofiler2), targeting the
human organism (organism = "hsapiens") and multiple
knowledge bases (sources = c("GO:BP","REAC")). Layers are
analysed in the supplied order
(layer_order = names(cons_list)), and the top enriched
terms per layer are visualised (show_terms = 8). Beyond
per‑layer tables and dot‑style summaries, it can also render a
term‑only overlap network (TON) across layers
(ton_show = TRUE) with display limits
(ton_top_terms, ton_min_overlap) and layout
(ton_layout = "fr"). Figures are exported with the chosen
graphics settings (format, width,
height, dpi) and previewed in RStudio
(show_in_rstudio = TRUE). The return typically includes
per‑layer enrichment tables, generated figures, the TON visualization,
and a manifest of file paths under a timestamped output folder.
out <- gp_enrich_multinet(
net,
organism = "hsapiens",
sources = c("GO:BP","REAC"),
layer_order = names(cons_list),
show_terms = 8,
ton_show = TRUE,
ton_top_terms = 12,
ton_min_overlap = 1,
ton_layout = "fr",
format = "png",
width = 10, height = 7, dpi = 300,
show_in_rstudio = TRUE
)
\(~\)
You can see gprofiler enrichment results for developing proximal tubule cells as follow:
\(~\)
\(~\)
gsea_msigdbr_layer_pairs() — pairwise GSEA based
on msigdbr + fgsea This routine carries out
rank‑based gene set enrichment analysis between
consecutive layer pairs, drawing gene‑set definitions
from MSigDB via msigdbr
(msigdb_species = "Homo sapiens",
msigdb_collections = c("H","C2","C5"),
msigdb_subcategories = c("CP:KEGG","GO:BP")) and performing
GSEA with fgsea (by default nperm = NULL triggers
fgseaMultilevel). Per‑gene rankings are computed from
vertex attributes (e.g., cor_attr = c("cor","r"),
padj_attr = c("padj","FDR","q")) using a specified score,
here pair_rank_by = "delta_signed_log10_padj". The function
writes a timestamped run directory
(results_dir, run_name) containing per‑pair
dotplots, enrichment curves,
cnet plots (cnet_show = TRUE, with layout
and label size controls), and an optional term‑only overlap
network (ton_show = TRUE).
out <- gsea_msigdbr_layer_pairs(
net,
layer_order = names(cons_list), # keep your preferred order
# where to read per-gene stats from vertex attributes:
cor_attr = c("cor","r"),
padj_attr = c("padj","FDR","q"),
# ranking choice:
pair_rank_by = "delta_signed_log10_padj",
# MSigDB settings (start light; you can add more collections later)
msigdb_species = "Homo sapiens",
msigdb_collections = c("H","C2","C5"),
msigdb_subcategories = c("CP:KEGG","GO:BP"),
minSize = 15,
maxSize = 500,
nperm = NULL, # <- set to NULL for fgseaMultilevel or set e.g. 10000 to use classic permutations
# output container (timestamp appended automatically)
results_dir = getOption("mlnet.results_dir","omicsDNA_results"),
run_name = "gsea_pairs_demo",
# plots: sizes & margins
dotplot_width = 8, dotplot_height = 6,
enrich_width = 8, enrich_height = 6,
cnet_width = 10, cnet_height = 8,
dotplot_margin_cm = c(0.8,0.8,0.8,0.8),
enrich_margin_cm = c(0.8,0.8,0.8,0.8),
cnet_margin_cm = c(0.8,0.8,0.8,0.8),
# cNET: terms + gene symbols (forced by the function)
cnet_show = TRUE,
cnet_top_terms = 12,
cnet_layout = "fr",
cnet_term_label_size = 3.2,
cnet_gene_label_size = 2.4,
cnet_gene_label_color = "grey25",
# term-only overlap network
ton_show = TRUE,
ton_top_terms = 12,
ton_min_overlap = 1,
ton_layout = "fr",
show_in_rstudio = TRUE,
seed = 1,
verbose = TRUE
)
# Where outputs went:
out$run_dir
read.csv(file.path(out$run_dir, "SUMMARY.csv"))
\(~\)
You can see GSEA enrichment results for comparison between maturing proximal tubule versus developing proximal tubule cells as follow:
\(~\)
\(~\)