Computational Protocol To Quantify Artesunate-Associated Transcriptomic Shifts In Experimental Cerebral Malaria Brain Tissue Using DESeq2

Summary

Experimental cerebral malaria (ECM) causes neuroinflammation and blood–brain barrier dysfunction. Using GSE162535, this reproducible R-based RNA-seq workflow compares control (CB), ECM (MB), and artesunate-treated (AB) brains, performing QC, PCA, and DESeq2 analyses. It identifies ECM-driven transcriptional changes and evaluates artesunate-mediated modulation of inflammatory and neurovascular pathways.

Abstract

RNA sequencing (RNA-seq) is widely used to define disease-associated transcriptional programs, but consistent, end-to-end workflows are required to ensure reproducible comparisons across experimental groups and to facilitate biologically interpretable outputs. Here, a complete RNA-seq analysis protocol is presented for evaluating brain transcriptomic changes in experimental cerebral malaria (ECM) and following artesunate treatment using the public dataset GSE162535. The workflow analyzes three groups—control brain (CB), ECM brain (MB), and artesunate-treated ECM brain (AB)—starting from an HTSeq-count matrix. After importing and formatting raw counts, the protocol constructs a DESeq2 dataset with CB as the reference level, filters low-count genes, and performs quality control using library-size visualization, principal component analysis, and sample distance clustering. Differential expression is then calculated for three primary contrasts (MB vs CB, AB vs MB, AB vs CB), with log2 fold-change shrinkage applied for stable effect-size estimation. The protocol exports full and significant differential-expression tables, generates volcano plots and MA plots, and produces heatmaps of the most variable genes and top-ranked differentially expressed genes per contrast. To support mechanistic interpretation, the workflow includes targeted immune marker extraction (e.g., cytokines, chemokines, microglial activation markers, BBB/endothelial genes) and conducts Gene Ontology and KEGG enrichment analyses for significantly up- and downregulated gene sets. This protocol provides a reproducible template for characterizing ECM-associated neuroinflammatory transcriptional programs and for quantifying artesunate-associated modulation of these signatures.

Introduction

Cerebral malaria (CM) is a life-threatening neurological complication of Plasmodium falciparum infection and remains a major contributor to malaria mortality despite advances in case management. CM is characterized by acute encephalopathy, microvascular dysfunction, endothelial activation, and blood–brain barrier (BBB) disruption, with downstream neuroinflammation that can drive coma and, among survivors, persistent neurocognitive sequelae¹. The pathogenesis of CM is multifactorial and involves host inflammatory responses interacting with parasite- and host-derived factors at the neurovascular interface, making it difficult to infer causal mechanisms from clinical endpoints alone¹.

Experimental cerebral malaria (ECM) models, particularly Plasmodium berghei ANKA infection in C57BL/6 mice, provide a tractable platform to interrogate brain-specific immunopathology, BBB injury, and neuroinflammatory signaling in a controlled setting^2,3. These models have been used to map cellular and molecular responses across disease stages and to test adjunctive interventions in vivo^2,3. However, ECM pathobiology is complex and highly dynamic, and targeted assays may miss coordinated, pathway-level shifts that occur across multiple immune and neurovascular programs.

Artesunate is the recommended first-line parenteral therapy for severe malaria and has demonstrated a substantial survival benefit compared with quinine across key evidence bases⁴. Although rapid parasite clearance is central to artesunate efficacy, neurological outcomes likely reflect both parasite reduction and secondary modulation of inflammatory and neurovascular pathways^1,4. Understanding how artemisinin treatment reshapes brain transcriptional programs during ECM can therefore provide mechanistic insights that complement clinical efficacy data and may identify candidate pathways for adjunctive neuroprotective strategies.

RNA sequencing (RNA-seq) enables unbiased, genome-wide profiling of transcriptional responses in disease and treatment states, supporting differential expression analysis and downstream functional interpretation. Public repositories such as the NCBI Gene Expression Omnibus (GEO) provide curated datasets suitable for reproducible re-analysis, including GSE162535, which contains brain RNA-seq from control brains (CB), ECM brains (MB), and artesunate-treated ECM brains (AB)⁵. To support reproducible discovery from such datasets, robust statistical frameworks are required for count-based differential expression, and enrichment tools are needed to interpret gene-level changes in terms of biological pathways and processes.

This study presents a reproducible, end-to-end RNA-seq analysis workflow for brain tissue across control (CB), experimental cerebral malaria (MB), and artesunate-treated (AB) groups. The novelty of this article lies in its standardized DESeq2-based pipeline, incorporating predefined biologically relevant contrasts (MB vs CB, AB vs MB, and AB vs CB), rigorous quality control outputs (library size assessment, principal component analysis, and sample distance heatmaps), and integrated downstream interpretation through Gene Ontology (GO) and KEGG pathway enrichment using clusterProfiler^6,7. In addition, the workflow implements structured immune-panel–based interpretation, enabling systematic characterization of neuroinflammatory, immune, and neurovascular transcriptional responses. By combining statistical rigor, transparency, and publication-ready outputs, this protocol provides a robust and reusable framework for analyzing ECM-associated transcriptomic dysregulation and assessing treatment-driven modulation in preclinical malaria studies.

Protocol

This study utilized publicly available RNA-seq data and did not involve any new human or animal experiments. Therefore, ethical approval and informed consent were not required (Table of Materials).

1. Prepare the computing environment and folder structure

1. Set up hardware and operating system
   1. Use a workstation or laptop with at least 8 GB RAM (16 GB recommended) and ≥10 GB free disk space for downloads and outputs.
   2. Use Windows, macOS, or Linux with permission to install R packages and write files to the working directory.
2. Install required software
   1. Install R (version 4.2 or later). Install RStudio Desktop (recommended) to run the workflow interactively and to manage the project directory.
3. Create a project directory and output folders
   1. Create a new folder for the analysis (e.g., GSE162535_RNAseq_DESeq2). Set this folder as the R working directory.
   2. Create the output folders exactly as follows:
      1. Create results/. Create results/fig/. Create results/tables/.
4. Install required R packages
   1. Install CRAN packages: tidyverse, pheatmap, and RColorBrewer. Install Bioconductor packages: DESeq2, apeglm (optional), clusterProfiler, and org.Mm.eg.db.
   2. Load the required packages at the start of the script: DESeq2, tidyverse, pheatmap, RColorBrewer, clusterProfiler, org.Mm.eg.db, and ggplot2. Set a seed for reproducibility by running set.seed(123).
5. Apply data integrity and privacy safeguards
   1. Store downloaded files in the project directory only. Restrict access to the project directory if the environment contains sensitive or embargoed data. Record software versions by exporting sessionInfo() to results/sessionInfo.txt at the end of the workflow.
      NOTE: This protocol performs an in silico analysis of public RNA-seq count data and does not involve handling biological specimens.

2. Obtain the RNA-seq count matrix and define the experimental groups

1. Download the GEO dataset files for GSE162535
   1. Download the HTSeq count file for GSE162535 from GEO and save it in the project directory as GSE162535_All.HTSeq.counts.txt.gz. Verify the file integrity by confirming that the file opens without error using gzfile() and read.delim().
2. Load the HTSeq counts into R
   1. Import the compressed count table using read.delim(gzfile(...), header = TRUE, check.names = FALSE, quote = "", comment.char = ""). Inspect the imported object using str(), head(), and colnames() to confirm that:
      1. The first column contains gene identifiers (e.g., AccID). The remaining columns contain sample-level counts.
3. Standardize data types and resolve duplicated gene identifiers
   1. Define sample columns as all columns except the gene identifier column (AccID). Coerce all sample columns to integer counts. Collapse duplicated gene identifiers by summing counts across rows sharing the same AccID.
   2. Convert the collapsed table to a standard data frame. Set row names to the gene identifier and remove the identifier column from the count matrix.
4. Verify the count matrix structure
   1. Confirm that the count matrix has 12 sample columns. Confirm that the sample column names follow the format AB_1..AB_4, CB_1..CB_4, and MB_1.MB_4.
   2. Stop execution if the count matrix does not have 12 columns.
5. Create and validate sample metadata
   1. Create a sample metadata table (colData) containing the following columns:
      1. sample: unique sample identifiers corresponding to the column names of the count matrix. group: biological condition assigned to each sample.
   2. Retrieve the corresponding sample annotations for dataset GSE162535 from the Gene Expression Omnibus (GEO).
   3. Verify each sample identity using GEO accession information and sample annotation fields describing experimental condition. Map each sample identifier in the count matrix to its verified GEO annotation.
   4. Assign each verified sample to one of the following groups: control brain (CB), experimental cerebral malaria brain (MB), or artesunate-treated experimental cerebral malaria brain (AB).
   5. Arrange the metadata table so that the sample order matches the column order of the count matrix.
   6. Compare sample identifiers in the metadata table against the count matrix to confirm one-to-one correspondence. Encode the group variable as a factor with levels ordered as CB, MB, and AB.
   7. Set CB as the reference level to define the control brain group as the baseline for downstream differential expression analysis.
   8. Set the row names of the metadata table equal to the sample identifiers.
      NOTE:Group labels were not assigned solely from the order of samples in the count matrix. Sample identities were independently verified using GEO metadata and accession information before alignment with the expression matrix to improve reproducibility and reduce the risk of sample misclassification.

3. Construct the DESeq2 dataset and perform basic quality control

1. Sanitize the count matrix
   1. Replace missing values in the count matrix with zero. Confirm the absence of missing values using column-wise is.na() summaries.
2. Create the DESeq2 dataset
   1. Create a DESeqDataSet using DESeqDataSetFromMatrix() with: countData = counts; colData = sample_info; design = ~ group.
3. Filter low-count genes
   1. Remove genes with total counts <10 across all samples using dds <- dds[rowSums(counts(dds)) >= 10, ]. Record the number of retained genes by printing the object summary.
4. Evaluate library sizes
   1. Compute library sizes as the column sums of the filtered count matrix. Generate a library size barplot and save it as results/fig/library_sizes.pdf.
   2. Inspect the library size plot and confirm that no sample has an extreme sequencing depth inconsistent with the study design.

4. Run DESeq2 and generate transformation objects for visualization

1. Fit the DESeq2 model
   1. Run differential expression modeling using dds <- DESeq(dds). Retain the fitted DDS object for all downstream results extraction.
2. Create transformed expression matrices
   1. Compute the regularized log transformation using rld <- rlog(dds, blind = FALSE). Compute the variance-stabilizing transformation using vsd <- vst(dds, blind = FALSE). Use rld and vsd for PCA, clustering, and heatmaps.
      NOTE: Use blind = FALSE to preserve group-dependent variance structure.

5. Perform global QC using PCA and sample-to-sample distance clustering

1. Generate principal component analysis (PCA)
   1. Compute PCA using plotPCA(rld, intgroup = "group", returnData = TRUE). Extract percent variance explained for PC1 and PC2. Plot PC1 vs PC2 using ggplot2, labeling points by sample name and coloring by group.
   2. Save the PCA plot as results/fig/PCA_samples.pdf. Confirm that biological replicates cluster by group and that no sample separates as an outlier.
2. Generate a sample distance heatmap
   1. Compute pairwise sample distances using dist(t(assay(vsd))). Convert the distance object to a matrix for visualization. Create a column annotation table containing the group factor for each sample.
   2. Plot the distance matrix using pheatmap() and save as results/fig/sample_distance_heatmap.pdf. Confirm that samples cluster primarily by group.

6. Compute the differential expression for the three primary contrasts

1. Define the contrasts
   1. Define the disease effect as ECM vs control: MB vs CB. Define the treatment effect within ECM as artesunate-treated vs ECM: AB vs MB. Define treatment vs baseline as artesunate-treated vs control: AB vs CB.
2. Extract DESeq2 results with log2 fold-change shrinkage
   1. Extract raw results for each contrast using results(dds, contrast = c("group", groupA, groupB)). Shrink log2 fold changes using lfcShrink(dds, contrast = c("group", groupA, groupB), res = res, type = "normal").
   2. Convert results to a data frame and store the gene identifier as a column named gene_id. Order results by nominal p-value for stable reporting.
   3. Save each full results table to results/tables/ as:
      1. DESeq2_MB_vs_CB_all_genes.csv, DESeq2_AB_vs_MB_all_genes.csv, DESeq2_AB_vs_CB_all_genes.csv
         NOTE: If shrinkage fails due to package configuration, rerun with an alternative shrinkage type supported in the local installation.
3. Define significance thresholds and export significant gene sets
   1. Define differentially expressed genes (DEGs) using: adjusted p-value (FDR) < 0.05 and absolute log2 fold change ≥ 1. Filter each contrast to exclude genes with missing adjusted p-values.
   2. Export significant DEG tables to: results/tables/DESeq2_MB_vs_CB_sig.csv, results/tables/DESeq2_AB_vs_MB_sig.csv, results/tables/DESeq2_AB_vs_CB_sig.csv
4. Summarize DEG counts per contrast
   1. Compute the number of significant genes per contrast. Save the summary table as results/tables/DE_summary_counts.csv.

7. Generate volcano plots for each contrast.

1. Create a volcano plot function.
   1. Compute -log10(adjusted p-value) for each gene. Classify each gene as Up, Down, or Not Significant using the thresholds: FDR < 0.05 and |log2FC| ≥ 1.
   2. Plot log2FC (x-axis) vs -log10(FDR) (y-axis) using ggplot2. Add dashed threshold lines at log2FC = ±1 and -log10(0.05).
2. Export volcano plots
   1. Save each plot as a PDF in results/fig/: volcano_MB_vs_CB.pdf, volcano_AB_vs_MB.pdf, volcano_AB_vs_CB.pdf.
      NOTE: Use consistent axis limits across contrasts to support visual comparison across figures.

8. Generate MA plots for each contrast

1. Export MA plots to a single PDF.
   1. Open a PDF device named results/fig/MA_plots.pdf. Plot MA plots for each raw DESeq2 result object using plotMA(). Label each plot with the contrast name. Close the PDF device.

9. Generate expression heatmaps for global variability and contrast-specific DE genes.

1. Plot the top variable genes across all samples.
   1. Compute row-wise variance across samples from the vsd assay matrix. Select the top 100 most variable genes. Mean-center each gene across samples.
   2. Generate a heatmap using pheatmap() with sample group annotations. Save the figure as results/fig/heatmap_top100_variable_genes.pdf.
2. Plot the top differentially expressed genes for each contrast.
   1. Select the top 50 genes by adjusted p-value for each contrast. Extract their VSD expression matrix and mean-center by gene. Generate a heatmap for each contrast and save as:
   2. results/fig/heatmap_top50_MB_vs_CB.pdf, results/fig/heatmap_top50_AB_vs_MB.pdf, results/fig/heatmap_top50_AB_vs_CB.pdf
      NOTE: Increase PDF width and height if row labels are enabled.

10. Perform targeted immune marker analysis.

1. Summarize selected immune genes across contrasts.
   1. Define an immune marker set: Il6, Il1b, Il10, Tnf, Ifng, Il21, and Icam1. Extract log2FC, p-value, and adjusted p-value for these genes from each contrast result table.
   2. Merge the three contrast summaries by gene identifier. Save the merged table as results/tables/immune_genes_summary.csv.
2. Generate an immune marker heatmap (rlog).
   1. Identify immune genes present in the transformed assay matrix. Extract the rlog expression matrix for present immune genes.
   2. Mean-center expression by gene. Generate and save the heatmap to results/fig/heatmap_immune_genes.pdf.

11. Perform extended immune panel analysis by functional category

1. Define immune panels.
   1. Define immune marker panels by category, including: Pro-inflammatory cytokines and interferon response genes, anti-inflammatory and regulatory genes, chemokines, microglial activation markers, astrocyte markers, BBB and endothelial activation markers, T-cell markers and exhaustion markers, monocyte/macrophage markers, complement pathway genes, and oxidative stress and cell death genes.
2. Extract the immune panel DE results for each contrast.
   1. Create a mapping table from gene symbols to immune categories. Filter each contrast result table for the immune panel genes. Join the category mapping to each filtered result table.
   2. Concatenate the three contrast-specific immune tables into a single table. Save the combined table as results/tables/immune_panels_DE_all_contrasts.csv. Generate a classified immune panel expression heatmap.
   3. Identify immune panel genes present in the rld assay matrix. Append category labels to row names to preserve panel classification in the figure. Plot the rlog heatmap with sample group annotations.
3. Save outputs.
   1. Save as: results/fig/heatmap_immune_panels_all.pdf, results/fig/heatmap_immune_panels_all.png

12. Generate immune-overlay volcano plots

1. Create immune-overlay volcano plots.
   1. Mark genes as immune or non-immune by joining the immune panel mapping to each contrast result table. Plot all non-immune genes as grey background points. Plot immune genes colored by immune category in the foreground. Add threshold lines at log2FC = ±1 and -log10(0.05).
2. Export immune-overlay volcano plots.
   1. Save MB vs CB immune-overlay volcano plots as:
      1. results/fig/volcano_MB_vs_CB_immune_overlay.pdf
      2. results/fig/volcano_MB_vs_CB_immune_overlay.png
   2. Save AB vs MB immune-overlay volcano plots as:
      1. results/fig/volcano_AB_vs_MB_immune_overlay.pdf
      2. results/fig/volcano_AB_vs_MB_immune_overlay.png
         NOTE: Curated immune gene panels were predefined and grouped into functional categories; the complete list of gene symbols with annotations is provided in Supplementary Table S1. This table ensures transparency and reproducibility, allowing direct reuse and validation of panel-based analyses across studies.

13. Perform functional enrichment analysis (GO and KEGG)

1. Define gene sets for enrichment
   1. For each contrast (MB vs CB and AB vs MB), extract significantly differentially expressed genes (DEGs) using an adjusted p-value threshold of < 0.05 and an absolute log₂ fold-change threshold ≥ 1.
   2. Split DEGs into upregulated genes (log₂FC > 0) and downregulated genes (log₂FC < 0).
2. Define the background (universe) gene set
   1. Use all genes retained after count filtering in the DESeq2 dataset as the background universe. Extract gene symbols from the filtered DESeq2 object. Convert background gene symbols to Entrez identifiers using the bitr() function with org.Mm.eg.db.
      NOTE: The use of a consistent background gene set ensures unbiased enrichment results.
3. Map gene identifiers for enrichment
   1. Convert gene symbols of each DEG set to Entrez identifiers using bitr(). Retain only successfully mapped genes for KEGG enrichment. Skip enrichment analysis for a gene set if no genes are successfully mapped.
      NOTE: Gene Ontology (GO) enrichment uses gene symbols, whereas KEGG enrichment requires Entrez identifiers.
4. Perform gene ontology (GO) enrichment (Biological Process)
   1. Run GO enrichment using enrichGO() with OrgDb = org.Mm.eg.db. Set keyType = "SYMBOL" and ontology (ont) = "BP". Specify the universe as all filtered gene symbols from the DESeq2 dataset.
   2. Use pAdjustMethod = "BH" for multiple-testing correction. Apply significance thresholds of pvalueCutoff = 0.05 and qvalueCutoff = 0.05. Restrict gene set sizes using minGSSize = 10 and maxGSSize = 500.
   3. Export GO enrichment results to the results/tables/ directory. Generate barplots of top enriched GO terms and save to results/fig/.
5. Perform KEGG pathway enrichment
   1. Run KEGG enrichment using enrichKEGG() with organism = "mmu". Provide Entrez identifiers of DEG sets as input genes. Use the mapped Entrez background as the universe.
   2. Apply pAdjustMethod = "BH". Use significance thresholds of pvalueCutoff = 0.05 and qvalueCutoff = 0.05. Restrict gene set sizes using minGSSize = 10 and maxGSSize = 500.
   3. Export KEGG enrichment tables to results/tables/. Generate barplots of enriched KEGG pathways and save to results/fig/.
      NOTE: KEGG enrichment may not return results if gene mapping is insufficient; such cases are handled without interrupting the workflow.

14. Perform immune-only enrichment analysis (optional module)

1. Construct immune-specific gene sets
   1. Define immune panel genes based on curated functional categories. Intersect DEG lists with immune panel genes. Split immune-specific DEGs into upregulated and downregulated sets for each contrast.
2. Perform enrichment on immune-specific gene sets.
   1. Apply the same GO and KEGG enrichment procedures described in Section 13. Use identical background universe and parameter settings. Save outputs using filenames containing the label “immune” to distinguish them from global enrichment results.

15. Generate enrichment dot plots (optional module)

1. Load enrichment results.
   1. Import GO or KEGG enrichment tables from the results/tables/ directory.
2. Transform enrichment metrics.
   1. Convert GeneRatio values from fractional format (x/y) to numeric ratios. Compute −log₁₀(adjusted p-values) for visualization.
3. Generate dot plots.
   1. Plot gene ratio on the x-axis and enriched term descriptions on the y-axis. Scale dot size by gene count and color by −log₁₀(adjusted p-value).
   2. Select top enriched terms based on adjusted p-value ranking. Save dot plots as PDF and PNG files in the results/fig/ directory.

16. Save session information and finalize the run

1. Export session information.
   1. Save sessionInfo() output to results/sessionInfo.txt to document R version and package versions.
2. Confirm successful completion
   1. Confirm that the workflow generated: Figures in results/fig/, tables in results/tables/.
   2. Confirm that the three primary contrast outputs exist and are non-empty: DESeq2_MB_vs_CB_all_genes.csv, DESeq2_AB_vs_MB_all_genes.csv, DESeq2_AB_vs_CB_all_genes.csv

Results

Data quality assessment and global transcriptomic structure

RNA-seq data from 12 brain samples (CB, MB, AB; n = 4 per group) were processed using the standardized workflow. After filtering low-count genes (≥10 total counts), the dataset was retained for downstream analyses. Principal component analysis (PCA) of rlog-transformed counts demonstrated separation of samples by experimental group, with replicates clustering consistently within groups (Figure 1A<...

Discussion

This study presents a reproducible, DESeq2-based workflow for the analysis of bulk RNA-seq data across predefined biological contrasts. The protocol integrates standardized preprocessing, verified metadata assignment, normalization, differential expression testing, and structured downstream outputs, enabling consistent and transparent generation of transcriptomic results. By explicitly defining analytical steps and parameters, the workflow is designed to improve reproducibility and comparability across studies

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
clusterProfiler (R package)	Bioconductor	RRID:SCR_016884	Functional enrichment analysis (GO and KEGG pathways)
DESeq2 (R package)	Bioconductor	RRID:SCR_015687	Differential expression analysis of count-based RNA-seq data
ggplot2 (R package)	CRAN	RRID:SCR_014601	Visualization of PCA plots, volcano plots, and summary figures
GitHub (optional)	GitHub Inc.	RRID:SCR_002630	Version control and sharing of reproducible scripts
HTSeq-count RNA-seq dataset (GSE162535)	NCBI Gene Expression Omnibus (GEO)	RRID:SCR_005012	Bulk RNA-seq count matrix used as input for analysis
matrixStats (R package)	CRAN	RRID:SCR_016361	Efficient computation of row/column statistics (e.g., variance)
openxlsx (R package)	CRAN	RRID:SCR_019215	Export of results tables to Excel format
Operating system	Microsoft / Apple / Linux	N/A	Windows 10+, macOS, or Linux supported
org.Mm.eg.db (R package)	Bioconductor	RRID:SCR_002815	Mouse gene annotation database for gene ID mapping
PDF viewer	Any	N/A	Viewing output figures (PCA, heatmaps, volcano plots)
Personal computer or workstation	Any	N/A	Minimum 16 GB RAM recommended for RNA-seq analysis
pheatmap (R package)	CRAN	RRID:SCR_016418	Heatmap visualization of gene expression and clustering
R Statistical Software (version ≥ 4.2)	R Foundation for Statistical Computing	RRID:SCR_001905	Core computational environment for all RNA-seq analysis
RColorBrewer (R package)	CRAN	RRID:SCR_015742	Color palettes for heatmaps and plots
RStudio Desktop	Posit Software	RRID:SCR_000432	Integrated development environment (IDE) for scripting and reproducibility
Sample metadata file (CSV format)	Generated / GEO annotations	N/A	Curated sample annotation linking samples to CB, MB, and AB groups
stringr (R package)	CRAN	RRID:SCR_019195	String processing for enrichment visualization and formatting
tibble (R package)	CRAN	RRID:SCR_019186	Data frame handling and tidy data structures
tidyverse (R package suite)	CRAN	RRID:SCR_019186	Data manipulation, transformation, and visualization