MurineMatchTx Manual

Author: Dian Wang, Guiwen Guan

Introduction

The Cross-Species Transcriptomic Similarity Analysis R package (MurineMatchTx) provides a comprehensive and user-friendly framework for evaluating transcriptional similarities between mouse cancer models and human cancer samples. Designed to streamline cross-species transcriptomic analysis, MurineMatchTx enables rapid, multidimensional exploration with minimal preprocessing—simply input standardized files following the provided guidelines.

The package integrates a wide array of functionalities, including differential expression gene integration and visualization, graphical display of PCA, transcriptional similarity scoring (TROM scores), immune infiltration profiling, and functional enrichment analysis (GO terms). Each function is highly customizable through flexible parameter settings, allowing researchers to tailor the analysis to specific experimental needs.

With MurineMatchTx, users can rigorously assess how closely animal models recapitulate the molecular features of human cancers, thereby facilitating translational research at the transcriptomic level.

Installation

You can install the package directly from GitHub:

# Install devtools if necessary
if (!require("devtools")) install.packages("devtools")

# Install the package
devtools::install_github("wangdian-PKU/MurineMatchTx")

General input files to be prepared in advance

We strongly suggest that all user - provided input files should be in the tsv format!!!

1. Expression matrix: Differentially expressed gene matrix between the cancer group and the normal group in mouse models.
2. Expression matrix: Differentially expressed gene matrix between human cancer group and the normal group (TCGA as examples).
3. Clinical Data: human clinical metadata with patient IDs and different clinical conditions data.
4. CNV Data: human CNV alterations of specified genes (e.g., TP53, PTEN).

Main Functions

1. Merge_DEG_datasets()— Merge differentially expressed data from human samples or mouse models

1.1 Purpose

The merge_DEG_datasets() function is designed to combine differential gene expression (DEG) results from TCGA (human) and multiple mouse cancer models into a unified, tidy data frame. This merged data can then be used to draw volcano plots of DEG downstream.

1.2 Input Data Format

This function requires DEG result files from TCGA and mouse models in .tsv or .csv format. It can be derived from the DEG matrix generated by DESeq2 or edgeR. These files must:

1. Be tab- or comma-separated.
2. Include at least the following columns:
   • logFC: Log2 fold change of gene expression.
   • FDR: False discovery rate (adjusted p-value).
3. Row names should be gene symbols or gene identifiers.

Example .tsv or .csv format after being read:

1.3 Function Parameters

1.4 Output

A merged data frame with the following columns:

1.5 DEG Classification Logic

The function uses the following rule to classify gene expression change:

if (logFC > threshold & FDR < threshold) → "Up"
if (logFC < -threshold & FDR < threshold) → "Down"
otherwise → "Stable"

1.6 Run Example

# Define file paths
tcga_file <- "./data/TCGA_DEG.tsv"

mouse_files <- list(
  "Our_Model"   = "./data/Our_Model_DEG.csv",
  "GSE172629"   = "./data/GSE172629_DEG.tsv",
  "GSE208279"   = "./data/GSE208279_DEG.tsv"
)

# Merge datasets with default thresholds
merged_DEG <- merge_DEG_datasets(
  human_file_path = tcga_file,
  mouse_files = mouse_files
)

# Optional: custom thresholds
merged_DEG_custom <- merge_DEG_datasets(
  human_file_path = tcga_file,
  mouse_files = mouse_files,
  logFC_value = 1.5,
  FDR_value = 0.01
)

# View result
head(merged_DEG)

1.7 Downstream Use Cases

1. As input for plot_DEG_volcano() functions.
2. For cross-species DEG pattern comparison.

1.8 Notes

1. The gene column is automatically added from rownames of each DEG file.
2. If the input files contain missing or invalid values, the function will throw warnings or errors.
3. File format must be .tsv or .csv, other formats (e.g., Excel .xlsx) are not supported directly.

2. plot_DEG_volcano()— Volcano Plot for Differential Expression Analysis

2.1 Purpose

The plot_DEG_volcano() function generates a volcano plot for visualizing differential expression results, highlighting significantly upregulated, downregulated, and stable genes, and also used to illustrate transcriptomic differences between experimental groups.

2.2 Input Data Format

This function requires a data frame (output from merge_DEG_datasets()) with the following mandatory columns:

Tip: Use merge_DEG_datasets() to generate this standardized input data frame.

2.3 Function Parameters

2.4 Output

1. Visual Output: A ggplot object showing genes grouped by contrast, with log2 fold change on the Y-axis and group on the X-axis. The color and size of each point reflects gene regulation status and significance.
2. File Output: A PDF file named volcano_plot.pdf will be saved to the specified output_path.

2.5 Run Example

# Prepare merged DEG data (e.g. from TCGA and mouse models)
merged_data <- merge_DEG_datasets(
  human_file_path = "./data/TCGA_DEG.tsv",
  mouse_files = list(
    "GSE172629" = "./data/GSE172629_DEG.tsv",
    "GSE208279" = "./data/GSE208279_DEG.tsv"
  )
)

# Generate volcano plot with default settings
p <- plot_DEG_volcano(deg_data = merged_data)

# Display the plot
p

# Further customization using ggplot2
p + ggplot2::theme_minimal()

2.6 Notes

1. Dot color is determined by the change column (e.g., “Up”, “Down”, “Stable”).
2. Transparency and dot size are proportional to -log10(FDR).

3. prepare_DEG_heatmap_data()— Process DEGs data for Heatmap Visualization

3.1 Purpose

The prepare_DEG_heatmap_data() function processes differential expression gene (DEG) results from TCGA and multiple mouse models to prepare a standardized dataset suitable for heatmap visualization. It performs significance filtering, homologous gene conversion, and logFC simplification, returning cleaned and comparable list of matrices for downstream visualization.

3.2 Input Data Format

See the file format requirements in the merge_DEG_datasets() documentation.

Each DEG file (both TCGA and mouse) must contain at least the following columns:

• logFC – log2 fold change
• FDR – false discovery rate (adjusted p-value)

Files must be in .tsv format and have gene identifiers as row names.

3.3 Function Parameters

3.4 What the Function Does

1. Load and Filter TCGA Data:

   • Reads the TCGA DEG file.

   • Selects genes with |logFC| ≥ 1 and FDR < 0.05.

2. Convert Human Genes to Mouse Homologs:

   • Uses homologene::homologene() to map human genes to mouse homologs.

   • Only genes with successful mappings are retained.

3. Count Up/Downregulated Genes:

   • Separates positive (logFC > 0) and negative (logFC < 0) genes.

4. Format TCGA logFC Values:

   • Upregulated: logFC set to 1

   • Downregulated: logFC set to -1

   • Others: set to NA

5. Process Mouse Model Files:

   • For each file:
     • Subsets rows to match TCGA homologs.
     • Retains only logFC and FDR columns.

3.5 Output

Returns a named list with the following components:

3.6 Run Example

# Define file paths
tcga_path <- "./data/TCGA_DEG.tsv"
mouse_files <- list(
  "Our_Model" = "./data/Our_Model_DEG.tsv",
  "GSE172629" = "./data/GSE172629_DEG.tsv"
)

# Run preprocessing
processed_data <- prepare_DEG_heatmap_data(tcga_path, mouse_files)

# View output
head(processed_data$processed_tcga)
head(processed_data$processed_mouse$Our_Model)

3.7 Note

1. The returned output is ready for use in the plot_DEG_heatmap() functions.
2. The logFC simplification logic ensures compatibility across datasets and makes interpretation of heatmaps more intuitive.

4. plot_DEG_heatmap()— Visualize Differential Expression Patterns in a Heatmap

4.1 Purpose

The plot_DEG_heatmap() function generates a binary expression heatmap (values of 1, -1, or NA) across multiple mouse models based on TCGA-derived significant genes. It is intended to highlight genes with high or low expression in TCGA tumors and examine their cross-species expression trends in mouse models.

This function is typically used after preparing input data with prepare_DEG_heatmap_data().

4.2 Input Data Format

1. This function expects the deg_data to be the output from prepare_DEG_heatmap_data().

4.3 Function Parameters

4.4 What the Function Does

1. Validates Input:
   • Ensures that deg_data contains the expected structure.
2. Combines Mouse Model Data:
   • Extracts logFC information from each mouse model dataset.
   • Combines them into a single binary matrix (merge_edger) using the TCGA gene set.
3. Assigns Row/Column Labels:
   • Gene symbols are used as rownames.
   • Column names are either taken from mouse_files or user-defined via col_names.
4. Plots Binary Expression Matrix:
   • 1 → gene highly expressed in TCGA tumor samples.
   • -1 → gene lowly expressed in TCGA tumor samples.
   • NA → not significant or not conserved.
   • Color-coded as red (up), blue (down), and white (NA).
   • Drawn using ComplexHeatmap::Heatmap() and saved to PDF.
5. Returns Processed Matrix:
   • Output matrix can be reused in downstream plots (e.g., plot_DEG_barplot()).

4.5 Output

1. A PDF heatmap saved to the provided output_path.
2. The function also returns the merged matrix (merge_edger) used for plot_DEG_barplot().

4.6 Color Legend

Legend labels:

• "Highly expression in tumor" = 1
• "Low expression in tumor" = -1

4.7 Run Example

# Process DEG data
processed_data <- prepare_DEG_heatmap_data(tcga_file, mouse_files)

# Generate heatmap
merge_edger <- plot_DEG_heatmap(
  deg_data = processed_data,
  mouse_files = mouse_files
)

4.8 Notes

1. This function is tightly coupled with prepare_DEG_heatmap_data() — don’t use it with arbitrary input.
2. You can reuse merge_edger for barplot visualization via plot_DEG_barplot().

5. plot_DEG_barplot()— Barplot of DEG Counts in Mouse Models

5.1 Purpose

The plot_DEG_barplot() function generates a stacked barplot to compare the number of upregulated and downregulated genes in each mouse model, based on TCGA-derived gene categories. It uses the DEG matrix generated by the plot_DEG_heatmap() function.

5.2 Input Data Format

• The function requires a matrix (merge_edger) as input, which must be the output of plot_DEG_heatmap() and contain:

  • Rows: Filtered genes shared with TCGA.

  • Columns: Mouse models.

  • Values:

    • 1 for genes upregulated in TCGA tumors.
    • -1 for downregulated genes.
    • NA for genes not meeting the threshold or unmatched.

5.3 Function Parameters

5.4 What the Function Does

1. Validates Input:
   • Checks that the input is a matrix with numeric values (1, -1, NA).
2. Transforms Matrix to Long Format:
   • Converts the wide-format DEG matrix into long format for ggplot2.
3. Assigns Group Labels:
   • Converts values to "up" (for 1) or "down" (for -1) gene categories.
   • NA values are excluded.
4. Generates Barplot:
   • X-axis: Mouse model names (contrast).
   • Y-axis: Count of up/down genes per model.
   • Colors: Customizable (default: red for up, blue for down).
   • Uses ggplot2::geom_bar() to generate a clean stacked barplot.
5. Saves Output:
   • Barplot is saved as "barplot.pdf" inside the specified output_path.

5.5 Output

5.6 Color Legend

You can change colors via the colors argument, for example:

colors = list(up = "red", down = "blue")

5.7 Run Example

# Generate barplot from merge_edger matrix
barplot <- plot_DEG_barplot(
  heatmap_data = merge_edger,
  colors = list(up = "#ff7676", down = "#66d4ff")
)

# Display in RStudio
print(barplot)

5.8 Notes

1. This function is intended to be used after plot_DEG_heatmap().
2. Genes not matched or not significantly regulated (NA) are automatically excluded.
3. The returned ggplot object allows additional customization using the full ggplot2 ecosystem.

6. prepare_PCA_data()- prepare data for PCA plotting

6.1 Purpose

The prepare_PCA_data() function merges TCGA and mouse model expression data, assigns group and batch labels, and optionally performs batch correction. This prepares the data for PCA analysis and plotting of transcriptional similarity in downstream function plot_PCA_tinyarray().

6.2 Input Data Format

Before you wanna using this fuction, you should first obtain the homologous gene-counts matrix files for TCGA and each mouse model, which could be obtained via the homologene package. The specific code and the expected data is as follows:

# For TCGA
count <- read.table("./TCGA.rawcounts")
library(homologene)
homolo_gene <- homologene(rownames(count), inTax = 9606, outTax = 10090)[,c(1:2)]
count <- count[rownames(count) %in% homolo_gene[, 1], ]
rownames(count) <- toupper(rownames(count))
count <- rownames_to_column(count,var = "X") # The var parameter must set as "X"


# For mouse model
count <- read.table("./GSE.rawcounts")
library(homologene)
homolo_gene <- homologene(rownames(count), inTax = 10090, outTax = 9606)[,c(1:2)]
count <- count[rownames(count) %in% homolo_gene[, 1], ]
rownames(count) <- toupper(rownames(count))
count <- rownames_to_column(count,var = "X") # The var parameter must set as "X"

The input expression files (for both TCGA and mouse models) must:

1. Be in .tsv or .csv format.
2. Contain homologene symbols as the first column (colnames = X).

6.3 Function Parameters

6.4 What the Function Does

1. Reads TCGA and mouse model expression files.
2. Merges all datasets using gene symbols.
3. Generates:
   • group_all: sample group labels like TCGA_normal, GSExxx_tumor.
   • batch: batch labels for each dataset.
4. Applies batch correction using sva::ComBat() if enabled.

6.5 Output

Returns a list containing:

6.6 Run Example

mouse_files <- list(
  "GSE172629" = "./data/GSE172629.homo.tsv",
  "GSExxxxxx" = "./data/GSExxxxxx.homo.csv"
)

sample_counts <- list(
  "TCGA" = c(normal = 50, tumor = 374),
  "GSE172629" = c(normal = 3, tumor = 3)
)

pca_data <- prepare_PCA_data(
  tcga_file = "./data/TCGA.tsv",
  mouse_files = mouse_files,
  sample_counts = sample_counts,
  batch_correction = TRUE
)

6.7 Notes

1. File formats are automatically detected (.tsv or .csv).
2. All gene sets across datasets must be aligned (same gene order after merge).
3. If sample counts do not match actual column numbers, an error is thrown.
4. The output corrected_data is ready for PCA, clustering, or heatmap visualization.

7. plot_PCA_tinyarray()- PCA visualization

7.1 Purpose

The plot_PCA_tinyarray() function performs principal component analysis (PCA) and visualizes the result using the tinyarray::draw_pca() function. It supports both raw and batch-corrected expression matrices and highlights sample grouping with customizable colors.

7.2 Input Data Format

The function expects a list object returned from prepare_PCA_data(), which includes:

1. A merged expression matrix (merged_data)
2. Batch labels (batch)
3. Sample group labels (group_all)

• An optional batch-corrected matrix (corrected_data)

7.3 Function Parameters

7.4 What the Function Does

1. Selects batch-corrected or raw merged data depending on the batch_correction flag.
2. Converts sample group vector into a factor.
3. Automatically assigns or validates colors for each sample group.
4. Performs PCA using tinyarray::draw_pca(), generating a 2D PCA plot with group labels and coloring.

7.5 Output

This function does not return an object, but directly draws the PCA plot to the active graphic device (RStudio Plots pane, PDF, etc.) for personalized adjustment.

7.6 Run Example

# Assume you have run `prepare_PCA_data()` and stored the result in `pca_data`

plot_PCA_tinyarray(
  pca_data = pca_data,
  batch_correction = TRUE  # Use ComBat-corrected data
)

# Customize colors manually (optional)
custom_colors <- c("#1b9e77", "#d95f02", "#7570b3")
plot_PCA_tinyarray(
  pca_data = pca_data,
  batch_correction = FALSE,
  colors = custom_colors
)

7.7 Notes

1. If batch_correction = TRUE, but no corrected_data exists in the input list, the function falls back to using the original merged data.
2. Ensure that group_all matches the number and order of samples in the expression matrix.
3. Default colors support up to 15 groups. Beyond that, manually define a longer colors vector.

8. prepare_similarity_data()- combine CNV and clinical file for similarity analysis

8.1 Purpose

This function prepares TCGA patient-level clinical and CNV (copy number variation) information for use in cross-species similarity analysis with mouse models. It extracts user-defined clinical risk factors(condition) and CNV data for selected genes and merges them into a clean, analysis-ready data frame.

8.2 Input Data Format

• The input clinical_file must be separated by tabs and contain the column named “PATIENT_ID”. And the character value in “PATIENT ID” should be unique!!! The table read using read.table(clinic_file_path, header = TRUE, sep = "\t", stringsAsFactors = FALSE) is presented in the figure.

• The input cnv_file should be a tab-separated dataframe with rows representing gene symbols and columns representing TCGA PATIENT_IDs. The table read using read.table (cnv_file_path, header=TRUE, sep="\ t", stringsAsFactors=False) is presented in the figure.

8.3 Function Parameters

8.4 What the Function Does

1. Reads the clinical and CNV data from TCGA .tsv files.
2. Extracts only the PATIENT_ID column and the user-specified clinical variables (via pattern matching).
3. Selects only CNV data rows corresponding to the user-specified cnv_genes.
4. Merges CNV and clinical variables into a unified matrix using PATIENT_ID as the key.
5. Encodes CNV values as ordered factors with levels -2, -1, 0, 1.

8.5 Output

A data frame containing:

• PATIENT_ID: Unique patient ID (standardized with . instead of -)
• One column per clinical_vars, with matched status (or NA if not present)
• One column per cnv_genes, representing CNV status (-2, -1, 0, or 1) as factors

This merged output is used as input in downstream similarity scoring functions (e.g., compute_similarity_matrix() → plot_TROMscore_heatmap()).

8.6 Run Example

merge_clinic_CNV <- prepare_similarity_data(
  clinic_file = "data_bcr_clinical_data_patient.tsv",
  cnv_file = "CNA.tsv",
  clinical_vars = c("Hepatitis B", "Non-Alcoholic"),
  cnv_genes = c("TP53", "PTEN"),
  risk_factor_column = "HISTORY_HEPATO_CARCINOMA_RISK_FACTORS"
)

8.7 Notes

1. The clinical_vars will be matched by partial string search (case-insensitive) in the column defined by risk_factor_column.
2. The resulting data frame is critical for annotating heatmaps and stratifying patients based on underlying genetic and clinical profiles.

9. compute_similarity_matrix()- Calculate similarity score

9.1 Purpose

The compute_similarity_matrix() function calculates a transcriptional similarity matrix between TCGA samples and mouse models, using the ws.trom() algorithm from the TROM package. It filters the similarity matrix to retain only biologically meaningful comparisons (tumor-tumor and normal-normal) and merges the filtered matrix with TCGA clinical and CNV data for downstream analysis.

9.2 Input Data Format

The input expression_matrix should be the batch-corrected expression matrix (typically the output from prepare_PCA_data()), where:

• Rows represent genes.
• Columns represent TCGA and mouse model samples.
• Row names are gene symbols.

The input merge_clinic_CNV is the processed clinical + CNV data result from prepare_similarity_data()

9.3 Function Parameters

9.4 What the Function Does

1. Calculates transcriptional similarity matrix using TROM::ws.trom() between TCGA and mouse models.
2. Filters out irrelevant comparisons:
   • Removes tumor-normal and normal-tumor comparisons.
   • Retains only tumor-tumor and normal-normal similarities.
3. Appends clinical annotations and CNV data to the filtered matrix.
4. Sorts the final matrix by tissue type and clinical/CNV variables.

9.5 Output

The function returns a named list with two components:

• similarity_matrix: The full similarity matrix (unfiltered), as computed by ws.trom().
• dm_trom_del_backgroud_clinic: The filtered similarity matrix (excluding background noise), merged with clinical and CNV data for use in plot_TROMscore_heatmap() and plot_TROMscore_boxplot().

9.6 Run Example

similarity_results <- compute_similarity_matrix(
  expression_matrix = pca_data$corrected_data,
  group_all = pca_data$group_all,
  tcga_sample_size = 424,
  tcga_normal_count = 50,
  merge_clinic_CNV = merge_clinic_CNV,
  cnv_genes = c("TP53", "PTEN"),
  clinical_vars = c("Hepatitis B", "Non-Alcoholic")
)

9.7 Notes

1. ws.trom() performs cross-species comparison based on gene expression similarity and returns a matrix of TROM scores.
2. The filtering step ensures that only biologically comparable samples (tumor vs. tumor, normal vs. normal) are retained.
3. The function requires accurate group_all labeling, where TCGA samples are listed first.
4. The PATIENT_ID is inferred from the first 12 characters of TCGA sample IDs.

10. plot_TROMscore_heatmap()- TROM similarity score heatmap visualization

10.1 Purpose

The plot_TROMscore_heatmap() function generates a heatmap visualization of transcriptional similarity scores (TROM scores) between TCGA samples and mouse cancer models. It provides a high-level overview of how closely mouse models resemble different human tumor subtypes at the transcriptomic level. It also integrates CNV and clinical information in TCGA samples.

10.2 Input Data Format

1. dm_trom_del_backgroud_clinic: The data frame from compute_similarity_matrix(), containing filtered similarity scores merged with clinical and CNV data.
2. trom_matrix: The data frame from compute_similarity_matrix(), containing the raw TROM similarity matrix before background filtering.
3. Sample groupings and counts must match those used in earlier PCA and similarity computations (e.g., from prepare_PCA_data() and compute_similarity_matrix()).

10.3 Function Parameters

10.4 What the Function Does

1. Extracts TROM scores for TCGA and mouse model samples.
2. Computes average TROM scores for each mouse model group (tumor/normal).
3. Formats the similarity matrix into a structure suitable for ComplexHeatmap.
4. Adds two layers of annotations:
   • Column annotations: whether each mouse model sample is tumor or normal.
   • Row annotations: TCGA clinical and CNV information.
5. Draws and saves the heatmap using the ComplexHeatmap package.

10.5 Output

The function returns a list with two elements:

1. heatmap_obj: The drawn heatmap object (from ComplexHeatmap::draw()).
2. merge_long: A long-format data frame containing:
   • x: Mouse model sample label.
   • y: Mean similarity score.
   • group: Tumor or normal.
   • condition: Experimental model condition.
   • score: Rounded average TROM score.

This merge_long object is used directly as input for plot_TROMscore_boxplot().

10.6 Run Example

heatmap_results <- plot_TROMscore_heatmap(
  dm_trom_del_backgroud_clinic = similarity_results$dm_trom_del_backgroud_clinic,
  trom_matrix = similarity_results$similarity_matrix,
  tcga_sample_size = 424,
  tcga_normal_count = 50,
  model_group_count = 52,
  condition = list("HBV pten p53 ko" = 5, "DEN Cl4" = 16),
  clinical_vars = c("Hepatitis B", "Non-Alcoholic"),
  cnv_genes = c("TP53", "PTEN"),
  group_all = pca_data$group_all,
  output_path = "./similarity/TROM_heatmap.pdf"
)

10.7 Notes

1. This heatmap emphasizes tumor-to-tumor and normal-to-normal transcriptional similarity across species.
2. Colors are mapped from low (white) to high (deep red) TROM scores using circlize::colorRamp2().
3. Clinical annotations are reversed (rev()) to align correctly with TCGA row order in the matrix.
4. The merge_long result can be reused for barplots or boxplots of similarity scores.
5. The function automatically creates directories if output_path does not exist.

11. plot_TROMscore_boxplot()- TROM similarity score boxplot visualization

11.1 Purpose

The plot_TROMscore_boxplot() function generates two types of plots to visualize transcriptional similarity scores (TROM scores):

1. A boxplot that compares TROM scores between tumor and normal groups across mouse models.
2. A tileplot that displays the average TROM score for each mouse model.

These visualizations help assess which mouse models best resemble specific TCGA tumor profiles.

11.2 Input Data Format

• The input merge_long should be the output from plot_TROMscore_heatmap(), a data frame in long format containing:
  • x: sample labels
  • y: TROM score
  • group: “tumor” or “normal”
  • condition: experimental condition/group
  • score: average score (used in tileplot)

11.3 Function Parameters

11.4 What the Function Does

1. Checks input: Ensures that merge_long contains the required columns.
2. Generates a boxplot: Displays similarity score distributions per condition (tumor/normal) using ggplot2::geom_boxplot().
3. Computes average scores per condition and generates a tileplot using ggplot2::geom_tile().
4. Combines plots vertically using patchwork::plot_layout() if combine_plots = TRUE.
5. Saves the final plot to the specified output directory as a .pdf.

11.5 Output

• If combine_plots = TRUE, returns a single combined ggplot object.
• If combine_plots = FALSE, returns a list of two ggplot objects: list(tileplot, boxplot).
• A PDF file named TROM_boxplot.pdf is saved to the output_path.

11.6 Run Example

plot_TROMscore_boxplot(
  merge_long = heatmap_results$merge_long,
  output_path = "./similarity/TROM_boxplot",
  combine_plots = TRUE
)

11.7 Notes

1. The tileplot highlights model-level mean similarity scores using a color gradient (blue → white → red).
2. The boxplot separates tumor/normal groups per condition for distributional comparison.
3. A vertical dashed line (xintercept = 2.6) is included for visual reference (e.g., model split).
4. You may modify the returned ggplot object(s) using + theme() or other ggplot2 functions.

12. GO_enrichment_analysis()

12.1 Purpose

This function performs Gene Ontology (GO) enrichment analysis on differentially expressed genes (DEGs) identified from edgeR results. It supports analysis across species (human, mouse, rat) and GO categories (BP, MF, CC). The enriched GO terms are visualized using GO term similarity clustering via simplifyGO.

12.2 Input Data Format

Please refer to the input data structure described in the merge_DEG_datasets() documentation.

12.3 Function Parameters

12.4 What the Function Does

1. Reads the DEG file and filters genes with |logFC| > 1 and FDR < 0.05.
2. Maps gene symbols to gene IDs using the specified OrgDb package (e.g., org.Hs.eg.db).
3. Performs GO enrichment using clusterProfiler::enrichGO().
4. Calculates GO term similarity using simplifyEnrichment::GO_similarity().
5. Saves a GO similarity heatmap to a PDF using simplifyGO().

12.5 Output

• A .pdf file showing GO term similarity clustering, saved to output_path.

12.6 Run Example

GO_enrichment_analysis(
  deg_file = "./data/TCGA_DEG.tsv",
  species = "hsa",
  org_db = "org.Hs.eg.db",
  ont = "BP",
  column_title = "TCGA_GO_BP_terms"
)

12.7 Notes

1. Be sure to install and load the corresponding organism annotation package (org.Hs.eg.db, org.Mm.eg.db, or org.Rn.eg.db) before running the function.
2. Input DEG file must contain row names as gene symbols, and columns logFC and FDR.
3. The output file will be overwritten if one with the same name already exists in the output path.

13. plot_TME_barplot()- Tumor Microenvironment cell composition visualization

13.1 Purpose

This function performs immune infiltration analysis from differential expression gene matrix using various deconvolution methods (e.g., CIBERSORT, xCell), and visualizes the relative abundance of immune cell types across experimental groups in a stacked barplot.

13.2 Input Data Format

Please refer to the input format described in the documentation for merge_DEG_datasets(). The input can be either:

• A raw count matrix (genes × samples), or
• A TPM matrix (if input_type = "TPM"), with matching gene identifiers.

13.3 Function Parameters

13.4 What the Function Does

1. Checks input validity and converts counts to TPM using IOBR::count2tpm() (if needed).
2. Runs immune deconvolution via IOBR::deconvo_tme() using the specified method.
3. Aggregates the infiltration results by sample group.
4. Generates a stacked barplot using IOBR::cell_bar_plot().

13.5 Output

• PDF Figure: A stacked barplot (method_cell_bar_plot.pdf) visualizing immune cell compositions.
• CSV Table: Cell abundance table (method_result.csv) for all samples.
• Returned List:
  • eset_tpm: TPM-normalized expression matrix (used in analysis)
  • method_model: Full immune infiltration results (before group merging)

13.6 Run Example

group_as <- c(
  rep("TCGA_normal", 50),
  rep("TCGA_tumor", 374),
  rep("HBV_Pten_KO_normal", 3),
  rep("HBV_Pten_KO_tumor", 3)
)

TME_barplot_result <- plot_TME_barplot(
  data_matrix = merge.data,
  input_type = "count",
  method = "cibersort",
  group_as = group_as
)

# Extract TPM for downstream analysis
eset_tpm <- TME_barplot_result$eset_tpm

13.7 Notes

1. Ensure group_as matches the sample order in data_matrix.
2. Supported species for TPM conversion (count2tpm) is currently "hsa" only. Modify org = parameter if needed.
3. Deconvolution methods depend on package IOBR; make sure it’s properly installed and loaded.
4. Output will be written to output_path as .pdf and .csv.

14. calculate_metabolism_score()

14.1 Purpose

This function calculates metabolism-related gene signature scores from TPM-normalized RNA-seq expression data.

14.2 Input Data Format

Same as described in the plot_TME_barplot() documentation. The input must be a TPM matrix (genes × samples) with gene identifiers as row names and sample names as column names. Typically generated from plot_TME_barplot().

14.3 Function Parameters

14.4 What the Function Does

1. Validates input format and parameter values.
2. Calls IOBR::calculate_sig_score() using a predefined metabolism signature ("signature_metabolism").
3. Computes metabolism signature scores per sample using the selected scoring method.

14.5 Output

• A data frame where:
  • Rows represent samples.
  • Columns contain the metabolism score (or multiple scores, depending on method).

This can be used for downstream comparison, group visualization, or correlation analysis.

14.6 Run Example

# Assume eset_tpm is output from plot_TME_barplot()
sig_meta <- calculate_metabolism_score(
  eset_tpm = eset_tpm,
  method = "pca",
  mini_gene_count = 2
)

14.7 Notes

1. This function uses a fixed internal metabolism gene signature from IOBR.
2. PCA-based scores reflect the first principal component of metabolism-related gene expression.
3. If too few genes are available for scoring (mini_gene_count not met), those samples may be excluded.
4. Make sure IOBR is properly installed and loaded before running this function.

15. plot_metabolism_heatmap()- metabolism score visualization

15.1 Purpose

This function generates a heatmap visualization of metabolism-related signature scores across sample groups.

15.2 Input Data Format

The input sig_meta must be a data frame of metabolism signature scores, usually the output from calculate_metabolism_score(). It must include a first column identifying group labels (e.g., group_as).

15.3 Function Parameters

15.4 What the Function Does

1. Aggregates metabolism scores by group using group_as (must be defined globally).
2. Transposes the matrix so that pathways are rows and sample groups are columns.
3. Uses ComplexHeatmap::Heatmap() to create a heatmap.
4. Applies custom clustering method to both rows and columns.
5. Saves the plot as a PDF file in the specified output_path.

15.5 Output

• A heatmap PDF file named heatmap.pdf saved to output_path.

15.6 Run Example

# Assuming `sig_meta` is the output of `calculate_metabolism_score()`
heatmap_obj <- plot_metabolism_heatmap(
  sig_meta = sig_meta,
  output_path = "./metabolism/",
  clustering_method = "canberra",
  width = 10,
  height = 22
)

15.7 Notes

1. If group_as is not found or doesn’t match the sample structure, the function will fail silently or return incorrect results.
2. Distance metric selection (manhattan vs canberra) can influence clustering results; try both to explore structure.
3. Heatmap is saved automatically; no need to call ggsave() or other output functions.

16. plot_TME_boxplot()- immune infiltration scores boxplot visualization

16.1 Purpose

This function computes immune infiltration scores using the CIBERSORT method and visualizes them via boxplots for each immune signature. It is primarily used to compare immune cell composition between groups across different mouse models or experimental conditions.

16.2 Input Data Format

• TPM Expression Matrix A data.frame of gene expression values (TPM), genes × samples format. Usually generated by plot_TME_barplot().
• Group Labels (group_as) A character vector of group labels in the format "condition_group" (e.g., "HBV_normal", "HBV_tumor", "DEN_normal").
• Condition List (condition_list) A named list specifying the number of samples per mouse model:

condition_list <- list("GSEXXX" = 6, "GSExxx" = 8)

16.3 Function Parameters

16.4 What the Function Does

1. Runs CIBERSORT via IOBR::deconvo_tme() to get immune cell fractions.
2. Iterates through each immune signature (e.g., B cells, macrophages, etc.).
3. Generates ggplot boxplots showing infiltration score distribution across groups.
4. Saves each immune signature boxplot as a PDF file in the specified output path.

16.5 Output

• A series of PDF files (one per immune signature) saved to output_path.

16.6 Run Example

# Define sample count per condition
condition_list <- list("GSEXXX" = 6, "GSExxx" = 8)

# Run boxplot function
plot_list <- plot_TME_boxplot(
  eset_tpm = TME_barplot_result$eset_tpm,
  group_as = group_as,
  condition_list = condition_list
)

16.7 Notes

• The function currently supports only the CIBERSORT method by default.
• You can customize colors and themes by post-editing the ggplot objects returned.

Visualization Outputs

All plots generated are publication-ready, saved in PDF format by default, and include:

1. DEG Volcano plots, Heatmaps, Barplots
2. PCA Plots
3. TROM Score Heatmaps, Boxplots, Barplots
4. Immune Infiltration Barplots, Boxplots
5. GO Enrichment Similarity Plots
6. Metabolism Heatmap

Contributions & Issues

We warmly welcome contributions and suggestions! Please feel free to open an issue or pull request.

Contact

• Author: Dian Wang, Guiwen Guan
• Email: wangdian@bjmu.edu.cn
• GitHub:MurineMatchTx

Thanks for very much for using our package📦!