Wrapper Function

The differential_analysis function runs the full differential analysis workflow in one call, covering three sequential steps: statistical testing, visualization, and pathway annotation. For users who prefer fine-grained control over parameters, please refer to the detailed parameter settings below.

First, count matrices across samples are compared condition-by-condition using a two-round limma-voom framework, identifying genomic regions with significant differences in accessibility. When column clusters are provided, counts within each cluster are aggregated across pairs before testing, and results are reported per cluster. Each pairwise comparison produces a set of significant regions with associated statistics.

Next, the significant regions from each comparison are visualized as log2 count heatmaps, with columns grouped by sample. This allows direct visual inspection of accessibility patterns across conditions.

Finally, the significant regions are submitted to GREAT for pathway enrichment analysis against a user-specified MSigDB gene set collection, revealing the biological processes associated with each set of differential regions. Results are summarised in per-comparison TSV tables and an optional bubble plot.

Parameters

Parameter Type Default Description Example
cm_path character vector — Paths to per-sample count matrices in .feather format, one per sample. Must have the same length as conditions. Each file must contain a pos column and one column per CRF pair. c(“C1.feather”, “C2.feather”, “T1.feather”, “T2.feather”)
conditions character vector — Condition labels for each sample in cm_path, in the same order. Each condition must have at least 2 replicates. All pairwise comparisons between conditions are tested. c(“C”, “C”, “T”, “T”)
sample_names character vector or NULL NULL Optional sample names used in output file names and heatmap column labels. Must be unique and the same length as conditions. If NULL, names are auto-generated as <condition><replicate_index> (e.g. C1, C2, T1). c(“Control_1”, “Control_2”, “Treat_1”, “Treat_2”)
out_dir character “./” Root output directory for all results. Differential region tables and heatmaps are written here; pathway enrichment outputs are written to <out_dir>/pathway_enrichment/. Created recursively if it does not exist. “./differential”
col_cluster_file_path character or NULL NULL Path to a TSV file with columns pair and cluster, mapping CRF pairs to column clusters. When provided, counts within each cluster are summed across pairs before testing, and heatmaps expand back to pair-level columns grouped by sample. When NULL, each pair is tested independently. “bicluster/col_table.tsv”
ref_genome character “hg38” Reference genome. Controls the TSS source for GREAT and the species for MSigDB gene set retrieval. Supported: “hg38”, “mm10”. “mm10”
min_support integer 2 Minimum number of non-zero replicates required within at least one condition for a genomic region to pass the initial non-zero support filter before differential testing. Only regions passing this filter are carried forward to the limma-voom modeling step and considered for significance testing. For example, with min_support = 2, a region is retained if it has non-zero counts in at least 2 replicates in Control or at least 2 replicates in Treatment. 1
lfc_threshold numeric 0.5 Minimum absolute log2 fold-change required for a region to be called significant after model fitting. A region is considered significant only if it satisfies both abs(logFC) > lfc_threshold and adj.P.Val < fdr_threshold. This parameter controls the effect size requirement in the final significance decision. 1
fdr_threshold numeric 0.05 Adjusted p-value cutoff used to define statistically significant regions after multiple-testing correction. Together with lfc_threshold, this parameter determines the final set of significant differential regions. 0.01
apply_annotation logical TRUE Whether to perform GREAT-based pathway enrichment annotation on significant differential regions. Set to FALSE to skip the pathway annotation step entirely. FALSE
plot logical TRUE Whether to generate visualizations. Controls both heatmap generation and the pathway enrichment bubble plot. Set to FALSE to suppress all plot outputs. FALSE

Example Usage

library(multiEpiCore)
# Test Data
cm_path    <- c("count_matrix/C1_Count_Matrix_800.feather",
                "count_matrix/C2_Count_Matrix_800.feather",
                "count_matrix/T1_Count_Matrix_800.feather",
                "count_matrix/T2_Count_Matrix_800.feather")
conditions <- c("C", "C", "T", "T")
col_cluster_path <- "bicluster/col_table.tsv"
out_dir <- "./differential"

differential_analysis(cm_path = cm_path, conditions = conditions, col_cluster_file_path = col_cluster_path, out_dir = out_dir)

1. Differential Region Detection

The differential_regions() function identifies differential genomic regions between experimental conditions using limma-voom. It supports two distinct analysis modes depending on whether a column cluster file is provided: 1. Cluster-aware mode (recommended for grouped CRF analysis)
- The function reads a two-column table (pair, cluster) that defines which CRF pairs belong to each column cluster. - For each cluster: - All CRF pairs assigned to that cluster are aggregated by row-wise summation within each sample. - The aggregated matrix (regions × samples) is used for differential testing. - Differential analysis is therefore performed at the cluster level, integrating signal across biologically related CRF pairs.

  1. Per-pair mode (each CRF pair analyzed independently)
  • Each CRF pair is automatically treated as its own cluster.
  • No aggregation across pairs is performed.
  • Differential analysis is conducted independently for each CRF pair.

For each analysis group (cluster or CRF pair), differential analysis is performed using a standardized limma-voom workflow consisting of sequential filtering, normalization, and statistical modeling:

  1. Non-zero support filtering
    Regions are first filtered based on signal support across replicates.
    A region is retained only if it has non-zero counts in at least min_support replicates within at least one condition.
    This step removes regions with insufficient signal for reliable downstream modeling.

  2. Initial normalization and voom transformation
    The count matrix is normalized using TMM (calcNormFactors) and transformed to log2-CPM values via voom, which also estimates the mean–variance relationship and assigns observation-level precision weights.

  3. Mean expression filtering
    Based on the initial voom output, regions with low overall expression are removed.
    Specifically, regions with mean log2-CPM below the mean_quantile cutoff are discarded to improve statistical power and reduce noise.

  4. Re-normalization and voom re-fitting
    The filtered matrix is re-normalized and passed through voom again to obtain updated weights for the final model fitting.

  5. Linear modeling and empirical Bayes moderation (limma)
    A design matrix is constructed from the input conditions, and a linear model is fitted for each region.
    Empirical Bayes moderation (eBayes) is applied to stabilize variance estimates across regions.

  6. Differential region identification
    Regions are defined as differential if they satisfy both:

    • absolute log2 fold-change > lfc_threshold
    • adjusted p-value (FDR) < fdr_threshold

Parameters

Parameter Type Default Description Example
cm_path character vector — Paths to per-sample count matrices in .feather format. Each file should represent one sample, with rows as genomic regions (stored in a required pos column) and columns as CRF pairs. All files are harmonized to the common set of regions and CRF pairs before analysis. c(“C1.feather”, “C2.feather”, “T1.feather”, “T2.feather”)
conditions character vector — Condition labels corresponding to cm_path, in the same order. Each unique condition must have at least 2 replicates, because limma-voom is performed on replicated group comparisons. c(“Control”, “Control”, “Treatment”, “Treatment”)
sample_names character vector NULL Optional unique sample names corresponding to cm_path. If NULL, names are automatically generated from conditions in the form of condition + replicate index (for example, Control1, Control2). c(“C1”, “C2”, “T1”, “T2”)
out_dir character “./” Output directory for all result files, including full differential tables, significant-region tables, summary tables, summary plots, and significant-region count matrices. “diff_results”
col_cluster_file_path character NULL Optional path to a two-column TSV file with columns pair and cluster. If provided, the function runs in cluster-aware mode: CRF pairs assigned to the same cluster are summed within each sample and tested as one aggregated group. If NULL, the function runs in per-pair mode, where each CRF pair is tested independently. “col_cluster.tsv”
min_support integer 2 Minimum number of non-zero replicates required within at least one condition for a genomic region to pass the initial non-zero support filter. For example, with min_support = 2, a region is retained if it has non-zero counts in at least 2 replicates in Control or at least 2 replicates in Treatment. 1
mean_quantile numeric 0 Quantile cutoff applied to the mean log2-CPM values after the first voom transformation. Regions with mean expression below this cutoff are removed before the final limma fit. Larger values make the mean-expression filter more stringent. 0.1
lfc_threshold numeric 0.5 Minimum absolute log2 fold-change required for a region to be called significant after model fitting. A region must satisfy both abs(logFC) > lfc_threshold and adj.P.Val < fdr_threshold. 1
fdr_threshold numeric 0.05 Adjusted p-value cutoff used to define significant regions after multiple-testing correction. 0.01

Output Files

The function automatically derives all pairwise comparisons from the unique values in conditions (sorted alphabetically, in the form <test>_vs_<ref>), and creates one subdirectory per comparison under out_dir. All output files are written into the corresponding subdirectory:

  1. All regions - <comparison>/<cluster|pair>_all.tsv
    • Complete limma-voom differential analysis results for all tested genomic regions within the specified cluster and comparison.
    • Each row represents one genomic region and the table includes the following columns:
      • pos: Genomic region identifier corresponding to the original pos field in the count matrix.
      • logFC: Log2 fold change between the test and reference conditions (positive = higher in test).
      • AveExpr: Average log2-CPM expression level across all samples in the comparison.
      • t: Moderated t-statistic from limma’s empirical Bayes model.
      • P.Value: Raw p-value from the moderated t-test.
      • adj.P.Val: Benjamini-Hochberg FDR-adjusted p-value.
      • B: Log-odds that the region is truly differentially expressed.
pos logFC AveExpr t P.Value adj.P.Val B
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
chr14_77181601_77182400 4.6299 1.0508 9.6915 3.28e-22 1.80e-16 38.9807
chr1_205504001_205504800 3.9690 0.7699 7.9780 1.49e-15 4.07e-10 24.3462
chr16_67982401_67983200 3.9468 0.5336 7.8500 4.16e-15 6.90e-10 23.3585
chr3_12883201_12884000 3.9314 0.3721 7.8259 5.04e-15 6.90e-10 23.1788
chr16_2856001_2856800 3.8655 0.1653 7.5056 6.12e-14 6.03e-09 20.7914
…

This table shows differential analysis results for cluster 4 in the test dataset.

  1. Significant regions - <comparison>/<cluster|pair>_sig.tsv
    • Filtered results containing only significant regions (FDR < threshold & |logFC| > threshold)
    • Only generated when at least one significant region is detected under the specified thresholds.
  2. Significant region counts - <comparison>/<cluster|pair>_sig_counts.tsv
    • Raw fragment counts for significant regions, with one row per region and one column per sample (no clustering) or one column per sample:pair combination (with clustering), preserving pair-level resolution for downstream visualization.
    • Only generated alongside _sig.tsv when at least one significant region is detected.
pos C1 C2 T1 T2
C1:H3K4me3-H3K9me3 C1:H3K27me3-H3K4me3 C1:H3K4me1-H3K4me3 C1:H3K4me3-H3K4me3 C2:H3K4me3-H3K9me3 C2:H3K27me3-H3K4me3 C2:H3K4me1-H3K4me3 C2:H3K4me3-H3K4me3 T1:H3K4me3-H3K9me3 T1:H3K27me3-H3K4me3 T1:H3K4me1-H3K4me3 T1:H3K4me3-H3K4me3 T2:H3K4me3-H3K9me3 T2:H3K27me3-H3K4me3 T2:H3K4me1-H3K4me3 T2:H3K4me3-H3K4me3
chr14_77181601_77182400 0.288288288288288 0 0 0.288288288288288 0 0 0 0 1.20578778135048 10.6053869384119 5.71643061399974 10.0646506873736 2 4.00456585667853 3.11594961401672 4.08021515668152
chr1_205504001_205504800 0 0 0 0 0 0.523255813953488 0 0 1 4.36625846297515 2 3.8 1 2 3.62058405337417 8.80857631815813
chr16_67982401_67983200 0 0 0 0.320689655172414 0 0 0 0 6 0 4 3.86699507389163 0 1 3 3.87676767676768
chr3_12883201_12884000 0.062962962962963 0.062962962962963 0 0 0 0 0 0 1.96317280453258 2.02685520361991 1 3.552084637536 2 0.412213740458015 3 5.23860220144396
chr16_2856001_2856800 0 0 0 0 0 0 0 0 0 2 3.38220600617457 2.83513513513514 1 2.66066838046272 3 1
…

This table shows differential analysis results for cluster 4 in the test dataset.

  1. Summary table per comparison - <comparison>/summary.tsv Includes:
    • cluster ID
    • number of regions before filtering
    • after non-zero filter
    • after mean filter
    • total significant
    • up-regulated
    • down-regulated
cluster n_rows_before n_rows_after_nonzero n_rows_after_mean n_sig n_up n_down
cluster1 3860350 453271 453271 0 0 0
cluster2 3860350 970126 970126 14092 8085 6007
cluster3 3860350 365883 365883 0 0 0
cluster4 3860350 547339 547339 8250 5701 2549
  1. Summary plot per comparison - <comparison>/summary.pdf
  • Mirror bar chart visualizing the number of differentially accessible regions per group, with up-regulated regions shown in red (above x-axis) and down-regulated regions in blue (below x-axis).

Example Usage

library(multiEpiCore)

# Test data
cm_path <- c("count_matrix/C1_Count_Matrix_800.feather", "count_matrix/C2_Count_Matrix_800.feather", "count_matrix/T1_Count_Matrix_800.feather", "count_matrix/T2_Count_Matrix_800.feather")
conditions <- c("C", "C", "T", "T")
col_cluster_file_path <- "bicluster/col_table.tsv"
out_dir <- "differential"
differential_regions(cm_path = cm_path, conditions = conditions, col_cluster_file_path = col_cluster_file_path, out_dir = out_dir)

2. Differential Accessibility Heatmaps

The differential_heatmap() function generates per-comparison, per-group heatmaps for all significant differential regions produced by differential_regions(). It iterates over the full result object returned by differential_regions() and dispatches each comparison-group combination to differential_heatmap_single(), which renders a region × sample log2 count heatmap as a PDF file.

When the underlying analysis used column clustering (col_cluster_file_path), columns in the heatmap represent individual sample:pair combinations within the cluster, grouped and titled by sample. Without clustering, each column represents one sample directly.

Parameters

Parameter Type Default Description Example
dr_result nested list — The return value of differential_regions(). A two-level nested list structured as dr_result[[grp_name]][[cmp_tag]], where the first level keys are cluster or pair names (e.g. “cluster4”, “H3K27ac_peak1”) and the second level keys are comparison tags (e.g. “T_vs_C”). Each leaf is the file path to the corresponding _sig_counts.tsv. dr_result <- differential_regions(…)
out_dir character “./” Output directory for heatmap PDFs. Created recursively if it does not exist. One PDF is written per comparison-group combination, named <cmp_tag>_<grp_name>.pdf. “heatmaps”
show_colnames logical FALSE Whether to display individual column names on the heatmap. When FALSE, only the sample-level column group titles are shown. Set to TRUE to additionally show individual pair identifiers within each sample group. TRUE
col_width_mm numeric 40 Width per individual pair column in mm. Total heatmap width is computed as n_pairs × col_width_mm + (n_samples − 1) × 8 mm, where the second term accounts for the fixed 8 mm gap between sample groups. In cluster mode (when col_cluster_file_path was supplied to differential_regions()), each sample group contains multiple pair columns; col_width_mm still refers to a single pair column, so wider clusters automatically occupy proportionally more space. In non-cluster mode each sample has exactly one column, so col_width_mm equals the per-sample width directly. 30
row_height_mm numeric 0.5 Height per region row in mm. Total heatmap body height is n_sig_regions × row_height_mm mm, so the PDF scales automatically with the number of significant regions. Note that values below ~0.35 mm (one pixel at 72 dpi) will be visually indistinguishable; for large region counts (>1000) a value of 0.13–0.5 mm produces dense stripe patterns comparable to a fixed-height layout. 1.0
random_seed integer 42 Random seed passed to set.seed() before each raster rendering call to ensure pixel-level reproducibility across runs. 123

Output Files

Output is organized by comparison tag (cmp_tag) read from dr_result:

  1. Heatmap - <comparison>/<cluster|pair>.pdf
  • Rows: Significant genomic regions
  • Columns: All CRF pairs across all samples
  • Column blocks: Grouped by sample

This table shows differential analysis heatmap for cluster 4 in the test dataset.

Example Usage

library(multiEpiCore)

# Test data
dr_result <- list(
  cluster2 = list(
    T_vs_C = "differential/T_vs_C_cluster2_sig_counts.tsv"
  ),
  cluster4 = list(
    T_vs_C = "differential/T_vs_C_cluster4_sig_counts.tsv"
  )
)

differential_heatmap(
  dr_result = dr_result,
  out_dir   = "differential"
)

3. Differential Pathway Annotation

The differential_pathway_annotation() function performs pathway enrichment analysis on the significant differential regions from each comparison-group combination using GREAT (Genomic Regions Enrichment of Annotations Tool). GREAT assigns genomic regions to nearby genes via basal-plus-extension TSS domains, then tests for over-representation of gene sets using a binomial test. Gene sets are sourced from MSigDB via the msigdbr package, and can be configured to use any available collection such as Hallmark, Reactome, or GO terms.

Each comparison-group combination (e.g. T vs C within cluster 4) is tested independently. Results are written as per-comparison TSV tables and optionally visualised as a bubble plot that summarises enrichment patterns across all comparisons at a glance.

Parameters

Parameter Type Default Description Example
dr_result nested list — The return value of differential_regions(). A two-level nested list structured as dr_result[[grp_name]][[cmp_tag]], where each leaf is the file path to the corresponding _sig_counts.tsv. Each comparison-group combination becomes one GREAT query set. dr_result <- differential_regions(…)
out_dir character “./” Output directory for enrichment TSV tables and the bubble plot PDF. Created recursively if it does not exist. One TSV is written per query set, named pathway_annotation_<cmp_tag>_<grp_name>.tsv. “pathway”
ref_genome character “hg38” Reference genome. Controls both the TSS source used by GREAT (“hg38” → TxDb.Hsapiens.UCSC.hg38.knownGene, “mm10” → TxDb.Mmusculus.UCSC.mm10.knownGene) and the species passed to msigdbr(). Must be one of “hg38” or “mm10”. “mm10”
msigdb_collection character “H” MSigDB collection abbreviation. Common values: “H” (Hallmark, 50 curated biological states), “C2” (curated gene sets including KEGG and Reactome), “C5” (GO gene sets). See msigdbr::msigdbr_collections() for the full list. “C2”
msigdb_subcollection character or NULL NULL MSigDB sub-collection abbreviation. Required when msigdb_collection has subcollections. For example, “CP:REACTOME” or “CP:KEGG_LEGACY” under “C2”, and “GO:BP”, “GO:MF”, or “GO:CC” under “C5”. Leave NULL for collections without subcollections such as “H”. “CP:REACTOME”
plot logical TRUE Whether to generate a bubble plot PDF (pathway_annotation.pdf) summarising enrichment across all query sets. Each row is a pathway ordered by best adjusted p-value across query sets; bubble size encodes log2(1 + fold_enrichment) and color encodes −log10(padj). FALSE

Output Files

Output is organized by comparison tag (cmp_tag) read from dr_result, with a pathway_annotation/ subdirectory created under each:

  1. Enrichment tables — <cmp_tag>/pathway_annotation/pathway_annotation_<grp_name>.tsv
  • pathway : Pathway or gene set name, with collection-specific prefixes removed (e.g. HALLMARK_ stripped for collection "H")
  • hits_region: Number of input regions associated with at least one gene in the gene set
  • fold: Fold enrichment of observed over expected region hits
  • p: Binomial test p-value
  • padj: BH-adjusted p-value across all tested gene sets
  • hits_gene: Number of genes in the gene set that were hit by at least one input region
pathway hits_region fold p padj hits_gene
<chr> <int> <dbl> <dbl> <dbl> <int>
ESTROGEN_RESPONSE_LATE 264 1.74372272845793 0 0 103
ESTROGEN_RESPONSE_EARLY 309 1.66813846271819 0 0 112
HEME_METABOLISM 210 1.73048512649459 1.20237153566904e-13 2.00395255944841e-12 112
P53_PATHWAY 214 1.66542359524271 2.27673435659881e-12 2.84591794574851e-11 99
ADIPOGENESIS 180 1.47667160397026 4.18145192693231e-07 4.18145192693231e-06 88
APICAL_JUNCTION 202 1.41764920800881 1.24266646217563e-06 1.03555538514636e-05 94
UV_RESPONSE_UP 160 1.47287922146693 2.02699119755678e-06 1.4056939592888e-05 72
TNFA_SIGNALING_VIA_NFKB 255 1.34870029318323 2.24911033486208e-06 1.4056939592888e-05 109
XENOBIOTIC_METABOLISM 153 1.40379208431073 3.55465706951552e-05 0.000197480948306418 77
KRAS_SIGNALING_DN 207 1.28886497608787 0.000218809635529782 0.00109404817764891 82
MYOGENESIS 189 1.29612613123138 0.00030605164274955 0.00139114383067977 94
APICAL_SURFACE 68 1.51450373497597 0.000763245234423948 0.00318018847676645 22
IL6_JAK_STAT3_SIGNALING 71 1.44814605261822 0.00182060936451056 0.00700234370965601 38
COMPLEMENT 182 1.24466014163988 0.00218293749453657 0.00779620533763062 85
APOPTOSIS 153 1.26034481113274 0.00298056189986595 0.00993520633288651 78
UNFOLDED_PROTEIN_RESPONSE 86 1.36360512135354 0.00337278018971299 0.0101202710319365 46
COAGULATION 106 1.31851728359683 0.00344089215085841 0.0101202710319365 55
HYPOXIA 208 1.20775475781321 0.00408508124078466 0.0113474478910685 99
PEROXISOME 83 1.30331932988064 0.0112023518361346 0.0294798732529858 43
DNA_REPAIR 80 1.30205773702884 0.0128044845936148 0.0320112114840371 55
ALLOGRAFT_REJECTION 148 1.18870915581885 0.0210657421315557 0.049948131180146 68
IL2_STAT5_SIGNALING 203 1.15594723983505 0.0219771777192642 0.049948131180146 103
NOTCH_SIGNALING 41 1.352751222284 0.0365638577136234 0.079486647203529 22
MYC_TARGETS_V2 34 1.38116190236135 0.0417371750854787 0.0869524480947473 19
HEDGEHOG_SIGNALING 70 1.23591001049887 0.0467664045631269 0.0923888270525026 23
INTERFERON_GAMMA_RESPONSE 153 1.14886463973552 0.0480421900673014 0.0923888270525026 82
TGF_BETA_SIGNALING 69 1.21829631773623 0.0602666642784661 0.111604933849011 27
ANDROGEN_RESPONSE 108 1.15527635776966 0.0748601474518107 0.133678834735376 45
GLYCOLYSIS 171 1.11257302908581 0.0873055600289641 0.150526827636145 91
REACTIVE_OXYGEN_SPECIES_PATHWAY 28 1.27497045636489 0.12046585317067 0.200776421951118 20
MITOTIC_SPINDLE 159 1.08293783515106 0.165188855776895 0.26643363834983 87
PI3K_AKT_MTOR_SIGNALING 84 1.09632874745783 0.212732093921961 0.332393896753064 45
ANGIOGENESIS 40 1.14144142997914 0.221601898425801 0.335760452160304 13
INTERFERON_ALPHA_RESPONSE 56 1.11265243274659 0.229055593039617 0.336846460352377 29
MTORC1_SIGNALING 136 1.06364076045101 0.245581984444253 0.350831406348933 73
FATTY_ACID_METABOLISM 103 1.05800821660232 0.295711315823185 0.410710160865535 55
OXIDATIVE_PHOSPHORYLATION 105 1.04256584296824 0.346946038906058 0.4688459985217 62
BILE_ACID_METABOLISM 87 1.04364869883548 0.358943149554168 0.472293617834432 39
SPERMATOGENESIS 111 1.02500662584779 0.409614166640853 0.525146367488273 56
E2F_TARGETS 107 1.02208107671591 0.423229269291039 0.529036586613799 65
PROTEIN_SECRETION 77 1.0191626411716 0.448989029511285 0.547547596964982 39
CHOLESTEROL_HOMEOSTASIS 53 1.01810588612598 0.466391526722058 0.55522800800245 34
WNT_BETA_CATENIN_SIGNALING 46 0.979754174443012 0.574991270377173 0.668594500438573 21
UV_RESPONSE_DN 219 0.965380875723917 0.71059762307301 0.807497298946602 80
INFLAMMATORY_RESPONSE 159 0.936597083120077 0.808037021503829 0.897249773334346 82
MYC_TARGETS_V1 85 0.908889894979314 0.825469791467598 0.897249773334346 56
G2M_CHECKPOINT 136 0.885947916228726 0.931017173094037 0.990443801163869 78
PANCREAS_BETA_CELLS 33 0.738326975220715 0.971080305832514 0.994386024218582 18
EPITHELIAL_MESENCHYMAL_TRANSITION 206 0.859020240102634 0.989082044671641 0.994386024218582 91
KRAS_SIGNALING_UP 185 0.837898877967919 0.994386024218582 0.994386024218582 83

This table shows pathway annotation result for cluster 4 in the test dataset.

  1. Bubble plot — <cmp_tag>/pathway_annotation/pathway_annotation.pdf (when plot = TRUE)
  • Rows are pathways ordered by best adjusted p-value across all clusters/pairs.
  • Columns are clusters/pairs.
  • Bubble size encodes log2 fold enrichment and color encodes -log10(padj).

Example Usage

library(multiEpiCore)

# Test data
dr_result <- list(
  cluster2 = list(
    T_vs_C = "differential/T_vs_C_cluster2_sig_counts.tsv"
  ),
  cluster4 = list(
    T_vs_C = "differential/T_vs_C_cluster4_sig_counts.tsv"
  )
)

differential_pathway_annotation(
  dr_result = dr_result,
  out_dir   = "differential/pathway_annotation"
)