The short-chain dehydrogenase gene DHRS11 is involved in cantharidin synthesis in the blister beetle Epicauta impressicornis

This repository contains the complete analysis pipeline for weighted gene co-expression network analysis (WGCNA), differential expression analysis, and integrated gene–metabolite association analysis of the blister beetle Epicauta impressicornis (Coleoptera: Meloidae) RNA-seq and untargeted metabolomics data. The study focuses on the mevalonate (MVA) pathway and cantharidin biosynthesis, and identifies the short-chain dehydrogenase gene DHRS11 as a candidate involved in cantharidin synthesis.

Pipeline overview

01.data_preprocessing.R     →  gene_raw_counts.tsv + gene_all_fpkm.tsv
02.DESeq2_analysis.R        →  DEG lists + normalized counts
03.WGCNA.R                  →  module identification + GS/MM + KEGG enrichment
04.DEG_module_intersect.R   →  module ∩ DEG intersection
05.intersection_network.R   →  core gene-centered network plot
06.metab_qc.R               →  metabolomics QC + normalization (V2 data)
07.gene_metab_cor.R         →  gene–metabolite correlation + enrichment + figures
utils.R                     →  shared functions (save_plot, WGCNA settings)

Dependencies

• R ≥ 4.3
• WGCNA (≥ 1.72)
• DESeq2 (≥ 1.42)
• clusterProfiler (≥ 4.10)
• ggraph (≥ 2.2), tidygraph, igraph, ggrepel, ggnewscale
• MetaboAnalystR (optional; pathway analysis can be run via the web server)
• dplyr, tidyr, readr, readxl, ggplot2, scales

Run order

1. 01.data_preprocessing.R — Extract gene-level counts (StringTie) and FPKM (ballgown) from raw data; only g-prefix reference genes are retained (12,094 genes)
2. 02.DESeq2_analysis.R — Differential expression (7 comparisons vs M_0d, padj<0.001, |log2FC|≥2) → outputs 03_DEG_union_up.tsv / 03_DEG_union_down.tsv
3. 03.WGCNA.R — Signed co-expression network (bicor, β=12 auto-selected, R²=0.80), blockwiseModules (minModuleSize=30, deepSplit=2, mergeCutHeight=0.25); GS/MM analysis used to verify module–trait association (GS>0.6, kME>0.8); module KEGG/GO enrichment
4. 04.DEG_module_intersect.R — Intersection of yellow+pink module genes with up-regulated DEGs (369 candidate genes)
5. 05.intersection_network.R — TOM sub-network → two-step filtering → circular network plot
6. 06.metab_qc.R — Metabolomics QC (vendor V2 peak table, level-1 annotation filter, ≥50% detection, KNN imputation) + PQN normalization + log2
7. 07.gene_metab_cor.R — Group-mean Pearson correlation (369 genes × 837 metabolites, |r|>0.9, FDR<0.05), MetaboAnalyst pathway enrichment, and figure generation (network, heatmap, butterfly plot, treemap)

Key outputs

Output	Description
02.WGCNA/03.all_gene_module_assignment.tsv	All gene module assignments + kME
02.WGCNA/07.module_trait_relationships.pdf	Module–trait correlation heatmap
02.WGCNA/11.kegg_enrichment.pdf	KEGG enrichment dotplot
03.intersection/02.intersection_yellow_pink.tsv	yellow+pink module ∩ DEG (369 genes)
05.intersection_subnetwork/intersection_network_circular.pdf	Core gene-centered circular network
06.metab_association/04.gene_metab_cor_long.tsv	Significant gene–metabolite correlations
06.metab_association/16.pathway_enrichment_heatmap.pdf	Pathway enrichment + hit-metabolite heatmap
06.metab_association/18.gene_correlation_quadrant.pdf	Per-gene positive/negative association butterfly plot

Data availability

Raw RNA-seq reads and LC-MS/MS metabolomics data will be deposited in public repositories (accession numbers to be provided upon submission). Processed input matrices required to reproduce the analyses are included under 00.rawdata/, with a full per-file and per-column description in 00.rawdata/DATA_README.md. A template Data Availability statement for submission is provided in DATA_AVAILABILITY.md.

Note

• Only g-prefix reference genes are used (novel N. transcripts excluded)
• The 369 candidate genes are the full yellow+pink module gene sets intersected with up-regulated DEGs. GS/MM (GS>0.6, kME>0.8) was used to verify module–trait association within step 03, not to filter the intersection set
• Metabolomics uses the vendor-supplied V2 peak table (sample-data.xlsx); level-1 high-confidence annotations only
• Gene–metabolite correlations are computed on group-mean trajectories (10 time points) to match the destructive-sampling design
• utils.R must be sourced by 03.R, 05.R and 07.R (handled automatically)