Monthly Archives: June 2026

Generating report using R (Data_Tam_Metagenomics_2026_Wastewater)

To analyze your metagenomics data with greater control, you can use the R packages phyloseq, ggplot2, vegan, and MaAsLin2. The biobakery workflow has already generated the merged tables you need in the results/ directory.

Below is the complete R code to:

Load the data and reproduce the figures from wmgx_report.pdf (PCoA, Heatmap, Stacked Barplot).
Perform differential abundance analysis between Pre-treatment and Post-treatment.
Statistically test if the differences between time points are significant.

Prerequisites

Make sure you have the following R packages installed:

install.packages(c("phyloseq", "ggplot2", "vegan", "dplyr", "tidyr", "pheatmap"))
# MaAsLin2 is installed via Bioconductor
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install("MaAsLin2")

Step 1: Load Data and Create Metadata

The biobakery workflow merged your samples into single tables. We will load the MetaPhlAn species table and HUMAnN pathway table, and map your experimental design.

library(phyloseq)
library(ggplot2)
library(dplyr)
library(tidyr)

# 1. Load MetaPhlAn taxonomic profiles
mpa_file <- "results/metaphlan/merged/metaphlan_taxonomic_profiles.tsv"
mpa_raw <- read.table(mpa_file, header=TRUE, sep="\t", comment.char="", check.names=FALSE)
colnames(mpa_raw)[1] <- "clade_name"
rownames(mpa_raw) <- mpa_raw$clade_name
mpa_otu <- mpa_raw[, -1]

# Filter for species level (s__) and clean names
species_idx <- grep("s__", rownames(mpa_otu))
species_otu <- mpa_otu[species_idx, ]
species_names <- gsub(".*s__", "", rownames(species_otu))
species_names <- gsub("\\|.*", "", species_names)
species_names <- gsub("_", " ", species_names) # Replace underscores with spaces
rownames(species_otu) <- species_names

# 2. Create Metadata based on your file names
sample_names <- colnames(species_otu)
metadata <- data.frame(
  row.names = sample_names,
  TimePoint = case_when(
    grepl("Nov", sample_names) ~ "Nov",
    grepl("Jan", sample_names) ~ "Jan",
    grepl("Mar", sample_names) ~ "Mar",
    grepl("May", sample_names) ~ "May"
  ),
  Treatment = ifelse(grepl("Pre", sample_names), "Pre", "Post")
)

# 3. Create phyloseq object for Species
tax_mat <- matrix(rownames(species_otu), ncol=1, dimnames=list(rownames(species_otu), "Species"))
species_ps <- phyloseq(
  otu_table(species_otu, taxa_are_rows=TRUE),
  sample_data(metadata),
  tax_table(tax_mat)
)

# (Optional) Load HUMAnN pathways using the same logic
# path_file <- "results/humann/merged/pathabundance_relab.tsv"
# ... (similar loading steps as above) ...

Step 2: Reproduce Figures from `wmgx_report.pdf`

1. Ordination (PCoA)

Reproduces the Bray-Curtis Principal Coordinate Analysis.

library(vegan)

# Calculate Bray-Curtis distance and perform PCoA
ord <- ordinate(species_ps, method="PCoA", distance="bray")

# Plot
plot_ordination(species_ps, ord, color="Treatment", shape="TimePoint") +
  geom_point(size=4) +
  theme_minimal() +
  labs(title="PCoA of Species Profiles (Bray-Curtis)", 
       subtitle="Colored by Treatment, Shaped by Time Point") +
  theme(legend.position = "right")

2. Heatmap (Top 25 Species)

Reproduces the hierarchical clustering heatmap using Euclidean distance.

# Get top 25 most abundant species
top25 <- names(sort(taxa_sums(species_ps), decreasing=TRUE))[1:25]
ps_top25 <- prune_taxa(top25, species_ps)

# Plot heatmap (method=NULL defaults to hierarchical clustering)
plot_heatmap(ps_top25, method=NULL, distance="euclidean", 
             taxa.labels=TRUE, label.size=3, na.value=NA) +
  scale_fill_viridis_c() +
  labs(title="Heatmap of Top 25 Most Abundant Species")

3. Stacked Barplot (Top 15 Species)

Reproduces the stacked barplot, sorted by total abundance.

# Get top 15 species
top15 <- names(sort(taxa_sums(species_ps), decreasing=TRUE))[1:15]
ps_top15 <- prune_taxa(top15, species_ps)

# Melt data for ggplot2
df <- psmelt(ps_top15)

# Sort samples by the total abundance of these top 15 species (decreasing order)
sample_order <- df %>%
  group_by(Sample) %>%
  summarise(Total_Ab = sum(Abundance)) %>%
  arrange(desc(Total_Ab)) %>%
  pull(Sample)

df$Sample <- factor(df$Sample, levels=sample_order)

# Plot
ggplot(df, aes(x=Sample, y=Abundance, fill=Species)) +
  geom_bar(stat="identity", position="stack") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle=90, hjust=1, vjust=0.5)) +
  labs(title="Relative Abundance of Top 15 Most Abundant Species", 
       y="Relative Abundance") +
  scale_fill_discrete(name = "Species")

Step 3: Differential Abundance Analysis (Pre vs. Post)

To find which specific species or pathways are significantly different between Pre- and Post-treatment, we use MaAsLin2 (the official biobakery recommendation). It allows us to test the effect of Treatment while adjusting for TimePoint as a covariate.

library(MaAsLin2)

# MaAsLin2 requires samples as rows and features as columns
species_t <- t(otu_table(species_ps))
species_df <- as.data.frame(species_t)
species_df$sample_id <- rownames(species_df)

metadata_df <- data.frame(
  sample_id = rownames(metadata),
  metadata
)

# Run MaAsLin2
# We set "Pre" and "Nov" as reference groups. 
# The output coefficients will represent the difference of Post vs Pre.
fit_results <- MaAsLin2(
  input_data = species_df,
  input_metadata = metadata_df,
  output = "maaslin2_species_output",
  fixed_effects = c("Treatment", "TimePoint"),
  reference = c("Treatment,Pre", "TimePoint,Nov"), 
  random_effects = c(),
  correction = "FDR",
  standardize = FALSE,
  min_abundance = 0.001, 
  min_prevalence = 0.1   # Feature must be present in at least 10% of samples
)

# View significant results (q-value < 0.05)
sig_results <- subset(fit_results$results, qval < 0.05 & metadata == "Treatment")
if(nrow(sig_results) > 0) {
  print(sig_results[, c("feature", "coef", "qval")] %>% arrange(qval))
} else {
  print("No significantly different species found at FDR < 0.05.")
}

Note: You can run the exact same code for pathways by replacing species_df with your HUMAnN pathway data frame.

Step 4: Is the difference between time points significant?

To determine if the microbial community composition changes significantly across the different time points (Nov, Jan, Mar, May), we use PERMANOVA (via the adonis2 function in the vegan package). This tests the overall community structure rather than individual features.

library(vegan)

# Calculate Bray-Curtis distance matrix
dist_bc <- vegdist(t(otu_table(species_ps)), method="bray")

# Run PERMANOVA
# This tests the main effects of Treatment and TimePoint, and their interaction
adonis_res <- adonis2(dist_bc ~ Treatment * TimePoint, data=metadata, permutations=999)
print(adonis_res)

How to interpret the `adonis2` output:

Treatment row: Look at the Pr(>F) (p-value). If it is < 0.05, it means there is a statistically significant difference in the overall microbiome between Pre-treatment and Post-treatment.
TimePoint row: Look at the Pr(>F). If it is < 0.05, it means the microbial community composition differs significantly across the months (Nov, Jan, Mar, May).
Treatment:TimePoint row: If this interaction is significant (< 0.05), it means the effect of the treatment depends on the time of year (e.g., the treatment works differently in November compared to May).

Important Note on Statistical Power: Your dataset has 11 samples total. In the PERMANOVA model, this leaves only 3 residual degrees of freedom. While this mathematically allows for testing, the statistical power is low. Significant results are highly trustworthy, but non-significant results might simply be due to the small sample size.

Processing RNAseq (Data_Tam_RNAseq_2026_Dicl_Mero_Azith_Rifa_on_AYE)

Leave a reply

Preparing raw data

 mkdir raw_data; cd raw_data

 # control samples (8)
 ln -s ../X101SC26025981-Z01-J001/01.RawData/1/1_1.fq.gz AYE-WT_ctr_r1_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/1/1_2.fq.gz AYE-WT_ctr_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/2/2_1.fq.gz AYE-WT_ctr_r2_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/2/2_2.fq.gz AYE-WT_ctr_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/3/3_1.fq.gz AYE-WT_ctr_r3_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/3/3_2.fq.gz AYE-WT_ctr_r3_R2.fastq.gz

 ln -s ../X101SC26025981-Z01-J001/01.RawData/4/4_1.fq.gz AYE-T_ctr_r1_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/4/4_2.fq.gz AYE-T_ctr_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/5/5_1.fq.gz AYE-T_ctr_r2_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/5/5_2.fq.gz AYE-T_ctr_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/6/6_1.fq.gz AYE-T_ctr_r3_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/6/6_2.fq.gz AYE-T_ctr_r3_R2.fastq.gz

 ln -s ../X101SC26025981-Z01-J001/01.RawData/7/7_1.fq.gz AYE-O_ctr_r1_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/7/7_2.fq.gz AYE-O_ctr_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/8/8_1.fq.gz AYE-O_ctr_r2_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/8/8_2.fq.gz AYE-O_ctr_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/9/9_1.fq.gz AYE-O_ctr_r3_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/9/9_2.fq.gz AYE-O_ctr_r3_R2.fastq.gz

 ln -s ../X101SC26025981-Z01-J001/01.RawData/10/10_1.fq.gz O-Trans_ctr_r1_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/10/10_2.fq.gz O-Trans_ctr_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/11/11_1.fq.gz O-Trans_ctr_r2_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/11/11_2.fq.gz O-Trans_ctr_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/12/12_1.fq.gz O-Trans_ctr_r3_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/12/12_2.fq.gz O-Trans_ctr_r3_R2.fastq.gz

 ln -s ../X101SC26025981-Z01-J001/01.RawData/1new/1new_1.fq.gz WT-Trans_ctr_r1_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/1new/1new_2.fq.gz WT-Trans_ctr_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/2new/2new_1.fq.gz WT-Trans_ctr_r2_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/2new/2new_2.fq.gz WT-Trans_ctr_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/3new/3new_1.fq.gz WT-Trans_ctr_r3_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/3new/3new_2.fq.gz WT-Trans_ctr_r3_R2.fastq.gz

 ln -s ../X101SC26025981-Z01-J001/01.RawData/49/49_1.fq.gz AYE-WT_ctr_solid_r1_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/49/49_2.fq.gz AYE-WT_ctr_solid_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/50/50_1.fq.gz AYE-WT_ctr_solid_r2_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/50/50_2.fq.gz AYE-WT_ctr_solid_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/51/51_1.fq.gz AYE-WT_ctr_solid_r3_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/51/51_2.fq.gz AYE-WT_ctr_solid_r3_R2.fastq.gz

 ln -s ../X101SC26025981-Z01-J001/01.RawData/52/52_1.fq.gz AYE-O_ctr_solid_r1_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/52/52_2.fq.gz AYE-O_ctr_solid_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/53/53_1.fq.gz AYE-O_ctr_solid_r2_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/53/53_2.fq.gz AYE-O_ctr_solid_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/54/54_1.fq.gz AYE-O_ctr_solid_r3_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/54/54_2.fq.gz AYE-O_ctr_solid_r3_R2.fastq.gz

 ln -s ../X101SC26025981-Z01-J001/01.RawData/55/55_1.fq.gz AYE-T_ctr_solid_r1_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/55/55_2.fq.gz AYE-T_ctr_solid_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/56/56_1.fq.gz AYE-T_ctr_solid_r2_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/56/56_2.fq.gz AYE-T_ctr_solid_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/57/57_1.fq.gz AYE-T_ctr_solid_r3_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/57/57_2.fq.gz AYE-T_ctr_solid_r3_R2.fastq.gz

 # Diclofenac（双氯芬酸）treatment (6)
 ln -s ../X101SC26025981-Z01-J001/01.RawData/25/25_1.fq.gz AYE-WT_Diclo750_r1_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/25/25_2.fq.gz AYE-WT_Diclo750_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/26/26_1.fq.gz AYE-WT_Diclo750_r2_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/26/26_2.fq.gz AYE-WT_Diclo750_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/27/27_1.fq.gz AYE-WT_Diclo750_r3_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/27/27_2.fq.gz AYE-WT_Diclo750_r3_R2.fastq.gz

 ln -s ../X101SC26025981-Z01-J001/01.RawData/28/28_1.fq.gz AYE-T_Diclo375_r1_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/28/28_2.fq.gz AYE-T_Diclo375_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/29/29_1.fq.gz AYE-T_Diclo375_r2_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/29/29_2.fq.gz AYE-T_Diclo375_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/30/30_1.fq.gz AYE-T_Diclo375_r3_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/30/30_2.fq.gz AYE-T_Diclo375_r3_R2.fastq.gz

 ln -s ../X101SC26025981-Z01-J001/01.RawData/31/31_1.fq.gz AYE-O_Diclo375_r1_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/31/31_2.fq.gz AYE-O_Diclo375_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/32/32_1.fq.gz AYE-O_Diclo375_r2_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/32/32_2.fq.gz AYE-O_Diclo375_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/33/33_1.fq.gz AYE-O_Diclo375_r3_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/33/33_2.fq.gz AYE-O_Diclo375_r3_R2.fastq.gz

 ln -s ../X101SC26025981-Z01-J001/01.RawData/34/34_1.fq.gz O-Trans_Diclo375_r1_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/34/34_2.fq.gz O-Trans_Diclo375_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/35/35_1.fq.gz O-Trans_Diclo375_r2_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/35/35_2.fq.gz O-Trans_Diclo375_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/36/36_1.fq.gz O-Trans_Diclo375_r3_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/36/36_2.fq.gz O-Trans_Diclo375_r3_R2.fastq.gz

 ln -s ../X101SC26025981-Z01-J001/01.RawData/4new/4new_1.fq.gz WT-Trans_Diclo750_r1_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/4new/4new_2.fq.gz WT-Trans_Diclo750_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/5new/5new_1.fq.gz WT-Trans_Diclo750_r2_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/5new/5new_2.fq.gz WT-Trans_Diclo750_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/6new/6new_1.fq.gz WT-Trans_Diclo750_r3_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/6new/6new_2.fq.gz WT-Trans_Diclo750_r3_R2.fastq.gz

 ln -s ../X101SC26025981-Z01-J001/01.RawData/73/73_1.fq.gz AYE-WT_Diclo1250_solid_r1_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/73/73_2.fq.gz AYE-WT_Diclo1250_solid_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/74/74_1.fq.gz AYE-WT_Diclo1250_solid_r2_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/74/74_2.fq.gz AYE-WT_Diclo1250_solid_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/75/75_1.fq.gz AYE-WT_Diclo1250_solid_r3_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/75/75_2.fq.gz AYE-WT_Diclo1250_solid_r3_R2.fastq.gz

 # Rifampicin（利福平）treatment (4)
 ln -s ../X101SC26025981-Z01-J001/01.RawData/13/13_1.fq.gz AYE-WT_Rifampicin1.5_r1_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/13/13_2.fq.gz AYE-WT_Rifampicin1.5_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/14/14_1.fq.gz AYE-WT_Rifampicin1.5_r2_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/14/14_2.fq.gz AYE-WT_Rifampicin1.5_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/15/15_1.fq.gz AYE-WT_Rifampicin1.5_r3_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/15/15_2.fq.gz AYE-WT_Rifampicin1.5_r3_R2.fastq.gz

 ln -s ../X101SC26025981-Z01-J001/01.RawData/16/16_1.fq.gz AYE-T_Rifampicin2_r1_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/16/16_2.fq.gz AYE-T_Rifampicin2_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/17/17_1.fq.gz AYE-T_Rifampicin2_r2_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/17/17_2.fq.gz AYE-T_Rifampicin2_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/18/18_1.fq.gz AYE-T_Rifampicin2_r3_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/18/18_2.fq.gz AYE-T_Rifampicin2_r3_R2.fastq.gz

 ln -s ../X101SC26025981-Z01-J001/01.RawData/19/19_1.fq.gz AYE-O_Rifampicin2_r1_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/19/19_2.fq.gz AYE-O_Rifampicin2_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/20/20_1.fq.gz AYE-O_Rifampicin2_r2_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/20/20_2.fq.gz AYE-O_Rifampicin2_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/21/21_1.fq.gz AYE-O_Rifampicin2_r3_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/21/21_2.fq.gz AYE-O_Rifampicin2_r3_R2.fastq.gz

 ln -s ../X101SC26025981-Z01-J001/01.RawData/22/22_1.fq.gz O-Trans_Rifampicin2_r1_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/22/22_2.fq.gz O-Trans_Rifampicin2_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/23/23_1.fq.gz O-Trans_Rifampicin2_r2_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/23/23_2.fq.gz O-Trans_Rifampicin2_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/24/24_1.fq.gz O-Trans_Rifampicin2_r3_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/24/24_2.fq.gz O-Trans_Rifampicin2_r3_R2.fastq.gz

 # Meropenem（美罗培南）treatment (4)
 ln -s ../X101SC26025981-Z01-J001/01.RawData/37/37_1.fq.gz AYE-WT_Mero0.35-0.5_r1_R1.fastq.gz  #AYE-WT_Mero0.5_r1
 ln -s ../X101SC26025981-Z01-J001/01.RawData/37/37_2.fq.gz AYE-WT_Mero0.35-0.5_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/38/38_1.fq.gz AYE-WT_Mero0.35-0.5_r2_R1.fastq.gz  #AYE-WT_YX_Mero0.35_r2
 ln -s ../X101SC26025981-Z01-J001/01.RawData/38/38_2.fq.gz AYE-WT_Mero0.35-0.5_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/39/39_1.fq.gz AYE-WT_Mero0.35-0.5_r3_R1.fastq.gz  #AYE-WT_public_Mero0.35_r3
 ln -s ../X101SC26025981-Z01-J001/01.RawData/39/39_2.fq.gz AYE-WT_Mero0.35-0.5_r3_R2.fastq.gz

 ln -s ../X101SC26025981-Z01-J001/01.RawData/40/40_1.fq.gz AYE-T_Mero0.15_r1_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/40/40_2.fq.gz AYE-T_Mero0.15_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/41/41_1.fq.gz AYE-T_Mero0.15_r2_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/41/41_2.fq.gz AYE-T_Mero0.15_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/42/42_1.fq.gz AYE-T_Mero0.15_r3_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/42/42_2.fq.gz AYE-T_Mero0.15_r3_R2.fastq.gz

 ln -s ../X101SC26025981-Z01-J001/01.RawData/43/43_1.fq.gz AYE-O_Mero0.5_r1_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/43/43_2.fq.gz AYE-O_Mero0.5_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/44/44_1.fq.gz AYE-O_Mero0.5_r2_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/44/44_2.fq.gz AYE-O_Mero0.5_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/45/45_1.fq.gz AYE-O_Mero0.5_r3_R1.fastq.gz  #Mero0.45
 ln -s ../X101SC26025981-Z01-J001/01.RawData/45/45_2.fq.gz AYE-O_Mero0.5_r3_R2.fastq.gz

 ln -s ../X101SC26025981-Z01-J001/01.RawData/46/46_1.fq.gz O-Trans_Mero0.25_r1_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/46/46_2.fq.gz O-Trans_Mero0.25_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/47/47_1.fq.gz O-Trans_Mero0.25_r2_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/47/47_2.fq.gz O-Trans_Mero0.25_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/48/48_1.fq.gz O-Trans_Mero0.25_r3_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/48/48_2.fq.gz O-Trans_Mero0.25_r3_R2.fastq.gz

 # Azithromycin（阿奇霉素）treatment (5), among them, F_ctr_solid is clinical isolate.
 ln -s ../X101SC26025981-Z01-J001/01.RawData/58/58_1.fq.gz F_ctr_solid_r1_R1.fastq.gz  #clinical
 ln -s ../X101SC26025981-Z01-J001/01.RawData/58/58_2.fq.gz F_ctr_solid_r1_R2.fastq.gz  #clinical
 ln -s ../X101SC26025981-Z01-J001/01.RawData/59/59_1.fq.gz F_ctr_solid_r2_R1.fastq.gz  #clinical
 ln -s ../X101SC26025981-Z01-J001/01.RawData/59/59_2.fq.gz F_ctr_solid_r2_R2.fastq.gz  #clinical
 ln -s ../X101SC26025981-Z01-J001/01.RawData/60/60_1.fq.gz F_ctr_solid_r3_R1.fastq.gz  #clinical
 ln -s ../X101SC26025981-Z01-J001/01.RawData/60/60_2.fq.gz F_ctr_solid_r3_R2.fastq.gz  #clinical

 ln -s ../X101SC26025981-Z01-J001/01.RawData/61/61_1.fq.gz AYE-WT_Azi20_solid_r1_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/61/61_2.fq.gz AYE-WT_Azi20_solid_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/62/62_1.fq.gz AYE-WT_Azi20_solid_r2_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/62/62_2.fq.gz AYE-WT_Azi20_solid_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/63/63_1.fq.gz AYE-WT_Azi20_solid_r3_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/63/63_2.fq.gz AYE-WT_Azi20_solid_r3_R2.fastq.gz

 ln -s ../X101SC26025981-Z01-J001/01.RawData/67/67_1.fq.gz AYE-T_Azi20_solid_r1_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/67/67_2.fq.gz AYE-T_Azi20_solid_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/68/68_1.fq.gz AYE-T_Azi20_solid_r2_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/68/68_2.fq.gz AYE-T_Azi20_solid_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/69/69_1.fq.gz AYE-T_Azi20_solid_r3_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/69/69_2.fq.gz AYE-T_Azi20_solid_r3_R2.fastq.gz

 ln -s ../X101SC26025981-Z01-J001/01.RawData/64/64_1.fq.gz AYE-O_Azi20_solid_r1_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/64/64_2.fq.gz AYE-O_Azi20_solid_r1_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/65/65_1.fq.gz AYE-O_Azi20_solid_r2_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/65/65_2.fq.gz AYE-O_Azi20_solid_r2_R2.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/66/66_1.fq.gz AYE-O_Azi20_solid_r3_R1.fastq.gz
 ln -s ../X101SC26025981-Z01-J001/01.RawData/66/66_2.fq.gz AYE-O_Azi20_solid_r3_R2.fastq.gz

 ln -s ../X101SC26025981-Z01-J001/01.RawData/70/70_1.fq.gz F_Azi20_solid_r1_R1.fastq.gz  #clinical
 ln -s ../X101SC26025981-Z01-J001/01.RawData/70/70_2.fq.gz F_Azi20_solid_r1_R2.fastq.gz  #clinical
 ln -s ../X101SC26025981-Z01-J001/01.RawData/71/71_1.fq.gz F_Azi20_solid_r2_R1.fastq.gz  #clinical
 ln -s ../X101SC26025981-Z01-J001/01.RawData/71/71_2.fq.gz F_Azi20_solid_r2_R2.fastq.gz  #clinical
 ln -s ../X101SC26025981-Z01-J001/01.RawData/72/72_1.fq.gz F_Azi20_solid_r3_R1.fastq.gz  #clinical
 ln -s ../X101SC26025981-Z01-J001/01.RawData/72/72_2.fq.gz F_Azi20_solid_r3_R2.fastq.gz  #clinical

Preparing the directory trimmed

 mkdir trimmed trimmed_unpaired;
 for sample_id in AYE-O_Azi20_solid_r1 AYE-O_Azi20_solid_r2 AYE-O_Azi20_solid_r3 AYE-O_ctr_r1 AYE-O_ctr_r2 AYE-O_ctr_r3 AYE-O_ctr_solid_r1 AYE-O_ctr_solid_r2 AYE-O_ctr_solid_r3 AYE-O_Diclo375_r1 AYE-O_Diclo375_r2 AYE-O_Diclo375_r3 AYE-O_Mero0.5_r1 AYE-O_Mero0.5_r2 AYE-O_Mero0.5_r3 AYE-O_Rifampicin2_r1 AYE-O_Rifampicin2_r2 AYE-O_Rifampicin2_r3 AYE-T_Azi20_solid_r1 AYE-T_Azi20_solid_r2 AYE-T_Azi20_solid_r3 AYE-T_ctr_r1 AYE-T_ctr_r2 AYE-T_ctr_r3 AYE-T_ctr_solid_r1 AYE-T_ctr_solid_r2 AYE-T_ctr_solid_r3 AYE-T_Diclo375_r1 AYE-T_Diclo375_r2 AYE-T_Diclo375_r3 AYE-T_Mero0.15_r1 AYE-T_Mero0.15_r2 AYE-T_Mero0.15_r3 AYE-T_Rifampicin2_r1 AYE-T_Rifampicin2_r2 AYE-T_Rifampicin2_r3 AYE-WT_Azi20_solid_r1 AYE-WT_Azi20_solid_r2 AYE-WT_Azi20_solid_r3 AYE-WT_ctr_r1 AYE-WT_ctr_r2 AYE-WT_ctr_r3 AYE-WT_ctr_solid_r1 AYE-WT_ctr_solid_r2 AYE-WT_ctr_solid_r3 AYE-WT_Diclo1250_solid_r1 AYE-WT_Diclo1250_solid_r2 AYE-WT_Diclo1250_solid_r3 AYE-WT_Diclo750_r1 AYE-WT_Diclo750_r2 AYE-WT_Diclo750_r3 AYE-WT_Mero0.35-0.5_r1 AYE-WT_Mero0.35-0.5_r2 AYE-WT_Mero0.35-0.5_r3 AYE-WT_Rifampicin1.5_r1 AYE-WT_Rifampicin1.5_r2 AYE-WT_Rifampicin1.5_r3 F_Azi20_solid_r1 F_Azi20_solid_r2 F_Azi20_solid_r3 F_ctr_solid_r1 F_ctr_solid_r2 F_ctr_solid_r3 O-Trans_ctr_r1 O-Trans_ctr_r2 O-Trans_ctr_r3 O-Trans_Diclo375_r1 O-Trans_Diclo375_r2 O-Trans_Diclo375_r3 O-Trans_Mero0.25_r1 O-Trans_Mero0.25_r2 O-Trans_Mero0.25_r3 O-Trans_Rifampicin2_r1 O-Trans_Rifampicin2_r2 O-Trans_Rifampicin2_r3 WT-Trans_ctr_r1 WT-Trans_ctr_r2 WT-Trans_ctr_r3 WT-Trans_Diclo750_r1 WT-Trans_Diclo750_r2 WT-Trans_Diclo750_r3; do \
 for sample_id in AYE-T_Diclo375_r2; do \
         java -jar /home/jhuang/Tools/Trimmomatic-0.36/trimmomatic-0.36.jar PE -threads 100 raw_data/${sample_id}_R1.fastq.gz raw_data/${sample_id}_R2.fastq.gz trimmed/${sample_id}_R1.fastq.gz trimmed_unpaired/${sample_id}_R1.fastq.gz trimmed/${sample_id}_R2.fastq.gz trimmed_unpaired/${sample_id}_R2.fastq.gz ILLUMINACLIP:/home/jhuang/Tools/Trimmomatic-0.36/adapters/TruSeq3-PE-2.fa:2:30:10:8:TRUE LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 AVGQUAL:20; done 2> trimmomatic_pe.log;
 done

(Optional) using trinity to find the most closely reference

 #Trinity --seqType fq --max_memory 50G --left trimmed/wt_r1_R1.fastq.gz  --right trimmed/wt_r1_R2.fastq.gz --CPU 12

 #https://www.genome.jp/kegg/tables/br08606.html#prok
 acb     KGB     Acinetobacter baumannii ATCC 17978  2007    GenBank
 abm     KGB     Acinetobacter baumannii SDF     2008    GenBank
 aby     KGB     Acinetobacter baumannii AYE     2008    GenBank --> *
 abc     KGB     Acinetobacter baumannii ACICU   2008    GenBank
 abn     KGB     Acinetobacter baumannii AB0057  2008    GenBank
 abb     KGB     Acinetobacter baumannii AB307-0294  2008    GenBank
 abx     KGB     Acinetobacter baumannii 1656-2  2012    GenBank
 abz     KGB     Acinetobacter baumannii MDR-ZJ06    2012    GenBank
 abr     KGB     Acinetobacter baumannii MDR-TJ  2012    GenBank
 abd     KGB     Acinetobacter baumannii TCDC-AB0715     2012    GenBank
 abh     KGB     Acinetobacter baumannii TYTH-1  2012    GenBank
 abad    KGB     Acinetobacter baumannii D1279779    2013    GenBank
 abj     KGB     Acinetobacter baumannii BJAB07104   2013    GenBank
 abab    KGB     Acinetobacter baumannii BJAB0715    2013    GenBank
 abaj    KGB     Acinetobacter baumannii BJAB0868    2013    GenBank
 abaz    KGB     Acinetobacter baumannii ZW85-1  2013    GenBank
 abk     KGB     Acinetobacter baumannii AbH12O-A2   2014    GenBank
 abau    KGB     Acinetobacter baumannii AB030   2014    GenBank
 abaa    KGB     Acinetobacter baumannii AB031   2014    GenBank
 abw     KGB     Acinetobacter baumannii AC29    2014    GenBank
 abal    KGB     Acinetobacter baumannii LAC-4   2015    GenBank
 #Note that the Acinetobacter baumannii strain ATCC 19606 chromosome, complete genome (GenBank: CU459141.1) was choosen as reference!

Preparing samplesheet.csv

 sample,fastq_1,fastq_2,strandedness
 Urine_r1,Urine_r1_R1.fq.gz,Urine_r1_R2.fq.gz,auto
 ...

Downloading CU459141.fasta and CU459141.gff from GenBank and preparing CU459141_m.gff

 #Example1: http://xgenes.com/article/article-content/157/prepare-virus-gtf-for-nextflow-run/
 #Default NOT_WORKING: --gtf_group_features 'gene_id'  --gtf_extra_attributes 'gene_name' --featurecounts_group_type 'gene_biotype' --featurecounts_feature_type 'exon'
 #(host_env) !NOT_WORKING! jhuang@WS-2290C:~/DATA/Data_Tam_RNAseq_2024$ /usr/local/bin/nextflow run rnaseq/main.nf --input samplesheet.csv --outdir results    --fasta "/home/jhuang/DATA/Data_Tam_RNAseq_2024/CU459141.fasta" --gff "/home/jhuang/DATA/Data_Tam_RNAseq_2024/CU459141.gff"        -profile docker -resume  --max_cpus 55 --max_memory 512.GB --max_time 2400.h    --save_align_intermeds --save_unaligned --save_reference    --aligner 'star_salmon'    --gtf_group_features 'gene_id'  --gtf_extra_attributes 'gene_name' --featurecounts_group_type 'gene_biotype' --featurecounts_feature_type 'transcript'

 # -- DEBUG_1 (CDS --> exon in CP059040.gff) --
 #Checking the record (see below) in results/genome/CP059040.gtf
 #In ./results/genome/CP059040.gtf e.g. "CP059040.1      Genbank transcript      1       1398    .       +       .       transcript_id "gene-H0N29_00005"; gene_id "gene-H0N29_00005"; gene_name "dnaA"; Name "dnaA"; gbkey "Gene"; gene "dnaA"; gene_biotype "protein_coding"; locus_tag "H0N29_00005";"
 #--featurecounts_feature_type 'transcript' returns only the tRNA results
 #Since the tRNA records have "transcript and exon". In gene records, we have "transcript and CDS". replace the CDS with exon

 grep -P "\texon\t" CP059040.gff | sort | wc -l    #96
 grep -P "cmsearch\texon\t" CP059040.gff | wc -l    #=10  ignal recognition particle sRNA small typ, transfer-messenger RNA, 5S ribosomal RNA
 grep -P "Genbank\texon\t" CP059040.gff | wc -l    #=12  16S and 23S ribosomal RNA
 grep -P "tRNAscan-SE\texon\t" CP059040.gff | wc -l    #tRNA 74
 wc -l star_salmon/AUM_r3/quant.genes.sf  #--featurecounts_feature_type 'transcript' results in 96 records!

 grep -P "\tCDS\t" CU459141.gff3 | wc -l  #3659
 sed 's/\tCDS\t/\texon\t/g' CU459141.gff3 > CU459141_m.gff
 grep -P "\texon\t" CU459141_m.gff | sort | wc -l  #3760

 # -- DEBUG_2: combination of 'CU459141_m.gff' and 'exon' results in ERROR, using 'transcript' instead!
 --gff "/home/jhuang/DATA/Data_Tam_RNAseq_2026_on_AYE/CU459141_m.gff" --featurecounts_feature_type 'transcript'

 # -- DEBUG_3: make sure the header of fasta is the same to the *_m.gff file

nextflow run

 # ---- SUCCESSFUL with directly downloaded gff3 and fasta from NCBI using docker after replacing 'CDS' with 'exon' ----
 (host_env) mv trimmed/*.fastq.gz .
 (host_env) nextflow run nf-core/rnaseq -r 3.14.0 -profile docker \
 --input samplesheet.csv --outdir results    --fasta "/home/jhuang/DATA/Data_Tam_RNAseq_2026_on_AYE/CU459141.fasta" --gff "/home/jhuang/DATA/Data_Tam_RNAseq_2026_on_AYE/CU459141_m.gff"        -resume  --max_cpus 90 --max_memory 900.GB --max_time 2400.h    --save_align_intermeds --save_unaligned --save_reference    --aligner 'star_salmon'    --gtf_group_features 'gene_id'  --gtf_extra_attributes 'gene_name' --featurecounts_group_type 'gene_biotype' --featurecounts_feature_type 'transcript'

(OPTIONAL, since the R-script complete_deg_pipeline_custom_cutoff.R runs completely) Generate PCA-plot

 #mamba activate r_env

 #install.packages("ggfun")
 # Import the required libraries
 library("AnnotationDbi")
 library("clusterProfiler")
 library("ReactomePA")
 library(gplots)
 library(tximport)
 library(DESeq2)
 #library("org.Hs.eg.db")
 library(dplyr)
 library(tidyverse)
 #install.packages("devtools")
 #devtools::install_version("gtable", version = "0.3.0")
 library(gplots)
 library("RColorBrewer")
 #install.packages("ggrepel")
 library("ggrepel")
 # install.packages("openxlsx")
 library(openxlsx)
 library(EnhancedVolcano)
 library(DESeq2)
 library(edgeR)

 setwd("~/DATA/Data_Tam_RNAseq_2026_Dicl_Mero_Azith_Rifa_on_AYE/results/star_salmon")
 # Define paths to your Salmon output quantification files

 # Store sample names in a character vector
 samples <- c(
     "AYE-O_Azi20_solid_r1", "AYE-O_Azi20_solid_r2", "AYE-O_Azi20_solid_r3", "AYE-O_ctr_r1", "AYE-O_ctr_r2",
     "AYE-O_ctr_r3", "AYE-O_ctr_solid_r1", "AYE-O_ctr_solid_r2", "AYE-O_ctr_solid_r3",
     "AYE-O_Diclo375_r1", "AYE-O_Diclo375_r2", "AYE-O_Diclo375_r3", "AYE-O_Mero0.5_r1",
     "AYE-O_Mero0.5_r2", "AYE-O_Mero0.5_r3", "AYE-O_Rifampicin2_r1", "AYE-O_Rifampicin2_r2",
     "AYE-O_Rifampicin2_r3", "AYE-T_Azi20_solid_r1", "AYE-T_Azi20_solid_r2", "AYE-T_Azi20_solid_r3",
     "AYE-T_ctr_r1", "AYE-T_ctr_r2", "AYE-T_ctr_r3", "AYE-T_ctr_solid_r1", "AYE-T_ctr_solid_r2",
     "AYE-T_ctr_solid_r3", "AYE-T_Diclo375_r1", "AYE-T_Diclo375_r2", "AYE-T_Diclo375_r3",
     "AYE-T_Mero0.15_r1", "AYE-T_Mero0.15_r2", "AYE-T_Mero0.15_r3", "AYE-T_Rifampicin2_r1",
     "AYE-T_Rifampicin2_r2", "AYE-T_Rifampicin2_r3", "AYE-WT_Azi20_solid_r1", "AYE-WT_Azi20_solid_r2",
     "AYE-WT_Azi20_solid_r3", "AYE-WT_ctr_r1", "AYE-WT_ctr_r2", "AYE-WT_ctr_r3", "AYE-WT_ctr_solid_r1",
     "AYE-WT_ctr_solid_r2", "AYE-WT_ctr_solid_r3", "AYE-WT_Diclo1250_solid_r1", "AYE-WT_Diclo1250_solid_r2",
     "AYE-WT_Diclo1250_solid_r3", "AYE-WT_Diclo750_r1", "AYE-WT_Diclo750_r2", "AYE-WT_Diclo750_r3",
     "AYE-WT_Mero0.35-0.5_r1", "AYE-WT_Mero0.35-0.5_r2", "AYE-WT_Mero0.35-0.5_r3", "AYE-WT_Rifampicin1.5_r1",
     "AYE-WT_Rifampicin1.5_r2", "AYE-WT_Rifampicin1.5_r3", "F_Azi20_solid_r1", "F_Azi20_solid_r2",
     "F_Azi20_solid_r3", "F_ctr_solid_r1", "F_ctr_solid_r2", "F_ctr_solid_r3",
     "O-Trans_ctr_r1","O-Trans_ctr_r2","O-Trans_ctr_r3",  "O-Trans_Diclo375_r1","O-Trans_Diclo375_r2","O-Trans_Diclo375_r3",
     "O-Trans_Mero0.25_r1", "O-Trans_Mero0.25_r2", "O-Trans_Mero0.25_r3", "O-Trans_Rifampicin2_r1",
     "O-Trans_Rifampicin2_r2", "O-Trans_Rifampicin2_r3", "WT-Trans_ctr_r1", "WT-Trans_ctr_r2",
     "WT-Trans_ctr_r3", "WT-Trans_Diclo750_r1", "WT-Trans_Diclo750_r2", "WT-Trans_Diclo750_r3"
 )

 # Automatically generate the named vector
 files <- setNames(paste0("./", samples, "/quant.sf"), samples)
 # Import the transcript abundance data with tximport
 txi <- tximport(files, type = "salmon", txIn = TRUE, txOut = TRUE)

 # -----------------------------------------------------------------
 # ---- Step 1: Create Detailed Metadata from Your Sample Names ----

 # Extract metadata from sample names
 samples <- names(files)

 # Parse the complex sample names
 metadata <- data.frame(
 sample = samples,
 stringsAsFactors = FALSE
 )

 # Extract strain (everything before first underscore or hyphen treatment)
 metadata$strain <- sapply(strsplit(samples, "[-_]"), function(x) {
 if(x[1] %in% c("AYE", "O", "WT", "F")) {
     if(x[1] == "AYE" && length(x) > 1 && x[2] %in% c("WT", "T", "O")) {
     paste(x[1:2], collapse = "-")
     } else if(x[1] %in% c("O", "WT") && x[2] == "Trans") {
     paste(x[1:2], collapse = "-")
     } else {
     x[1]
     }
 } else {
     x[1]
 }
 })

 # Extract treatment type
 metadata$treatment <- sapply(samples, function(x) {
     if(grepl("_ctr", x)) return("ctrl")
     if(grepl("Diclo", x)) return("Diclo")
     if(grepl("Mero", x)) return("Mero")
     if(grepl("Azi", x)) return("Azi")
     if(grepl("Rifampicin", x)) return("Rifampicin")
     return("ctrl")
 })

 # Extract concentration
 metadata$concentration <- sapply(samples, function(x) {
     if(grepl("Diclo1250", x)) return("1250")
     if(grepl("Diclo750", x)) return("750")
     if(grepl("Diclo375", x)) return("375")
     if(grepl("Mero0.5", x)) return("0.5")
     if(grepl("Mero0.35", x)) return("0.35")
     if(grepl("Mero0.25", x)) return("0.25")
     if(grepl("Mero0.15", x)) return("0.15")
     if(grepl("Azi20", x)) return("20")
     if(grepl("Rifampicin2", x)) return("2")
     if(grepl("Rifampicin1.5", x)) return("1.5")
     return("0")
 })

 # Extract condition (solid vs liquid)
 metadata$condition <- ifelse(grepl("_solid", samples), "solid", "liquid")

 # Extract replicate
 metadata$replicate <- sapply(strsplit(samples, "_"), function(x) {
 rep_part <- x[length(x)]
 gsub("r", "", rep_part)
 })

 # Create combined group for easy comparisons
 metadata$group <- paste(metadata$strain, metadata$treatment, metadata$concentration, sep = "_")

 # Set row names
 rownames(metadata) <- metadata$sample

 # Reorder to match txi columns
 metadata <- metadata[colnames(txi$counts), ]

 # ---------------------------------------------
 # ---- Step 2: Choose Your Design Strategy ----

 # Strategy A: Full Factorial Design (if balanced)
 #dds <- DESeqDataSetFromTximport(txi, metadata,
 #                         design = ~ strain + treatment + condition)

 # --> Strategy B: Combined Group Factor ⭐ RECOMMENDED
 metadata$group <- factor(paste(metadata$strain,
                                 metadata$treatment,
                                 metadata$concentration,
                                 metadata$condition,
                                 sep = "_"))

 dds <- DESeqDataSetFromTximport(txi, metadata,
                                 design = ~ group)
 dds <- DESeq(dds)

 # See all available comparisons
 resultsNames(dds)

 # -------------------------------------------------------------
 # ---- Step 3: Set Up Specific Comparisons from Your Notes ----
 # ==========================================
 # 1. Define Exact Comparisons from Your Notes
 # ==========================================
 planned_comparisons <- list(
 # --- Baseline / Strain Controls ---
 AYE_T_ctr_vs_AYE_WT_ctr            = list(treat = "AYE-T_ctrl_0_liquid",   ctrl = "AYE-WT_ctrl_0_liquid"),
 AYE_O_ctr_vs_AYE_WT_ctr            = list(treat = "AYE-O_ctrl_0_liquid",   ctrl = "AYE-WT_ctrl_0_liquid"),
 O_Trans_ctr_vs_AYE_WT_ctr          = list(treat = "O-Trans_ctrl_0_liquid", ctrl = "AYE-WT_ctrl_0_liquid"),
 WT_Trans_ctr_vs_AYE_WT_ctr         = list(treat = "WT-Trans_ctrl_0_liquid",ctrl = "AYE-WT_ctrl_0_liquid"),
 AYE_O_ctr_vs_AYE_T                 = list(treat = "AYE-O_ctrl_0_liquid",   ctrl = "AYE-T_ctrl_0_liquid"),
 O_Trans_ctr_vs_AYE_T               = list(treat = "O-Trans_ctrl_0_liquid", ctrl = "AYE-T_ctrl_0_liquid"),
 WT_Trans_ctr_vs_AYE_T              = list(treat = "WT-Trans_ctrl_0_liquid",ctrl = "AYE-T_ctrl_0_liquid"),

 # --- Condition Effects (Solid vs Liquid) ---
 AYE_WT_ctr_solid_vs_AYE_WT_ctr     = list(treat = "AYE-WT_ctrl_0_solid",   ctrl = "AYE-WT_ctrl_0_liquid"),
 AYE_O_ctr_solid_vs_AYE_O_ctr       = list(treat = "AYE-O_ctrl_0_solid",    ctrl = "AYE-O_ctrl_0_liquid"),
 AYE_T_ctr_solid_vs_AYE_T_ctr       = list(treat = "AYE-T_ctrl_0_solid",    ctrl = "AYE-T_ctrl_0_liquid"),
 AYE_O_ctr_solid_vs_AYE_WT_ctr_solid= list(treat = "AYE-O_ctrl_0_solid",    ctrl = "AYE-WT_ctrl_0_solid"),
 AYE_T_ctr_solid_vs_AYE_WT_ctr_solid= list(treat = "AYE-T_ctrl_0_solid",    ctrl = "AYE-WT_ctrl_0_solid"),

 # --- Diclofenac ---
 AYE_WT_Diclo750_vs_AYE_WT_ctr      = list(treat = "AYE-WT_Diclo_750_liquid",   ctrl = "AYE-WT_ctrl_0_liquid"),
 AYE_T_Diclo375_vs_AYE_WT_ctr       = list(treat = "AYE-T_Diclo_375_liquid",    ctrl = "AYE-WT_ctrl_0_liquid"),
 AYE_O_Diclo375_vs_AYE_WT_ctr       = list(treat = "AYE-O_Diclo_375_liquid",    ctrl = "AYE-WT_ctrl_0_liquid"),
 O_Trans_Diclo375_vs_AYE_WT_ctr     = list(treat = "O-Trans_Diclo_375_liquid",  ctrl = "AYE-WT_ctrl_0_liquid"),
 WT_Trans_Diclo750_vs_AYE_WT_ctr    = list(treat = "WT-Trans_Diclo_750_liquid", ctrl = "AYE-WT_ctrl_0_liquid"),
 Diclo_AYE_WT_1250_solid_vs_solid_ctr = list(treat = "AYE-WT_Diclo_1250_solid", ctrl = "AYE-WT_ctrl_0_solid"),

 # --- Meropenem ---
 AYE_WT_Mero_vs_AYE_WT_ctr          = list(treat = "AYE-WT_Mero_0.35_liquid", ctrl = "AYE-WT_ctrl_0_liquid"),
 AYE_T_Mero_vs_AYE_WT_ctr           = list(treat = "AYE-T_Mero_0.15_liquid",      ctrl = "AYE-WT_ctrl_0_liquid"),
 AYE_O_Mero_vs_AYE_WT_ctr           = list(treat = "AYE-O_Mero_0.5_liquid",       ctrl = "AYE-WT_ctrl_0_liquid"),
 O_Trans_Mero_vs_AYE_WT_ctr         = list(treat = "O-Trans_Mero_0.25_liquid",    ctrl = "AYE-WT_ctrl_0_liquid"),
 AYE_T_Mero_vs_AYE_T_ctr            = list(treat = "AYE-T_Mero_0.15_liquid",      ctrl = "AYE-T_ctrl_0_liquid"),

 # --- Azithromycin (Solid) ---
 AYE_WT_Azi_vs_solid_ctr            = list(treat = "AYE-WT_Azi_20_solid", ctrl = "AYE-WT_ctrl_0_solid"),
 AYE_T_Azi_vs_solid_ctr             = list(treat = "AYE-T_Azi_20_solid",  ctrl = "AYE-T_ctrl_0_solid"),
 AYE_O_Azi_vs_solid_ctr             = list(treat = "AYE-O_Azi_20_solid",  ctrl = "AYE-O_ctrl_0_solid"),
 F_Azi_vs_F_solid_ctr               = list(treat = "F_Azi_20_solid",      ctrl = "F_ctrl_0_solid"),

 # --- Rifampicin ---
 AYE_WT_Rif_vs_AYE_WT_ctr           = list(treat = "AYE-WT_Rifampicin_1.5_liquid", ctrl = "AYE-WT_ctrl_0_liquid"),
 AYE_T_Rif_vs_AYE_T_ctr             = list(treat = "AYE-T_Rifampicin_2_liquid",    ctrl = "AYE-T_ctrl_0_liquid"),
 AYE_O_Rif_vs_AYE_O_ctr             = list(treat = "AYE-O_Rifampicin_2_liquid",    ctrl = "AYE-O_ctrl_0_liquid"),
 O_Trans_Rif_vs_O_Trans_ctr         = list(treat = "O-Trans_Rifampicin_2_liquid",  ctrl = "O-Trans_ctrl_0_liquid")
 )

 # Additional Analyses
 planned_comparisons <- list(
 # --- Diclofenac_2 ---
 # * AYE-T_Diclo_375_liquid vs AYE-T_ctrl_0_liquid
 # * AYE-O_Diclo_375_liquid vs AYE-O_ctrl_0_liquid
 # * O-Trans_Diclo_375_liquid vs AYE-O-Trans_ctrl_0_liquid
 # * WT-Trans_Diclo_750_liquid vs AYE-WT-Trans_ctrl_0_liquid
 AYE_T_Diclo375_vs_AYE_T_ctr        = list(treat = "AYE-T_Diclo_375_liquid",        ctrl = "AYE-T_ctrl_0_liquid"),
 AYE_O_Diclo375_vs_AYE_O_ctr        = list(treat = "AYE-O_Diclo_375_liquid",        ctrl = "AYE-O_ctrl_0_liquid"),
 O_Trans_Diclo375_vs_O_Trans_ctr    = list(treat = "O-Trans_Diclo_375_liquid",  ctrl = "O-Trans_ctrl_0_liquid"),
 WT_Trans_Diclo750_vs_WT_Trans_ctr    = list(treat = "WT-Trans_Diclo_750_liquid", ctrl = "WT-Trans_ctrl_0_liquid"),

 # --- Meropenem_2 ---
 # * AYE-T_Mero_0.15_liquid vs AYE-T_ctrl_0_liquid
 # * AYE-O_Mero_0.5_liquid vs AYE-O_ctrl_0_liquid
 # * O-Trans_Mero_0.25_liquid vs AYE-O-Trans_ctrl_0_liquid
 AYE_T_Mero_vs_AYE_T_ctr        = list(treat = "AYE-T_Mero_0.15_liquid",    ctrl = "AYE-T_ctrl_0_liquid"),
 AYE_O_Mero_vs_AYE_O_ctr        = list(treat = "AYE-O_Mero_0.5_liquid",     ctrl = "AYE-O_ctrl_0_liquid"),
 O_Trans_Mero_vs_O_Trans_ctr    = list(treat = "O-Trans_Mero_0.25_liquid",  ctrl = "O-Trans_ctrl_0_liquid")
 )

 # ==========================================
 # 2. Verification & Validation Script
 # ==========================================
 # Identify which column in colData holds your group names
 group_col <- if("group" %in% colnames(colData(dds))) "group" else
             if("treatment" %in% colnames(colData(dds))) "treatment" else
             stop("❌ Please specify the correct colData column containing group names.")

 actual_groups <- unique(colData(dds)[[group_col]])

 cat("\n", paste(rep("=", 85), collapse=""), "\n")
 cat("📋 VERIFICATION OF NOTE-DERIVED COMPARISONS\n")
 cat(paste(rep("=", 85), collapse=""), "\n\n")

 validation_results <- data.frame(
 Comparison_Name = character(),
 Treatment_String = character(),
 Control_String = character(),
 Status = character(),
 Suggested_Contrast = character(),
 stringsAsFactors = FALSE
 )

 for(name in names(planned_comparisons)) {
 trt <- planned_comparisons[[name]]$treat
 ctl <- planned_comparisons[[name]]$ctrl

 # Find closest matches in actual data
 trt_match <- actual_groups[grepl(trt, actual_groups, fixed = TRUE)]
 ctl_match <- actual_groups[grepl(ctl, actual_groups, fixed = TRUE)]

 status <- if(length(trt_match) > 0 && length(ctl_match) > 0) "✅ VALID" else "⚠️  CHECK"
 contrast_str <- if(status == "✅ VALID")
     paste0('c("', group_col, '", "', trt_match[1], '", "', ctl_match[1], '")') else "N/A"

 validation_results <- rbind(validation_results, data.frame(
     Comparison_Name = name,
     Treatment_String = trt,
     Control_String = ctl,
     Status = status,
     Suggested_Contrast = contrast_str,
     stringsAsFactors = FALSE
 ))

 cat(sprintf("%-45s | T:%-25s C:%-20s | %s\n", name, trt, ctl, status))
 if(status == "⚠️  CHECK") {
     if(length(trt_match) == 0) cat("   🔍 Treat not found. Closest: ", paste(head(actual_groups[grepl(strsplit(trt, "_")[[1]][1], actual_groups)], 3), collapse=", "), "\n")
     if(length(ctl_match) == 0) cat("   🔍 Ctrl not found.  Closest: ", paste(head(actual_groups[grepl(strsplit(ctl, "_")[[1]][1], actual_groups)], 3), collapse=", "), "\n")
 }
 }

 # ==========================================
 # 3. Auto-Generate DESeq2 results() Calls (Optional)
 # ==========================================
 valid_comparisons <- validation_results[validation_results$Status == "✅ VALID", ]
 if(nrow(valid_comparisons) > 0) {
 cat("\n📜 READY-TO-RUN DESeq2 CONTRASTS:\n")
 cat(paste(rep("-", 60), collapse=""), "\n")
 for(i in seq_len(nrow(valid_comparisons))) {
     cat(sprintf('res_%s <- results(dds, contrast = %s)\n',
                 gsub("[^A-Za-z0-9]", "_", valid_comparisons$Comparison_Name[i]),
                 valid_comparisons$Suggested_Contrast[i]))
 }
 } else {
 cat("\n⚠️  No exact matches found. Check your colData group naming convention.\n")
 }

 # -----------------------------
 # ---- Step 4: PCA figures ----

 # 🔍 What each figure shows:
 #
 #    01_PCA_by_Strain.png → Tests if genetic background (AYE-WT, AYE-T, AYE-O, Trans, F) is the dominant source of variation.
 #    02_PCA_by_Treatment.png → Shows clustering by antibiotic/drug exposure (ctrl, Diclo, Mero, Azi, Rifampicin).
 #    03_PCA_by_Condition.png → Reveals batch/growth media effects (solid vs liquid).
 #    04_PCA_CombinedGroups.png → Full experimental grouping with labeled sample names for quick outlier detection.
 #    05_PCA_Ellipses.png → Adds 95% confidence boundaries per strain to visualize group spread and overlap.
 #
 # ⚠️ Quick Checklist Before Running:
 #
 #    Ensure metadata columns (strain, treatment, condition, group) are attached to colData(dds).
 #    If ggrepel is missing, run install.packages("ggrepel").
 #    All PNGs will save to your current working directory (getwd()).

 # Install if missing: install.packages(c("ggplot2", "ggrepel"))
 library(DESeq2)
 library(ggplot2)
 library(ggrepel)

 # 1. Variance Stabilizing Transformation & Extract PCA Data
 vsd <- vst(dds, blind = FALSE)
 pca_data <- plotPCA(vsd, intgroup = c("strain", "treatment", "condition", "group"), returnData = TRUE)
 percent_var <- round(100 * attr(pca_data, "percentVar"))

 # Consistent theme for all plots
 base_theme <- theme_bw(base_size = 12) +
 theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 13),
         legend.position = "right",
         legend.title = element_text(face = "bold"),
         panel.grid.major = element_line(color = "grey90"),
         panel.grid.minor = element_blank())

 # --- Plot 1: Colored by Strain ---
 p1 <- ggplot(pca_data, aes(x = PC1, y = PC2, color = strain, shape = condition)) +
 geom_point(size = 3, alpha = 0.8) +
 geom_text_repel(aes(label = name), size = 2.5, max.overlaps = 20, show.legend = FALSE) +
 labs(x = paste0("PC1: ", percent_var[1], "% variance"),
     y = paste0("PC2: ", percent_var[2], "% variance"),
     title = "PCA: Samples Colored by Strain",
     color = "Strain", shape = "Condition") +
 base_theme
 ggsave("01_PCA_by_Strain.png", p1, width = 8, height = 6, dpi = 300)

 # --- Plot 2: Colored by Treatment ---
 p2 <- ggplot(pca_data, aes(x = PC1, y = PC2, color = treatment, shape = condition)) +
 geom_point(size = 3, alpha = 0.8) +
 labs(x = paste0("PC1: ", percent_var[1], "% variance"),
     y = paste0("PC2: ", percent_var[2], "% variance"),
     title = "PCA: Samples Colored by Treatment",
     color = "Treatment", shape = "Condition") +
 base_theme
 ggsave("02_PCA_by_Treatment.png", p2, width = 8, height = 6, dpi = 300)

 # --- Plot 3: Colored by Condition (Solid vs Liquid) ---
 p3 <- ggplot(pca_data, aes(x = PC1, y = PC2, color = condition, shape = strain)) +
 geom_point(size = 3, alpha = 0.8) +
 labs(x = paste0("PC1: ", percent_var[1], "% variance"),
     y = paste0("PC2: ", percent_var[2], "% variance"),
     title = "PCA: Samples Colored by Growth Condition",
     color = "Condition", shape = "Strain") +
 base_theme
 ggsave("03_PCA_by_Condition.png", p3, width = 8, height = 6, dpi = 300)

 # --- Plot 4: Combined Groups with Sample Labels ---
 p4 <- ggplot(pca_data, aes(x = PC1, y = PC2, color = group)) +
 geom_point(size = 3, alpha = 0.8) +
 geom_text_repel(aes(label = name), size = 2, max.overlaps = 30, box.padding = 0.3) +
 labs(x = paste0("PC1: ", percent_var[1], "% variance"),
     y = paste0("PC2: ", percent_var[2], "% variance"),
     title = "PCA: Combined Experimental Groups",
     color = "Group") +
 base_theme + theme(legend.position = "none")
 ggsave("04_PCA_CombinedGroups.png", p4, width = 9, height = 7, dpi = 300)

 # --- Plot 5: 95% Confidence Ellipses (by Strain) ---
 p5 <- ggplot(pca_data, aes(x = PC1, y = PC2, color = strain, fill = strain)) +
 geom_point(size = 3, alpha = 0.7) +
 stat_ellipse(level = 0.95, alpha = 0.2, geom = "polygon", show.legend = FALSE) +
 labs(x = paste0("PC1: ", percent_var[1], "% variance"),
     y = paste0("PC2: ", percent_var[2], "% variance"),
     title = "PCA: 95% Confidence Ellipses by Strain",
     color = "Strain", fill = "Strain") +
 base_theme
 ggsave("05_PCA_Ellipses.png", p5, width = 8, height = 6, dpi = 300)

 message("✅ All 5 PCA plots saved to working directory!")

Run Differential Expression & PCA Analysis Complete

 (r_env) cd ~/DATA/Data_Tam_RNAseq_2026_Dicl_Mero_Azith_Rifa_on_AYE/results/star_salmon/
 #(r_env) Rscript complete_deg_pipeline.R  #For standard cutoff in the project

 # Adapted the script to the following requests:
 # (a) Rifampicin: use genes with a cutoff of log2 fold change > 1.2 and < -1.2 for the KEGG and GO analyses.
 # (b) Baseline / Strain Controls: use genes with a cutoff of log2 fold change > 1.4 and < -1.4 for the KEGG and GO analyses.
 # (c) All other comparisons: please retain the same selection criteria as in the previous analysis you sent to me.

 # How it works:
 #   * Rifampicin: The script looks for "Rif" in the comparison name (e.g., 28_AYE_WT_Rif_vs_Ctrl) and applies |log2FC| >= 1.2.
 #   * Baseline/Strain Controls: The script looks for "_ctr_vs_" in the comparison name (e.g., 01_AYE_T_ctr_vs_AYE_WT_ctr) and applies |log2FC| >= 1.4.
 #   * All Others: Falls back to the original 2.0 cutoff.
 #   * The console output will now explicitly print which cutoff is being used for each specific comparison.

 (r_env) Rscript complete_deg_pipeline_custom_cutoff.R

LOG of Rscript complete_deg_pipeline_custom_cutoff.R

 -rw-rw-r-- 1 jhuang jhuang 355609 Jun 22 12:53 DEG_29_AYE_T_Rif_vs_Ctrl.xlsx
 -rw-rw-r-- 1 jhuang jhuang 253271 Jun 22 12:53 Volcano_29_AYE_T_Rif_vs_Ctrl.png
 -rw-rw-r-- 1 jhuang jhuang 368843 Jun 22 12:53 DEG_30_AYE_O_Rif_vs_Ctrl.xlsx
 -rw-rw-r-- 1 jhuang jhuang 358124 Jun 22 12:53 Volcano_30_AYE_O_Rif_vs_Ctrl.png
 -rw-rw-r-- 1 jhuang jhuang 347126 Jun 22 12:53 DEG_31_O_Trans_Rif_vs_Ctrl.xlsx
 -rw-rw-r-- 1 jhuang jhuang 283473 Jun 22 12:53 Volcano_31_O_Trans_Rif_vs_Ctrl.png

 -rw-rw-r-- 1 jhuang jhuang 352375 Jun 22 12:53 DEG_32_AYE_T_Diclo375_vs_AYE_T_ctr.xlsx
 -rw-rw-r-- 1 jhuang jhuang 257075 Jun 22 12:53 Volcano_32_AYE_T_Diclo375_vs_AYE_T_ctr.png

 -rw-rw-r-- 1 jhuang jhuang  355610 Jun  5 16:26 DEG_29_AYE_T_Rif_vs_Ctrl.xlsx
 -rw-rw-r-- 1 jhuang jhuang  254804 Jun  5 16:26 Volcano_29_AYE_T_Rif_vs_Ctrl.png
 -rw-rw-r-- 1 jhuang jhuang  368843 Jun  5 16:26 DEG_30_AYE_O_Rif_vs_Ctrl.xlsx
 -rw-rw-r-- 1 jhuang jhuang  358417 Jun  5 16:26 Volcano_30_AYE_O_Rif_vs_Ctrl.png
 -rw-rw-r-- 1 jhuang jhuang  347127 Jun  5 16:26 DEG_31_O_Trans_Rif_vs_Ctrl.xlsx
 -rw-rw-r-- 1 jhuang jhuang  284084 Jun  5 16:26 Volcano_31_O_Trans_Rif_vs_Ctrl.png
 -rw-rw-r-- 1 jhuang jhuang    1579 Jun  5 16:26 DEG_Summary_All_31.csv

 (r_env) jhuang@WS-2290C:/mnt/md1/DATA/Data_Tam_RNAseq_2026_Dicl_Mero_Azith_Rifa_on_AYE/results/star_salmon$ Rscript complete_deg_pipeline_custom_cutoff.R
 There were 22 warnings (use warnings() to see them)
 📖 Parsing annotation file to build tx2gene mapping...
 ✅ Created tx2gene mapping with 7520 entries.
 📥 Importing Salmon quantifications...
 reading in files with read_tsv
 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81
 removing duplicated transcript rows from tx2gene
 summarizing abundance
 summarizing counts
 summarizing length
   Note: levels of factors in the design contain characters other than
   letters, numbers, '_' and '.'. It is recommended (but not required) to use
   only letters, numbers, and delimiters '_' or '.', as these are safe characters
   for column names in R. [This is a message, not a warning or an error]
 using counts and average transcript lengths from tximport
 🚀 Running DESeq2 pipeline...
 estimating size factors
   Note: levels of factors in the design contain characters other than
   letters, numbers, '_' and '.'. It is recommended (but not required) to use
   only letters, numbers, and delimiters '_' or '.', as these are safe characters
   for column names in R. [This is a message, not a warning or an error]
 using 'avgTxLength' from assays(dds), correcting for library size
 estimating dispersions
 gene-wise dispersion estimates
 mean-dispersion relationship
   Note: levels of factors in the design contain characters other than
   letters, numbers, '_' and '.'. It is recommended (but not required) to use
   only letters, numbers, and delimiters '_' or '.', as these are safe characters
   for column names in R. [This is a message, not a warning or an error]
 final dispersion estimates
 fitting model and testing
 Available groups in dds:
 [1] "AYE-O_Azi_20_solid"           "AYE-O_ctrl_0_liquid"
 [3] "AYE-O_ctrl_0_solid"           "AYE-O_Diclo_375_liquid"
 [5] "AYE-O_Mero_0.5_liquid"        "AYE-O_Rifampicin_2_liquid"
 [7] "AYE-T_Azi_20_solid"           "AYE-T_ctrl_0_liquid"
 [9] "AYE-T_ctrl_0_solid"           "AYE-T_Diclo_375_liquid"
 [11] "AYE-T_Mero_0.15_liquid"       "AYE-T_Rifampicin_2_liquid"
 [13] "AYE-WT_Azi_20_solid"          "AYE-WT_ctrl_0_liquid"
 [15] "AYE-WT_ctrl_0_solid"          "AYE-WT_Diclo_1250_solid"
 [17] "AYE-WT_Diclo_750_liquid"      "AYE-WT_Mero_0.35-0.5_liquid"
 [19] "AYE-WT_Rifampicin_1.5_liquid" "F_Azi_20_solid"
 [21] "F_ctrl_0_solid"               "O-Trans_ctrl_0_liquid"
 [23] "O-Trans_Diclo_375_liquid"     "O-Trans_Mero_0.25_liquid"
 [25] "O-Trans_Rifampicin_2_liquid"  "WT-Trans_ctrl_0_liquid"
 [27] "WT-Trans_Diclo_750_liquid"
 📖 Parsing annotation file for gene names...

 🚀 PROCESSING  38  COMPARISONS
 Base Thresholds: padj <= 0.05, |log2FC| >= 2 (Dynamic adjustments apply per comparison)
 ================================================================================

 [01/38] 01_AYE_T_ctr_vs_AYE_WT_ctr
   -> Using dynamic LFC cutoff: |log2FC| >= 1.4

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 16, Down: 151)

 [02/38] 02_AYE_O_ctr_vs_AYE_WT_ctr
   -> Using dynamic LFC cutoff: |log2FC| >= 1.4

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 9, Down: 9)

 [03/38] 03_O_Trans_ctr_vs_AYE_WT_ctr
   -> Using dynamic LFC cutoff: |log2FC| >= 1.4

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 10, Down: 44)

 [04/38] 04_WT_Trans_ctr_vs_AYE_WT_ctr
   -> Using dynamic LFC cutoff: |log2FC| >= 1.4

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 15, Down: 153)

 [05/38] 05_AYE_O_ctr_vs_AYE_T_ctr
   -> Using dynamic LFC cutoff: |log2FC| >= 1.4

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 28, Down: 3)

 [06/38] 06_O_Trans_ctr_vs_AYE_T_ctr
   -> Using dynamic LFC cutoff: |log2FC| >= 1.4

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 11, Down: 10)

 [07/38] 07_WT_Trans_ctr_vs_AYE_T_ctr
   -> Using dynamic LFC cutoff: |log2FC| >= 1.4

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 4, Down: 14)

 [08/38] 08_AYE_WT_ctr_solid_vs_liquid
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 485, Down: 324)

 [09/38] 09_AYE_O_ctr_solid_vs_liquid
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 475, Down: 309)

 [10/38] 10_AYE_T_ctr_solid_vs_liquid
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 501, Down: 332)

 [11/38] 11_AYE_O_ctr_solid_vs_AYE_WT_solid
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 12, Down: 5)

 [12/38] 12_AYE_T_ctr_solid_vs_AYE_WT_solid
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 5, Down: 8)

 [13/38] 13_AYE_WT_Diclo750_vs_Ctrl
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 12, Down: 45)

 [14/38] 14_AYE_T_Diclo375_vs_Ctrl
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 75, Down: 181)

 [15/38] 15_AYE_O_Diclo375_vs_Ctrl
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 26, Down: 240)

 [16/38] 16_O_Trans_Diclo375_vs_Ctrl
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 17, Down: 82)

 [17/38] 17_WT_Trans_Diclo750_vs_Ctrl
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 25, Down: 230)

 [18/38] 18_AYE_WT_Diclo1250_solid_vs_Ctrl_solid
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 66, Down: 26)

 [19/38] 19_AYE_WT_Mero_vs_Ctrl
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 359, Down: 175)

 [20/38] 20_AYE_T_Mero_vs_Ctrl
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 93, Down: 110)

 [21/38] 21_AYE_O_Mero_vs_Ctrl
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 267, Down: 244)

 [22/38] 22_O_Trans_Mero_vs_Ctrl
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 183, Down: 153)

 [23/38] 23_AYE_T_Mero_vs_AYE_T_Ctrl
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 91, Down: 27)

 [24/38] 24_AYE_WT_Azi_solid_vs_Ctrl_solid
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 341, Down: 329)

 [25/38] 25_AYE_T_Azi_solid_vs_Ctrl_solid
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 510, Down: 433)

 [26/38] 26_AYE_O_Azi_solid_vs_Ctrl_solid
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 484, Down: 443)

 [27/38] 27_F_Azi_solid_vs_Ctrl_solid
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 147, Down: 271)

 [28/38] 28_AYE_WT_Rif_vs_Ctrl
   -> Using dynamic LFC cutoff: |log2FC| >= 1.2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 8, Down: 11)

 [29/38] 29_AYE_T_Rif_vs_Ctrl
   -> Using dynamic LFC cutoff: |log2FC| >= 1.2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 34, Down: 41)

 [30/38] 30_AYE_O_Rif_vs_Ctrl
   -> Using dynamic LFC cutoff: |log2FC| >= 1.2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 71, Down: 100)

 [31/38] 31_O_Trans_Rif_vs_Ctrl
   -> Using dynamic LFC cutoff: |log2FC| >= 1.2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 6, Down: 10)

 [32/38] 32_AYE_T_Diclo375_vs_AYE_T_ctr
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 47, Down: 6)

 [33/38] 33_AYE_O_Diclo375_vs_AYE_O_ctr
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 8, Down: 73)

 [34/38] 34_O_Trans_Diclo375_vs_O_Trans_ctr
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 2, Down: 17)

 [35/38] 35_WT_Trans_Diclo750_vs_WT_Trans_ctr
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 10, Down: 27)

 [36/38] 36_AYE_T_Mero_vs_AYE_T_ctr
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 91, Down: 27)

 [37/38] 37_AYE_O_Mero_vs_AYE_O_ctr
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 242, Down: 135)

 [38/38] 38_O_Trans_Mero_vs_O_Trans_ctr
   -> Using dynamic LFC cutoff: |log2FC| >= 2

 🔍 DEBUG: Columns in res_df  : baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean
 🔍 DEBUG: Columns in anno_df : gene_id, gene_id_clean, gene_name
 🔍 DEBUG: Columns after join:  baseMean, log2FoldChange, lfcSE, stat, pvalue, padj, GeneID, GeneID_clean, gene_id, gene_name
   ✅ Annotation step completed successfully.
   ✅ Excel + Volcano saved (Up: 177, Down: 50)

 ================================================================================
 📊 FINAL SUMMARY OF ALL 38 COMPARISONS
                                     name total  up down sig_total pct_sig
           # |log2FC| >= 1.4
               01_AYE_T_ctr_vs_AYE_WT_ctr  3609  16  151       167     4.6
               02_AYE_O_ctr_vs_AYE_WT_ctr  3609   9    9        18     0.5
              03_O_Trans_ctr_vs_AYE_WT_ctr  3609  10   44        54     1.5
           04_WT_Trans_ctr_vs_AYE_WT_ctr  3609  15  153       168     4.7
               05_AYE_O_ctr_vs_AYE_T_ctr  3609  28    3        31     0.9
             06_O_Trans_ctr_vs_AYE_T_ctr  3609  11   10        21     0.6
             07_WT_Trans_ctr_vs_AYE_T_ctr  3609   4   14        18     0.5

           # |log2FC| >= 2
           08_AYE_WT_ctr_solid_vs_liquid  3609 485  324       809    22.4
             09_AYE_O_ctr_solid_vs_liquid  3609 475  309       784    21.7
             10_AYE_T_ctr_solid_vs_liquid  3609 501  332       833    23.1
       11_AYE_O_ctr_solid_vs_AYE_WT_solid  3609  12    5        17     0.5
       12_AYE_T_ctr_solid_vs_AYE_WT_solid  3609   5    8        13     0.4
               13_AYE_WT_Diclo750_vs_Ctrl  3609  12   45        57     1.6
               14_AYE_T_Diclo375_vs_Ctrl  3609  75  181       256     7.1
               15_AYE_O_Diclo375_vs_Ctrl  3609  26  240       266     7.4
             16_O_Trans_Diclo375_vs_Ctrl  3609  17   82        99     2.7
             17_WT_Trans_Diclo750_vs_Ctrl  3609  25  230       255     7.1
 18_AYE_WT_Diclo1250_solid_vs_Ctrl_solid  3609  66   26        92     2.5
                   19_AYE_WT_Mero_vs_Ctrl  3609 359  175       534    14.8
                   20_AYE_T_Mero_vs_Ctrl  3609  93  110       203     5.6
                   21_AYE_O_Mero_vs_Ctrl  3609 267  244       511    14.2
                 22_O_Trans_Mero_vs_Ctrl  3609 183  153       336     9.3
             23_AYE_T_Mero_vs_AYE_T_Ctrl  3609  91   27       118     3.3
       24_AYE_WT_Azi_solid_vs_Ctrl_solid  3609 341  329       670    18.6
         25_AYE_T_Azi_solid_vs_Ctrl_solid  3609 510  433       943    26.1
         26_AYE_O_Azi_solid_vs_Ctrl_solid  3609 484  443       927    25.7
             27_F_Azi_solid_vs_Ctrl_solid  3609 147  271       418    11.6

           # |log2FC| >= 1.2
                   28_AYE_WT_Rif_vs_Ctrl  3609   8   11        19     0.5
                     29_AYE_T_Rif_vs_Ctrl  3609  34   41        75     2.1
                     30_AYE_O_Rif_vs_Ctrl  3609  71  100       171     4.7
                   31_O_Trans_Rif_vs_Ctrl  3609   6   10        16     0.4

           # |log2FC| >= 2
                                     name total  up down sig_total
           32_AYE_T_Diclo375_vs_AYE_T_ctr  3609  46    6        52
           33_AYE_O_Diclo375_vs_AYE_O_ctr  3609   8   73        81     2.2
       34_O_Trans_Diclo375_vs_O_Trans_ctr  3609   2   17        19     0.5
     35_WT_Trans_Diclo750_vs_WT_Trans_ctr  3609  10   27        37     1.0
               36_AYE_T_Mero_vs_AYE_T_ctr  3609  91   27       118     3.3
               37_AYE_O_Mero_vs_AYE_O_ctr  3609 242  135       377    10.4
           38_O_Trans_Mero_vs_O_Trans_ctr  3609 177   50       227     6.3

 ✨ All files saved to: /mnt/md1/DATA/Data_Tam_RNAseq_2026_Dicl_Mero_Azith_Rifa_on_AYE/results/star_salmon/DEG_Results_Complete
 📁 Each comparison contains:
     1️ DEG_

.xlsx (3 sheets: Complete_Results, Up_Regulated, Down_Regulated) 2️ Volcano_ .png (GeneName labels for top significant genes) 📤 EXPORTING COMPLETE RESULTS TO CSV FOR KEGG/GO… ✅ Successfully exported 38 CSV files to DEG_Results_Complete

KEGG and GO annotations in non-model organisms

(a) Rifampicin: use genes with a cutoff of log2 fold change > 1.2 and 1.4 and < -1.4 for the KEGG and GO analyses. (c) All other comparisons: please retain the same selection criteria as in the previous analysis you sent to me.

https://www.biobam.com/functional-analysis/

10.1. Assign KEGG and GO Terms (see diagram above)

Since your organism is non-model, standard R databases (org.Hs.eg.db, etc.) won’t work. You’ll need to manually retrieve KEGG and GO annotations.

* Preparing file 1 eggnog_out.emapper.annotations.txt for the R-code below: (KEGG Terms): EggNog based on orthology and phylogenies

    EggNOG-mapper assigns both KEGG Orthology (KO) IDs and GO terms.

    Install EggNOG-mapper:

        mamba create -n eggnog_env python=3.8 eggnog-mapper -c conda-forge -c bioconda  #eggnog-mapper_2.1.12
        mamba activate eggnog_env

    Run annotation:

        #diamond makedb --in eggnog6.prots.faa -d eggnog_proteins.dmnd
        mkdir /home/jhuang/mambaforge/envs/eggnog_env/lib/python3.8/site-packages/data/
        download_eggnog_data.py --dbname eggnog.db -y --data_dir /home/jhuang/mambaforge/envs/eggnog_env/lib/python3.8/site-packages/data/
        #NOT_WORKING: emapper.py -i CP059040_gene.fasta -o eggnog_dmnd_out --cpu 60 -m diamond[hmmer,mmseqs] --dmnd_db /home/jhuang/REFs/eggnog_data/data/eggnog_proteins.dmnd
        #Download CU459141_protein_.fasta from NCBI
        python ~/Scripts/update_fasta_header.py CU459141_protein_.fasta CU459141_protein.fasta
        emapper.py -i CU459141_protein.fasta -o eggnog_out --cpu 60 --resume
        #----> result annotations.tsv: Contains KEGG, GO, and other functional annotations.
        #---->  470.IX87_14445:
            * 470 likely refers to the organism or strain (e.g., Acinetobacter baumannii ATCC 19606 or another related strain).
            * IX87_14445 would refer to a specific gene or protein within that genome.

    Extract KEGG KO IDs from annotations.emapper.annotations.

* Preparing file 2 blast2go_annot.annot2_ for the R-code below:

  - Basic (GO Terms from 'Blast2GO 5 Basic', saved in blast2go_annot.annot): Using Blast/Diamond + Blast2GO_GUI based on sequence alignment + GO mapping

    * 'Load protein sequences' (Tags: NONE, generated columns: Nr, SeqName) -->
    * Buttons 'blast' (Tags: BLASTED, generated columns: Description, Length, #Hits, e-Value, sim mean),
    * Button 'mapping' (Tags: MAPPED, generated columns: #GO, GO IDs, GO Names), "Mapping finished - Please proceed now to annotation."
    * Button 'annot' (Tags: ANNOTATED, generated columns: Enzyme Codes, Enzyme Names), "Annotation finished."
            * Used parameter 'Annotation CutOff': The Blast2GO Annotation Rule seeks to find the most specific GO annotations with a certain level of reliability. An annotation score is calculated for each candidate GO which is composed by the sequence similarity of the Blast Hit, the evidence code of the source GO and the position of the particular GO in the Gene Ontology hierarchy. This annotation score cutoff select the most specific GO term for a given GO branch which lies above this value.
            * Used parameter 'GO Weight' is a value which is added to Annotation Score of a more general/abstract Gene Ontology term for each of its more specific, original source GO terms. In this case, more general GO terms which summarise many original source terms (those ones directly associated to the Blast Hits) will have a higher Annotation Score.

  - Advanced (GO Terms from 'Blast2GO 5 Basic'): Interpro based protein families / domains --> Button interpro

    * Button 'interpro' (Tags: INTERPRO, generated columns: InterPro IDs, InterPro GO IDs, InterPro GO Names) --> "InterProScan Finished - You can now merge the obtained GO Annotations."

  - MERGE the results of InterPro GO IDs (advanced) to GO IDs (basic) and generate final GO IDs, saved in blast2go_annot.annot2

    * Button 'interpro'/'Merge InterProScan GOs to Annotation' --> "Merge (add and validate) all GO terms retrieved via InterProScan to the already existing GO annotation." --> "Finished merging GO terms from InterPro with annotations. Maybe you want to run ANNEX (Annotation Augmentation)."
    * (NOT_USED) Button 'annot'/'ANNEX' --> "ANNEX finished. Maybe you want to do the next step: Enzyme Code Mapping."

  - PREPARING go_terms and ec_terms: annot_* file (NOTE that blast2go_annot.annot2 is after merging InterPro_GO_IDs and GO_IDs):

    cut -f1-2 -d$'\t' blast2go_annot.annot2 > blast2go_annot.annot2_

10.2. Perform KEGG and GO Enrichment in R

      (r_env) cd /mnt/md1/DATA/Data_Tam_RNAseq_2026_Dicl_Mero_Azith_Rifa_on_AYE/results/star_salmon/DEG_Results_Complete

      #For |deg_cutoff_log_foldchange| >=1.4
      sed "s/01_AYE_T_ctr_vs_AYE_WT_ctr/02_AYE_O_ctr_vs_AYE_WT_ctr/g" 1.R > 2.R
      ...

      #For |deg_cutoff_log_foldchange| >=2.0
      sed "s/08_AYE_WT_ctr_solid_vs_liquid/09_AYE_O_ctr_solid_vs_liquid/g" 8.R > 9.R
      sed "s/08_AYE_WT_ctr_solid_vs_liquid/10_AYE_T_ctr_solid_vs_liquid/g" 8.R > 10.R
      sed "s/08_AYE_WT_ctr_solid_vs_liquid/11_AYE_O_ctr_solid_vs_AYE_WT_solid/g" 8.R > 11.R
      sed "s/08_AYE_WT_ctr_solid_vs_liquid/12_AYE_T_ctr_solid_vs_AYE_WT_solid/g" 8.R > 12.R
      sed "s/08_AYE_WT_ctr_solid_vs_liquid/13_AYE_WT_Diclo750_vs_Ctrl/g" 8.R > 13.R
      sed "s/08_AYE_WT_ctr_solid_vs_liquid/14_AYE_T_Diclo375_vs_Ctrl/g" 8.R > 14.R
      sed "s/08_AYE_WT_ctr_solid_vs_liquid/15_AYE_O_Diclo375_vs_Ctrl/g" 8.R > 15.R
      sed "s/08_AYE_WT_ctr_solid_vs_liquid/16_O_Trans_Diclo375_vs_Ctrl/g" 8.R > 16.R
      sed "s/08_AYE_WT_ctr_solid_vs_liquid/17_WT_Trans_Diclo750_vs_Ctrl/g" 8.R > 17.R
      sed "s/08_AYE_WT_ctr_solid_vs_liquid/18_AYE_WT_Diclo1250_solid_vs_Ctrl_solid/g" 8.R > 18.R
      sed "s/08_AYE_WT_ctr_solid_vs_liquid/19_AYE_WT_Mero_vs_Ctrl/g" 8.R > 19.R
      sed "s/08_AYE_WT_ctr_solid_vs_liquid/20_AYE_T_Mero_vs_Ctrl/g" 8.R > 20.R
      sed "s/08_AYE_WT_ctr_solid_vs_liquid/21_AYE_O_Mero_vs_Ctrl/g" 8.R > 21.R
      sed "s/08_AYE_WT_ctr_solid_vs_liquid/22_O_Trans_Mero_vs_Ctrl/g" 8.R > 22.R
      sed "s/08_AYE_WT_ctr_solid_vs_liquid/23_AYE_T_Mero_vs_AYE_T_Ctrl/g" 8.R > 23.R
      sed "s/08_AYE_WT_ctr_solid_vs_liquid/24_AYE_WT_Azi_solid_vs_Ctrl_solid/g" 8.R > 24.R
      sed "s/08_AYE_WT_ctr_solid_vs_liquid/25_AYE_T_Azi_solid_vs_Ctrl_solid/g" 8.R > 25.R
      sed "s/08_AYE_WT_ctr_solid_vs_liquid/26_AYE_O_Azi_solid_vs_Ctrl_solid/g" 8.R > 26.R
      sed "s/08_AYE_WT_ctr_solid_vs_liquid/27_F_Azi_solid_vs_Ctrl_solid/g" 8.R > 27.R

      #For |deg_cutoff_log_foldchange| >=1.2
      sed "s/28_AYE_WT_Rif_vs_Ctrl/29_AYE_T_Rif_vs_Ctrl/g" 28.R > 29.R
      sed "s/28_AYE_WT_Rif_vs_Ctrl/30_AYE_O_Rif_vs_Ctrl/g" 28.R > 30.R
      sed "s/28_AYE_WT_Rif_vs_Ctrl/31_O_Trans_Rif_vs_Ctrl/g" 28.R > 31.R

      (r_env) jhuang@WS-2290C:/mnt/md1/DATA/Data_Tam_RNAseq_2026_Dicl_Mero_Azith_Rifa_on_AYE/results/star_salmon/DEG_Results_Complete$ Rscript 1.R
      #=== SUMMARY ===
      #Up-regulated genes: 16
      #  Valid KEGG IDs: 4
      #  Enriched pathways: 0
      #Down-regulated genes: 151
      #  Valid KEGG IDs: 50
      #  Enriched pathways: 4
      #'select()' returned 1:1 mapping between keys and columns
      #'select()' returned 1:1 mapping between keys and columns
      #'select()' returned 1:1 mapping between keys and columns
      #=== SUMMARY ===
      #Up-regulated genes: 16
      #  Valid GO IDs: 16
      #  Enriched GO-terms: 0
      #Down-regulated genes: 151
      #  Valid KEGG IDs: 151
      #  Enriched GO-terms: 3
      #...

10.3. Finalizing the KEGG and GO Enrichment table

      1. NOTE (Already realized in the code): geneIDs in KEGG_Enrichment have been already translated from ko to geneID in H0N29_*-format; If not, nachmachen using eggnog-res, 因为 eggnog里有1-1-mspping Info between ko-Name and GeneID.
      2. NEED_MANUAL_DELETION (Already setting the cutoff in the code): p.adjust values have been calculated, we have to filter all records in GO_Enrichment-results by |p.adjust|<=0.05. DON'T_NEED_ANY_MORE, since pvalueCutoff = 0.05 settings in enricher. Alternative using pvalueCutoff=1.0, marked the color as yellow if the p.adjusted <= 0.05 in GO_enrichment.
      3. NOTE (Not occuring in the new dataset): In rare case, the description is missing for some IDs, e.g. GO term: GO:0006807: replace GO:0006807    obsolete nitrogen compound metabolic process;  ko00975: Metabolism, Biosynthesis of other secondary metabolites

KEGG and GO enrichments via 1.R, 2.R, …

  #Update "01_AYE_T_ctr_vs_AYE_WT_ctr" with "02_AYE_O_ctr_vs_AYE_WT_ctr" in 1.R and save the updated script as 2.R.
  (r_env) jhuang@WS-2290C:/mnt/md1/DATA/Data_Tam_RNAseq_2026_Dicl_Mero_Azith_Rifa_on_AYE/results/star_salmon/DEG_Results_Complete$ Rscript 2.R
  === SUMMARY ===
  Up-regulated genes: 9
    Valid KEGG IDs: 1
    Enriched pathways: 0
  Down-regulated genes: 9
    Valid KEGG IDs: 6
    Enriched pathways: 2
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  === SUMMARY ===
  Up-regulated genes: 9
    Valid GO IDs: 9
    Enriched GO-terms: 0
  Down-regulated genes: 9
    Valid KEGG IDs: 9
    Enriched GO-terms: 2

  Rscript 3.R
  === SUMMARY ===
  Up-regulated genes: 10
    Valid KEGG IDs: 2
    Enriched pathways: 0
  Down-regulated genes: 44
    Valid KEGG IDs: 20
    Enriched pathways: 6
  === SUMMARY ===
  Up-regulated genes: 10
    Valid GO IDs: 10
    Enriched GO-terms: 0
  Down-regulated genes: 44
    Valid KEGG IDs: 44
    Enriched GO-terms: 0

  Rscript 4.R
  === SUMMARY ===
  Up-regulated genes: 15
    Valid KEGG IDs: 7
    Enriched pathways: 6
  Down-regulated genes: 152
    Valid KEGG IDs: 43
    Enriched pathways: 11
  === SUMMARY ===
  Up-regulated genes: 15
    Valid GO IDs: 15
    Enriched GO-terms: 0
  Down-regulated genes: 152
    Valid KEGG IDs: 152
    Enriched GO-terms: 0

  Rscript 5.R
  === SUMMARY ===
  Up-regulated genes: 28
    Valid KEGG IDs: 8
    Enriched pathways: 4
  Down-regulated genes: 3
    Valid KEGG IDs: 2
    Enriched pathways: 2
  'select()' returned 1:1 mapping between keys and columns
  === SUMMARY ===
  Up-regulated genes: 28
    Valid GO IDs: 28
    Enriched GO-terms: 0
  Down-regulated genes: 3
    Valid KEGG IDs: 3
    Enriched GO-terms: 1

  Rscript 6.R
  === SUMMARY ===
  Up-regulated genes: 28
    Valid KEGG IDs: 8
    Enriched pathways: 4
  Down-regulated genes: 3
    Valid KEGG IDs: 2
    Enriched pathways: 2
  'select()' returned 1:1 mapping between keys and columns
  === SUMMARY ===
  Up-regulated genes: 28
    Valid GO IDs: 28
    Enriched GO-terms: 0
  Down-regulated genes: 3
    Valid KEGG IDs: 3
    Enriched GO-terms: 1

  Rscript 7.R
  === SUMMARY ===
  Up-regulated genes: 4
    Valid KEGG IDs: 0
    Enriched pathways: 0
  Down-regulated genes: 14
    Valid KEGG IDs: 10
    Enriched pathways: 1
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  === SUMMARY ===
  Up-regulated genes: 4
    Valid GO IDs: 4
    Enriched GO-terms: 1
  Down-regulated genes: 14
    Valid KEGG IDs: 14
    Enriched GO-terms: 3

  Rscript 8.R
  === SUMMARY ===
  Up-regulated genes: 479
    Valid KEGG IDs: 246
    Enriched pathways: 30
  Down-regulated genes: 322
    Valid KEGG IDs: 186
    Enriched pathways: 13
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  === SUMMARY ===
  Up-regulated genes: 479
    Valid GO IDs: 479
    Enriched GO-terms: 4
  Down-regulated genes: 322
    Valid KEGG IDs: 322
    Enriched GO-terms: 6

  Rscript 9.R
  === SUMMARY ===
  Up-regulated genes: 472
    Valid KEGG IDs: 230
    Enriched pathways: 25
  Down-regulated genes: 306
    Valid KEGG IDs: 217
    Enriched pathways: 19
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  === SUMMARY ===
  Up-regulated genes: 472
    Valid GO IDs: 472
    Enriched GO-terms: 3
  Down-regulated genes: 306
    Valid KEGG IDs: 306
    Enriched GO-terms: 6

  Rscript 10.R
  === SUMMARY ===
  Up-regulated genes: 496
    Valid KEGG IDs: 236
    Enriched pathways: 28
  Down-regulated genes: 330
    Valid KEGG IDs: 250
    Enriched pathways: 26
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  === SUMMARY ===
  Up-regulated genes: 496
    Valid GO IDs: 496
    Enriched GO-terms: 2
  Down-regulated genes: 330
    Valid KEGG IDs: 330
    Enriched GO-terms: 6

  Rscript 11.R
  === SUMMARY ===
  Up-regulated genes: 12
    Valid KEGG IDs: 2
    Enriched pathways: 5
  Down-regulated genes: 5
    Valid KEGG IDs: 2
    Enriched pathways: 2
  === SUMMARY ===
  Up-regulated genes: 12
    Valid GO IDs: 12
    Enriched GO-terms: 0
  Down-regulated genes: 5
    Valid KEGG IDs: 5
    Enriched GO-terms: 0

  Rscript 12.R
  === SUMMARY ===
  Up-regulated genes: 5
    Valid KEGG IDs: 1
    Enriched pathways: 3
  Down-regulated genes: 8
    Valid KEGG IDs: 3
    Enriched pathways: 2
  === SUMMARY ===
  Up-regulated genes: 5
    Valid GO IDs: 5
    Enriched GO-terms: 0
  Down-regulated genes: 8
    Valid KEGG IDs: 8
    Enriched GO-terms: 0

  Rscript 13.R
  === SUMMARY ===
  Up-regulated genes: 12
    Valid KEGG IDs: 8
    Enriched pathways: 4
  Down-regulated genes: 45
    Valid KEGG IDs: 13
    Enriched pathways: 4
  'select()' returned 1:1 mapping between keys and columns
  === SUMMARY ===
  Up-regulated genes: 12
    Valid GO IDs: 12
    Enriched GO-terms: 1
  Down-regulated genes: 45
    Valid KEGG IDs: 45
    Enriched GO-terms: 0

  Rscript 14.R
  === SUMMARY ===
  Up-regulated genes: 73
    Valid KEGG IDs: 37
    Enriched pathways: 5
  Down-regulated genes: 180
    Valid KEGG IDs: 29
    Enriched pathways: 3
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  === SUMMARY ===
  Up-regulated genes: 73
    Valid GO IDs: 73
    Enriched GO-terms: 5
  Down-regulated genes: 180
    Valid KEGG IDs: 180
    Enriched GO-terms: 0

  Rscript 15.R
  === SUMMARY ===
  Up-regulated genes: 26
    Valid KEGG IDs: 14
    Enriched pathways: 10
  Down-regulated genes: 239
    Valid KEGG IDs: 38
    Enriched pathways: 3
  === SUMMARY ===
  Up-regulated genes: 26
    Valid GO IDs: 26
    Enriched GO-terms: 0
  Down-regulated genes: 239
    Valid KEGG IDs: 239
    Enriched GO-terms: 0

  Rscript 16.R
  === SUMMARY ===
  Up-regulated genes: 17
    Valid KEGG IDs: 8
    Enriched pathways: 11
  Down-regulated genes: 82
    Valid KEGG IDs: 14
    Enriched pathways: 3
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  === SUMMARY ===
  Up-regulated genes: 17
    Valid GO IDs: 17
    Enriched GO-terms: 0
  Down-regulated genes: 82
    Valid KEGG IDs: 82
    Enriched GO-terms: 2

  Rscript 17.R
  === SUMMARY ===
  Up-regulated genes: 25
    Valid KEGG IDs: 15
    Enriched pathways: 3
  Down-regulated genes: 229
    Valid KEGG IDs: 46
    Enriched pathways: 3
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  === SUMMARY ===
  Up-regulated genes: 25
    Valid GO IDs: 25
    Enriched GO-terms: 2
  Down-regulated genes: 229
    Valid KEGG IDs: 229
    Enriched GO-terms: 1

  Rscript 18.R
  === SUMMARY ===
  Up-regulated genes: 66
    Valid KEGG IDs: 18
    Enriched pathways: 8
  Down-regulated genes: 26
    Valid KEGG IDs: 19
    Enriched pathways: 4
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  === SUMMARY ===
  Up-regulated genes: 66
    Valid GO IDs: 66
    Enriched GO-terms: 5
  Down-regulated genes: 26
    Valid KEGG IDs: 26
    Enriched GO-terms: 3

  Rscript 19.R
  === SUMMARY ===
  Up-regulated genes: 355
    Valid KEGG IDs: 200
    Enriched pathways: 25
  Down-regulated genes: 175
    Valid KEGG IDs: 70
    Enriched pathways: 3
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  === SUMMARY ===
  Up-regulated genes: 355
    Valid GO IDs: 355
    Enriched GO-terms: 6
  Down-regulated genes: 175
    Valid KEGG IDs: 175
    Enriched GO-terms: 2

  Rscript 20.R
  === SUMMARY ===
  Up-regulated genes: 93
    Valid KEGG IDs: 48
    Enriched pathways: 21
  Down-regulated genes: 110
    Valid KEGG IDs: 32
    Enriched pathways: 4
  === SUMMARY ===
  Up-regulated genes: 93
    Valid GO IDs: 93
    Enriched GO-terms: 0
  Down-regulated genes: 110
    Valid KEGG IDs: 110
    Enriched GO-terms: 0

  Rscript 21.R
  === SUMMARY ===
  Up-regulated genes: 263
    Valid KEGG IDs: 154
    Enriched pathways: 26
  Down-regulated genes: 244
    Valid KEGG IDs: 72
    Enriched pathways: 3
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  === SUMMARY ===
  Up-regulated genes: 263
    Valid GO IDs: 263
    Enriched GO-terms: 8
  Down-regulated genes: 244
    Valid KEGG IDs: 244
    Enriched GO-terms: 5

  Rscript 22.R
  === SUMMARY ===
  Up-regulated genes: 182
    Valid KEGG IDs: 115
    Enriched pathways: 23
  Down-regulated genes: 153
    Valid KEGG IDs: 40
    Enriched pathways: 3
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  === SUMMARY ===
  Up-regulated genes: 182
    Valid GO IDs: 182
    Enriched GO-terms: 3
  Down-regulated genes: 153
    Valid KEGG IDs: 153
    Enriched GO-terms: 2

  Rscript 23.R
  === SUMMARY ===
  Up-regulated genes: 91
    Valid KEGG IDs: 37
    Enriched pathways: 17
  Down-regulated genes: 27
    Valid KEGG IDs: 14
    Enriched pathways: 1
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  === SUMMARY ===
  Up-regulated genes: 91
    Valid GO IDs: 91
    Enriched GO-terms: 1
  Down-regulated genes: 27
    Valid KEGG IDs: 27
    Enriched GO-terms: 4

  Rscript 24.R
  === SUMMARY ===
  Up-regulated genes: 333
    Valid KEGG IDs: 193
    Enriched pathways: 12
  Down-regulated genes: 328
    Valid KEGG IDs: 175
    Enriched pathways: 24
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  === SUMMARY ===
  Up-regulated genes: 333
    Valid GO IDs: 333
    Enriched GO-terms: 4
  Down-regulated genes: 328
    Valid KEGG IDs: 328
    Enriched GO-terms: 4

  Rscript 25.R
  === SUMMARY ===
  Up-regulated genes: 499
    Valid KEGG IDs: 264
    Enriched pathways: 19
  Down-regulated genes: 430
    Valid KEGG IDs: 239
    Enriched pathways: 36
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  === SUMMARY ===
  Up-regulated genes: 499
    Valid GO IDs: 499
    Enriched GO-terms: 2
  Down-regulated genes: 430
    Valid KEGG IDs: 430
    Enriched GO-terms: 0

  Rscript 26.R
  === SUMMARY ===
  Up-regulated genes: 474
    Valid KEGG IDs: 267
    Enriched pathways: 19
  Down-regulated genes: 440
    Valid KEGG IDs: 216
    Enriched pathways: 32
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  === SUMMARY ===
  Up-regulated genes: 474
    Valid GO IDs: 474
    Enriched GO-terms: 2
  Down-regulated genes: 440
    Valid KEGG IDs: 440
    Enriched GO-terms: 1

  Rscript 27.R
  === SUMMARY ===
  Up-regulated genes: 142
    Valid KEGG IDs: 69
    Enriched pathways: 1
  Down-regulated genes: 269
    Valid KEGG IDs: 148
    Enriched pathways: 27
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  === SUMMARY ===
  Up-regulated genes: 142
    Valid GO IDs: 142
    Enriched GO-terms: 5
  Down-regulated genes: 269
    Valid KEGG IDs: 269
    Enriched GO-terms: 3

  Rscript 28.R
  === SUMMARY ===
  Up-regulated genes: 8
    Valid KEGG IDs: 6
    Enriched pathways: 2
  Down-regulated genes: 11
    Valid KEGG IDs: 3
    Enriched pathways: 0
  === SUMMARY ===
  Up-regulated genes: 8
    Valid GO IDs: 8
    Enriched GO-terms: 0
  Down-regulated genes: 11
    Valid KEGG IDs: 11
    Enriched GO-terms: 0

  Rscript 29.R
  === SUMMARY ===
  Up-regulated genes: 34
    Valid KEGG IDs: 11
    Enriched pathways: 6
  Down-regulated genes: 41
    Valid KEGG IDs: 22
    Enriched pathways: 1
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  === SUMMARY ===
  Up-regulated genes: 34
    Valid GO IDs: 34
    Enriched GO-terms: 0
  Down-regulated genes: 41
    Valid KEGG IDs: 41

  Rscript 30.R
  === SUMMARY ===
  Up-regulated genes: 70
    Valid KEGG IDs: 47
    Enriched pathways: 16
  Down-regulated genes: 99
    Valid KEGG IDs: 40
    Enriched pathways: 3
  === SUMMARY ===
  Up-regulated genes: 70
    Valid GO IDs: 70
    Enriched GO-terms: 0
  Down-regulated genes: 99
    Valid KEGG IDs: 99
    Enriched GO-terms: 0

  Rscript 31.R
  === SUMMARY ===
  Up-regulated genes: 6
    Valid KEGG IDs: 3
    Enriched pathways: 0
  Down-regulated genes: 10
    Valid KEGG IDs: 5
    Enriched pathways: 1
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  === SUMMARY ===
  Up-regulated genes: 6
    Valid GO IDs: 6
    Enriched GO-terms: 0
  Down-regulated genes: 10
    Valid KEGG IDs: 10
    Enriched GO-terms: 5

  Rscript 32.R
  === SUMMARY ===
  Up-regulated genes: 46
    Valid KEGG IDs: 18
    Enriched pathways: 1
  Down-regulated genes: 6
    Valid KEGG IDs: 0
    Enriched pathways: 0
  No gene sets have size between 10 and 500 ...
  --> return NULL...
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns

  === SUMMARY ===
  Up-regulated genes: 46
    Valid GO IDs: 46
    Enriched GO-terms: 4
  Down-regulated genes: 6
    Valid KEGG IDs: 6
    Enriched GO-terms: 0

  Rscript 33.R
  === SUMMARY ===
  Up-regulated genes: 8
    Valid KEGG IDs: 9
    Enriched pathways: 8
  Down-regulated genes: 72
    Valid KEGG IDs: 7
    Enriched pathways: 3

  === SUMMARY ===
  Up-regulated genes: 8
    Valid GO IDs: 8
    Enriched GO-terms: 0
  Down-regulated genes: 72
    Valid KEGG IDs: 72
    Enriched GO-terms: 0

  Rscript 34.R
  === SUMMARY ===
  Up-regulated genes: 2
    Valid KEGG IDs: 2
    Enriched pathways: 3
  Down-regulated genes: 17
    Valid KEGG IDs: 3
    Enriched pathways: 0

  === SUMMARY ===
  Up-regulated genes: 2
    Valid GO IDs: 2
    Enriched GO-terms: 0
  Down-regulated genes: 17
    Valid KEGG IDs: 17
    Enriched GO-terms: 0

  Rscript 35.R
  === SUMMARY ===
  Up-regulated genes: 10
    Valid KEGG IDs: 5
    Enriched pathways: 5
  Down-regulated genes: 27
    Valid KEGG IDs: 3
    Enriched pathways: 0

  === SUMMARY ===
  Up-regulated genes: 10
    Valid GO IDs: 10
    Enriched GO-terms: 0
  Down-regulated genes: 27
    Valid KEGG IDs: 27
    Enriched GO-terms: 0

  Rscript 36.R
  === SUMMARY ===
  Up-regulated genes: 91
    Valid KEGG IDs: 37
    Enriched pathways: 17
  Down-regulated genes: 27
    Valid KEGG IDs: 14
    Enriched pathways: 1
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns

  === SUMMARY ===
  Up-regulated genes: 91
    Valid GO IDs: 91
    Enriched GO-terms: 1
  Down-regulated genes: 27
    Valid KEGG IDs: 27
    Enriched GO-terms: 4

  Rscript 37.R
  === SUMMARY ===
  Up-regulated genes: 238
    Valid KEGG IDs: 141
    Enriched pathways: 28
  Down-regulated genes: 134
    Valid KEGG IDs: 35
    Enriched pathways: 2
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns

  === SUMMARY ===
  Up-regulated genes: 238
    Valid GO IDs: 238
    Enriched GO-terms: 2
  Down-regulated genes: 134
    Valid KEGG IDs: 134
    Enriched GO-terms: 4

  Rscript 38.R
  === SUMMARY ===
  Up-regulated genes: 177
    Valid KEGG IDs: 104
    Enriched pathways: 22
  Down-regulated genes: 50
    Valid KEGG IDs: 13
    Enriched pathways: 4
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns
  'select()' returned 1:1 mapping between keys and columns

  === SUMMARY ===
  Up-regulated genes: 177
    Valid GO IDs: 177
    Enriched GO-terms: 2
  Down-regulated genes: 50
    Valid KEGG IDs: 50
    Enriched GO-terms: 4

(DEPRECATED due to be WRONG and TIME-CONSUMING) KEGG and GO enrichments by replacing the sample name in the R-code, and by copying and pasting the R-code directly on console

  # For your reference, here is the exact list of the 31 comparisons and their assigned groupings:

  # Baseline / Strain Controls (|log2FC| >= 1.4)
  "01_AYE_T_ctr_vs_AYE_WT_ctr"       = c("group", "AYE-T_ctrl_0_liquid",   "AYE-WT_ctrl_0_liquid"),
  "02_AYE_O_ctr_vs_AYE_WT_ctr"       = c("group", "AYE-O_ctrl_0_liquid",   "AYE-WT_ctrl_0_liquid"),
  "03_O_Trans_ctr_vs_AYE_WT_ctr"     = c("group", "O-Trans_ctrl_0_liquid", "AYE-WT_ctrl_0_liquid"),
  "04_WT_Trans_ctr_vs_AYE_WT_ctr"    = c("group", "WT-Trans_ctrl_0_liquid","AYE-WT_ctrl_0_liquid"),
  "05_AYE_O_ctr_vs_AYE_T_ctr"        = c("group", "AYE-O_ctrl_0_liquid",   "AYE-T_ctrl_0_liquid"),
  "06_O_Trans_ctr_vs_AYE_T_ctr"      = c("group", "O-Trans_ctrl_0_liquid", "AYE-T_ctrl_0_liquid"),
  "07_WT_Trans_ctr_vs_AYE_T_ctr"     = c("group", "WT-Trans_ctrl_0_liquid","AYE-T_ctrl_0_liquid"),

  # Condition Effects (|log2FC| >= 2.0)
  "08_AYE_WT_ctr_solid_vs_liquid"    = c("group", "AYE-WT_ctrl_0_solid",   "AYE-WT_ctrl_0_liquid"),
  "09_AYE_O_ctr_solid_vs_liquid"     = c("group", "AYE-O_ctrl_0_solid",    "AYE-O_ctrl_0_liquid"),
  "10_AYE_T_ctr_solid_vs_liquid"     = c("group", "AYE-T_ctrl_0_solid",    "AYE-T_ctrl_0_liquid"),
  "11_AYE_O_ctr_solid_vs_AYE_WT_solid"=c("group", "AYE-O_ctrl_0_solid",    "AYE-WT_ctrl_0_solid"),
  "12_AYE_T_ctr_solid_vs_AYE_WT_solid"=c("group", "AYE-T_ctrl_0_solid",    "AYE-WT_ctrl_0_solid"),

  # Diclofenac (|log2FC| >= 2.0)
  "13_AYE_WT_Diclo750_vs_Ctrl"       = c("group", "AYE-WT_Diclo_750_liquid",   "AYE-WT_ctrl_0_liquid"),
  "14_AYE_T_Diclo375_vs_Ctrl"        = c("group", "AYE-T_Diclo_375_liquid",    "AYE-WT_ctrl_0_liquid"),
  "15_AYE_O_Diclo375_vs_Ctrl"        = c("group", "AYE-O_Diclo_375_liquid",    "AYE-WT_ctrl_0_liquid"),
  "16_O_Trans_Diclo375_vs_Ctrl"      = c("group", "O-Trans_Diclo_375_liquid",  "AYE-WT_ctrl_0_liquid"),
  "17_WT_Trans_Diclo750_vs_Ctrl"     = c("group", "WT-Trans_Diclo_750_liquid", "AYE-WT_ctrl_0_liquid"),
  "18_AYE_WT_Diclo1250_solid_vs_Ctrl_solid" = c("group", "AYE-WT_Diclo_1250_solid", "AYE-WT_ctrl_0_solid"),

  # Meropenem (|log2FC| >= 2.0)
  "19_AYE_WT_Mero_vs_Ctrl"           = c("group", "AYE-WT_Mero_0.35-0.5_liquid", "AYE-WT_ctrl_0_liquid"),
  "20_AYE_T_Mero_vs_Ctrl"            = c("group", "AYE-T_Mero_0.15_liquid",      "AYE-WT_ctrl_0_liquid"),
  "21_AYE_O_Mero_vs_Ctrl"            = c("group", "AYE-O_Mero_0.5_liquid",       "AYE-WT_ctrl_0_liquid"),
  "22_O_Trans_Mero_vs_Ctrl"          = c("group", "O-Trans_Mero_0.25_liquid",    "AYE-WT_ctrl_0_liquid"),
  "23_AYE_T_Mero_vs_AYE_T_Ctrl"      = c("group", "AYE-T_Mero_0.15_liquid",      "AYE-T_ctrl_0_liquid"),

  # Azithromycin (Solid) (|log2FC| >= 2.0)
  "24_AYE_WT_Azi_solid_vs_Ctrl_solid"= c("group", "AYE-WT_Azi_20_solid", "AYE-WT_ctrl_0_solid"),
  "25_AYE_T_Azi_solid_vs_Ctrl_solid" = c("group", "AYE-T_Azi_20_solid",  "AYE-T_ctrl_0_solid"),
  "26_AYE_O_Azi_solid_vs_Ctrl_solid" = c("group", "AYE-O_Azi_20_solid",  "AYE-O_ctrl_0_solid"),
  "27_F_Azi_solid_vs_Ctrl_solid"     = c("group", "F_Azi_20_solid",      "F_ctrl_0_solid"),

  # Rifampicin (|log2FC| >= 1.2)
  "28_AYE_WT_Rif_vs_Ctrl"            = c("group", "AYE-WT_Rifampicin_1.5_liquid", "AYE-WT_ctrl_0_liquid"),
  "29_AYE_T_Rif_vs_Ctrl"             = c("group", "AYE-T_Rifampicin_2_liquid",    "AYE-T_ctrl_0_liquid"),
  "30_AYE_O_Rif_vs_Ctrl"             = c("group", "AYE-O_Rifampicin_2_liquid",    "AYE-O_ctrl_0_liquid"),
  "31_O_Trans_Rif_vs_Ctrl"           = c("group", "O-Trans_Rifampicin_2_liquid",  "O-Trans_ctrl_0_liquid")

  ./DEG_01_AYE_T_ctr_vs_AYE_WT_ctr.csv
  ./DEG_02_AYE_O_ctr_vs_AYE_WT_ctr.csv
  ./DEG_03_O_Trans_ctr_vs_AYE_WT_ctr.csv
  ./DEG_04_WT_Trans_ctr_vs_AYE_WT_ctr.csv
  ./DEG_05_AYE_O_ctr_vs_AYE_T_ctr.csv
  ./DEG_06_O_Trans_ctr_vs_AYE_T_ctr.csv
  ./DEG_07_WT_Trans_ctr_vs_AYE_T_ctr.csv
  ./DEG_08_AYE_WT_ctr_solid_vs_liquid.csv
  ./DEG_09_AYE_O_ctr_solid_vs_liquid.csv
  ./DEG_10_AYE_T_ctr_solid_vs_liquid.csv
  ./DEG_11_AYE_O_ctr_solid_vs_AYE_WT_solid.csv
  ./DEG_12_AYE_T_ctr_solid_vs_AYE_WT_solid.csv
  ./DEG_13_AYE_WT_Diclo750_vs_Ctrl.csv
  ./DEG_14_AYE_T_Diclo375_vs_Ctrl.csv
  ./DEG_15_AYE_O_Diclo375_vs_Ctrl.csv
  ./DEG_16_O_Trans_Diclo375_vs_Ctrl.csv
  ./DEG_17_WT_Trans_Diclo750_vs_Ctrl.csv
  ./DEG_18_AYE_WT_Diclo1250_solid_vs_Ctrl_solid.csv
  ./DEG_19_AYE_WT_Mero_vs_Ctrl.csv
  ./DEG_20_AYE_T_Mero_vs_Ctrl.csv
  ./DEG_21_AYE_O_Mero_vs_Ctrl.csv
  ./DEG_22_O_Trans_Mero_vs_Ctrl.csv
  ./DEG_23_AYE_T_Mero_vs_AYE_T_Ctrl.csv
  ./DEG_24_AYE_WT_Azi_solid_vs_Ctrl_solid.csv
  ./DEG_25_AYE_T_Azi_solid_vs_Ctrl_solid.csv
  ./DEG_26_AYE_O_Azi_solid_vs_Ctrl_solid.csv
  ./DEG_27_F_Azi_solid_vs_Ctrl_solid.csv
  ./DEG_28_AYE_WT_Rif_vs_Ctrl.csv
  ./DEG_29_AYE_T_Rif_vs_Ctrl.csv
  ./DEG_30_AYE_O_Rif_vs_Ctrl.csv
  ./DEG_31_O_Trans_Rif_vs_Ctrl.csv

    #BiocManager::install("GO.db")
    #BiocManager::install("AnnotationDbi")

    # Load required libraries
    library(openxlsx)  # For Excel file handling
    library(dplyr)     # For data manipulation
    library(tidyr)
    library(stringr)
    library(clusterProfiler)  # For KEGG and GO enrichment analysis
    #library(org.Hs.eg.db)  # Replace with appropriate organism database
    library(GO.db)
    library(AnnotationDbi)

    setwd("~/DATA/Data_Tam_RNAseq_2026_Dicl_Mero_Azith_Rifa_on_AYE/results/star_salmon/DEG_Results_Complete_manually")

    # 1. Blast2GO: Extract GO & EC terms (Primary source)
    annot_df <- read.table("/home/jhuang/b2gWorkspace_Tam_RNAseq_AYE/blast2go_annot.annot2_",
                        header = FALSE, sep = "\t", stringsAsFactors = FALSE, fill = TRUE)
    colnames(annot_df) <- c("GeneID", "Term")

    go_terms <- annot_df %>%
    filter(grepl("^GO:", Term)) %>%
    group_by(GeneID) %>%
    summarize(GOs = paste(Term, collapse = ","), .groups = "drop")

    ec_terms <- annot_df %>%
    filter(grepl("^EC:", Term)) %>%
    group_by(GeneID) %>%
    summarize(EC = paste(Term, collapse = ","), .groups = "drop")

    # Load the results
    res <- read.csv("DEG_01_AYE_T_ctr_vs_AYE_WT_ctr.csv")

    # Replace empty GeneName with modified GeneID
    res$GeneName <- ifelse(
        res$GeneName == "" | is.na(res$GeneName),
        gsub("gene-", "", res$GeneID),
        res$GeneName
    )

    # Remove duplicated genes by selecting the gene with the smallest padj
    duplicated_genes <- res[duplicated(res$GeneName), "GeneName"]

    res <- res %>%
    group_by(GeneName) %>%
    slice_min(padj, with_ties = FALSE) %>%
    ungroup()

    res <- as.data.frame(res)
    # Sort res first by padj (ascending) and then by log2FoldChange (descending)
    res <- res[order(res$padj, -res$log2FoldChange), ]
    # Read eggnog annotations
    eggnog_data <- read.delim("~/DATA/Data_Tam_RNAseq_2026_Dicl_Mero_Azith_Rifa_on_AYE/eggnog_out.emapper.annotations.txt", header = TRUE, sep = "\t")
    # Remove the "gene-" prefix from GeneID in res to match eggnog 'query' format
    res$GeneID <- gsub("gene-", "", res$GeneID)
    # Merge eggnog data with res based on GeneID
    res <- res %>%
    left_join(eggnog_data, by = c("GeneID" = "query"))

    # Merge with the res dataframe
    # Perform the left joins and rename columns
    res_updated <- res %>%
    left_join(go_terms, by = "GeneID") %>%
    left_join(ec_terms, by = "GeneID") %>% dplyr::select(-EC.x, -GOs.x) %>% dplyr::rename(EC = EC.y, GOs = GOs.y)

    # Filter up- and down-regulated genes (NEEDS_TO_BE_ADAPTED_BY_1.4_or_2.0_or_1.2)
    up_regulated <- res_updated[res_updated$log2FoldChange >= 1.4 & res_updated$padj <= 0.05, ]
    down_regulated <- res_updated[res_updated$log2FoldChange <= -1.4 & res_updated$padj <= 0.05, ]

    # Create a new workbook
    wb <- createWorkbook()
    addWorksheet(wb, "Complete_Data")
    writeData(wb, "Complete_Data", res_updated)
    addWorksheet(wb, "Up_Regulated")
    writeData(wb, "Up_Regulated", up_regulated)
    addWorksheet(wb, "Down_Regulated")
    writeData(wb, "Down_Regulated", down_regulated)
    saveWorkbook(wb, "Gene_Expression_with_Annotations_01_AYE_T_ctr_vs_AYE_WT_ctr.xlsx", overwrite = TRUE)

    # Set GeneName as row names after the join
    rownames(res_updated) <- res_updated$GeneName
    res_updated <- res_updated %>% dplyr::select(-GeneName)

    # ---------------------------------------------------------
    # ---- Perform KEGG enrichment analysis (up_regulated) ----
    gene_list_kegg_up <- up_regulated$KEGG_ko
    gene_list_kegg_up <- gsub("ko:", "", gene_list_kegg_up)
    kegg_enrichment_up <- enrichKEGG(gene = gene_list_kegg_up, organism = 'ko')

    # -- convert the GeneID (Kxxxxxx) to the true GeneID --
    # Step 0: Create KEGG to GeneID mapping
    kegg_to_geneid_up <- up_regulated %>%
    dplyr::select(KEGG_ko, GeneID) %>%
    filter(!is.na(KEGG_ko)) %>%  # Remove missing KEGG KO entries
    mutate(KEGG_ko = str_remove(KEGG_ko, "ko:"))  # Remove 'ko:' prefix if present

    # Step 1: Clean KEGG_ko values (separate multiple KEGG IDs)
    kegg_to_geneid_clean <- kegg_to_geneid_up %>%
    mutate(KEGG_ko = str_remove_all(KEGG_ko, "ko:")) %>%  # Remove 'ko:' prefixes
    separate_rows(KEGG_ko, sep = ",") %>%  # Ensure each KEGG ID is on its own row
    filter(KEGG_ko != "-") %>%  # Remove invalid KEGG IDs ("-")
    distinct()  # Remove any duplicate mappings

    # Step 2.1: Expand geneID column in kegg_enrichment_up
    expanded_kegg <- kegg_enrichment_up %>%
    as.data.frame() %>%
    separate_rows(geneID, sep = "/") %>%  # Split multiple KEGG IDs (Kxxxxx)
    left_join(kegg_to_geneid_clean, by = c("geneID" = "KEGG_ko"), relationship = "many-to-many") %>%  # Explicitly handle many-to-many
    distinct() %>%  # Remove duplicate matches
    group_by(ID) %>%
    summarise(across(everything(), ~ paste(unique(na.omit(.)), collapse = "/")), .groups = "drop")  # Re-collapse results
    #dplyr::glimpse(expanded_kegg)

    # Step 3.1: Replace geneID column in the original dataframe
    kegg_enrichment_up_df <- as.data.frame(kegg_enrichment_up)
    # Remove old geneID column and merge new one
    kegg_enrichment_up_df <- kegg_enrichment_up_df %>%
    dplyr::select(-geneID) %>%  # Remove old geneID column
    left_join(expanded_kegg %>% dplyr::select(ID, GeneID), by = "ID") %>%  # Merge new GeneID column
    dplyr::rename(geneID = GeneID)  # Rename column back to geneID

    # -----------------------------------------------------------
    # ---- Perform KEGG enrichment analysis (down_regulated) ----
    # Step 1: Extract KEGG KO terms from down-regulated genes
    gene_list_kegg_down <- down_regulated$KEGG_ko
    gene_list_kegg_down <- gsub("ko:", "", gene_list_kegg_down)
    # Step 2: Perform KEGG enrichment analysis
    kegg_enrichment_down <- enrichKEGG(gene = gene_list_kegg_down, organism = 'ko')
    # --- Convert KEGG gene IDs (Kxxxxxx) to actual GeneIDs ---
    # Step 3: Create KEGG to GeneID mapping from down_regulated dataset
    kegg_to_geneid_down <- down_regulated %>%
    dplyr::select(KEGG_ko, GeneID) %>%
    filter(!is.na(KEGG_ko)) %>%  # Remove missing KEGG KO entries
    mutate(KEGG_ko = str_remove(KEGG_ko, "ko:"))  # Remove 'ko:' prefix if present
    # Step 4: Clean KEGG_ko values (handle multiple KEGG IDs)
    kegg_to_geneid_down_clean <- kegg_to_geneid_down %>%
    mutate(KEGG_ko = str_remove_all(KEGG_ko, "ko:")) %>%  # Remove 'ko:' prefixes
    separate_rows(KEGG_ko, sep = ",") %>%  # Ensure each KEGG ID is on its own row
    filter(KEGG_ko != "-") %>%  # Remove invalid KEGG IDs ("-")
    distinct()  # Remove duplicate mappings
    # Step 5: Expand geneID column in kegg_enrichment_down
    expanded_kegg_down <- kegg_enrichment_down %>%
    as.data.frame() %>%
    separate_rows(geneID, sep = "/") %>%  # Split multiple KEGG IDs (Kxxxxx)
    left_join(kegg_to_geneid_down_clean, by = c("geneID" = "KEGG_ko"), relationship = "many-to-many") %>%  # Handle many-to-many mappings
    distinct() %>%  # Remove duplicate matches
    group_by(ID) %>%
    summarise(across(everything(), ~ paste(unique(na.omit(.)), collapse = "/")), .groups = "drop")  # Re-collapse results
    # Step 6: Replace geneID column in the original kegg_enrichment_down dataframe
    kegg_enrichment_down_df <- as.data.frame(kegg_enrichment_down) %>%
    dplyr::select(-geneID) %>%  # Remove old geneID column
    left_join(expanded_kegg_down %>% dplyr::select(ID, GeneID), by = "ID") %>%  # Merge new GeneID column
    dplyr::rename(geneID = GeneID)  # Rename column back to geneID
    # View the updated dataframe
    head(kegg_enrichment_down_df)

    # Create a new workbook
    wb <- createWorkbook()
    # Save enrichment results to the workbook
    addWorksheet(wb, "KEGG_Enrichment_Up")
    writeData(wb, "KEGG_Enrichment_Up", as.data.frame(kegg_enrichment_up_df))
    # Save enrichment results to the workbook
    addWorksheet(wb, "KEGG_Enrichment_Down")
    writeData(wb, "KEGG_Enrichment_Down", as.data.frame(kegg_enrichment_down_df))
    #saveWorkbook(wb, "KEGG_Enrichment.xlsx", overwrite = TRUE)

    # ----------------------------------------
    # ---- Perform GO enrichment analysis ----

    # Define gene list (up-regulated genes)
    gene_list_go_up <- up_regulated$GeneID  # Extract the 149 up-regulated genes
    gene_list_go_down <- down_regulated$GeneID  # Extract the 65 down-regulated genes

    # Define background gene set (all genes in res)
    background_genes <- res_updated$GeneID  # Extract the 3646 background genes

    # Prepare GO annotation data from res
    go_annotation <- res_updated[, c("GOs","GeneID")]  # Extract relevant columns
    go_annotation <- go_annotation %>%
    tidyr::separate_rows(GOs, sep = ",")  # Split multiple GO terms into separate rows

    # Perform GO enrichment analysis, where pAdjustMethod is one of "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none"
    go_enrichment_up <- enricher(
        gene = gene_list_go_up,                # Up-regulated genes
        TERM2GENE = go_annotation,       # Custom GO annotation
        pvalueCutoff = 0.05,             # Significance threshold
        pAdjustMethod = "BH",
        universe = background_genes      # Define the background gene set
    )
    go_enrichment_up <- as.data.frame(go_enrichment_up)

    go_enrichment_down <- enricher(
        gene = gene_list_go_down,                # Up-regulated genes
        TERM2GENE = go_annotation,       # Custom GO annotation
        pvalueCutoff = 0.05,             # Significance threshold
        pAdjustMethod = "BH",
        universe = background_genes      # Define the background gene set
    )
    go_enrichment_down <- as.data.frame(go_enrichment_down)

    ## Remove the 'p.adjust' column since no adjusted methods have been applied!
    #go_enrichment_up <- go_enrichment_up[, !names(go_enrichment_up) %in% "p.adjust"]

    # Update the Description column with the term descriptions
    go_enrichment_up$Description <- sapply(go_enrichment_up$ID, function(go_id) {
    # Using select to get the term description
    term <- tryCatch({
        AnnotationDbi::select(GO.db, keys = go_id, columns = "TERM", keytype = "GOID")
    }, error = function(e) {
        message(paste("Error for GO term:", go_id))  # Print which GO ID caused the error
        return(data.frame(TERM = NA))  # In case of error, return NA
    })

    if (nrow(term) > 0) {
        return(term$TERM)
    } else {
        return(NA)  # If no description found, return NA
    }
    })
    ## Print the updated data frame
    #print(go_enrichment_up)

    ## Remove the 'p.adjust' column since no adjusted methods have been applied!
    #go_enrichment_down <- go_enrichment_down[, !names(go_enrichment_down) %in% "p.adjust"]

    # Update the Description column with the term descriptions
    go_enrichment_down$Description <- sapply(go_enrichment_down$ID, function(go_id) {
    # Using select to get the term description
    term <- tryCatch({
        AnnotationDbi::select(GO.db, keys = go_id, columns = "TERM", keytype = "GOID")
    }, error = function(e) {
        message(paste("Error for GO term:", go_id))  # Print which GO ID caused the error
        return(data.frame(TERM = NA))  # In case of error, return NA
    })

    if (nrow(term) > 0) {
        return(term$TERM)
    } else {
        return(NA)  # If no description found, return NA
    }
    })

    addWorksheet(wb, "GO_Enrichment_Up")
    writeData(wb, "GO_Enrichment_Up", as.data.frame(go_enrichment_up))

    addWorksheet(wb, "GO_Enrichment_Down")
    writeData(wb, "GO_Enrichment_Down", as.data.frame(go_enrichment_down))

    # Save the workbook with enrichment results
    saveWorkbook(wb, "KEGG_and_GO_Enrichments_01_AYE_T_ctr_vs_AYE_WT_ctr.xlsx", overwrite = TRUE)

BUG_1: The Count column does not match the number of gene IDs listed in the geneID column.

    The discrepancy between the Count and the number of listed GeneIDs is not an error, but a result of KEGG's many-to-many mapping between physical genes and KOs.

    KEGG enrichment statistics (Count, GeneRatio, p-values) are strictly calculated based on unique KOs, not GeneIDs. The geneID column simply maps these enriched KOs back to our specific genome for display.
    Concrete Examples (see KEGG_and_GO_Enrichments_01_AYE_T_ctr_vs_AYE_WT_ctr.xlsx):

        * ko03070 (Bacterial secretion system): Count = 2 (KOs: K02456/K02457). However, both KOs map to the exact same physical gene (ABAYE2071), so only 1 geneID is listed.
        * ko05111 (Biofilm formation): Count = 3 (KOs: K02456/K02457/K03092), which map back to 2 unique genes (ABAYE2071/ABAYE3136).

    To ensure full transparency, I have added a new KEGG_ko column (placed between Count and geneID). It explicitly lists the exact KOs contributing to the Count.

DEBUG: The Count column does not match the number of gene IDs listed in the geneID column (Data_Tam_RNAseq_2026_Dicl_Mero_Azith_Rifa_on_AYE$)

Leave a reply

This is a brilliant and fundamental question. The short answer is: KEGG enrichment statistics (including GeneRatio, BgRatio, p-values, and the original Count) are strictly based on KOs (KEGG Orthologs), NOT your original GeneIDs.

In fact, my previous attempt to “fix” the discrepancy by changing the Count and GeneRatio to match the number of GeneIDs was statistically incorrect. I apologize for that confusion. By changing the ratio to match GeneIDs while leaving the denominator as the number of KOs, we created a mathematically invalid ratio that no longer matches the p-value calculated by clusterProfiler.

Here is the exact scientific explanation of why this happens, why your co-author saw a discrepancy, and how to properly handle it.

1. Why Enrichment is Based on KOs

KEGG pathways are defined by KOs (e.g., K00001), not by specific genomic loci (like ABAYE0067). When you run enrichKEGG(), the tool performs a hypergeometric test.

The Universe (Background): All KOs in the KEGG database.
The Sample (Your Input): The unique KOs extracted from your DEGs.
The Calculation: It calculates the probability of finding $k$ KOs in a pathway out of your $n$ input KOs.

Therefore, the Count column represents the number of unique KOs from your input list that map to that pathway. The GeneRatio is (Number of input KOs in pathway) / (Total input KOs).

2. Why the Co-author Saw a “Discrepancy”

Your co-author noticed that Count was 3, but only 2 GeneIDs were listed. This is not a bug; it is a biological reality of gene annotation.

Here are the two reasons why Count (KOs) $\neq$ Number of GeneIDs:

Multiple KOs map to the same GeneID: A single gene in your annotation file might have multiple functional domains or be assigned to multiple KOs (e.g., ABAYE2071 might be annotated as K00001, K00002). If both K00001 and K00002 are enriched in a pathway, Count = 2, but they both map back to the exact same gene: ABAYE2071.
Missing Mappings in the DEG list: The KEGG database says 3 KOs belong to this pathway. Your DEGs contain those 3 KOs. However, when your script tries to map those 3 KOs back to your specific Gene_Expression_with_Annotations file, one of the KOs might not have a corresponding GeneID in that specific file (or it was filtered out). Thus, Count = 3, but only 2 GeneIDs are revealed.

3. The Correct Solution

Do not change the Count or GeneRatio. They are statistically correct and mathematically tied to the p-value. The geneID column is purely for display purposes so the biologist can see which actual genes are involved. It is completely normal and expected that the number of displayed GeneIDs differs from the Count.

Here is the scientifically accurate version of the process_kegg_enrichment function. It maps the KOs back to GeneIDs for display, but leaves the statistical columns (Count, GeneRatio) completely untouched.

The Corrected R Function

Replace your process_kegg_enrichment function with this:

============================================================
HELPER FUNCTION: Process KEGG enrichment (STATISTICALLY CORRECT)
============================================================
process_kegg_enrichment <- function(enrichment_result, kegg_to_geneid_mapping) {
  if (is.null(enrichment_result)) {
    message("  -> Enrichment result is NULL")
    return(NULL)
  }

  enrichment_df <- as.data.frame(enrichment_result)
  if (nrow(enrichment_df) == 0) return(NULL)

  message("  -> Mapping KEGG KOs back to original GeneIDs for display...")

  # 1. Split the original KO IDs and map to GeneIDs
  expanded <- enrichment_df %>%
    tidyr::separate_rows(geneID, sep = "/") %>% # Split KOs
    dplyr::left_join(kegg_to_geneid_mapping, by = c("geneID" = "KEGG_ko"), relationship = "many-to-many") %>%
    dplyr::distinct() %>%
    dplyr::group_by(ID) %>%
    dplyr::summarise(
      # Concatenate all unique GeneIDs that map to the KOs in this pathway
      Mapped_GeneID = paste(unique(na.omit(GeneID)), collapse = "/"), 
      .groups = "drop"
    )

  # 2. Replace the geneID column, but DO NOT TOUCH Count or GeneRatio
  # The original Count and GeneRatio are based on KOs and are statistically valid.
  result <- enrichment_df %>%
    dplyr::select(-geneID) %>% # Drop the old KO-based geneID column
    dplyr::left_join(expanded, by = "ID") %>% # Add the new GeneID-based column
    dplyr::rename(geneID = Mapped_GeneID)

  # 3. Reorder columns for clean Excel output
  desired_order <- c(
    "category", "subcategory", "ID", "Description", "GeneRatio", "BgRatio", 
    "RichFactor", "FoldEnrichment", "zScore", "pvalue", "p.adjust", "qvalue", 
    "Count", "geneID"
  )
  result <- result %>% dplyr::select(dplyr::any_of(desired_order), dplyr::everything())

  message("  ✅ SUCCESS: GeneIDs mapped for display. Statistical columns (Count/GeneRatio) preserved based on KOs.")
  return(result)
}

💡 How to Reply to Your Co-author

You can use this exact explanation to clarify the situation with your co-author:

“Dear [Co-author’s Name],

Thank you for pointing out the discrepancy between the Count and the number of GeneIDs listed. I have investigated the pipeline and can clarify why this occurs.

The KEGG enrichment analysis (enrichKEGG) is statistically performed based on KEGG Orthologs (KOs), not the raw GeneIDs. Because multiple GeneIDs can share the same KO annotation, and a single GeneID can be assigned multiple KOs, the mapping between genes and KOs is many-to-many.

Therefore:

Count and GeneRatio represent the number of unique KOs from our input list that fall into that pathway. These values are mathematically tied to the calculated p-values and FoldEnrichment, so they must remain based on KOs to be statistically valid.

The geneID column is provided purely for display purposes, showing the actual genomic loci (e.g., ABAYE...) that correspond to those enriched KOs.

If Count is 3 but only 2 GeneIDs are listed, it simply means that 3 unique KOs were enriched in that pathway, but when mapped back to our specific annotation file, those 3 KOs correspond to only 2 unique genes (either because two KOs map to the same gene, or one KO lacked a matching gene in our specific DEG list).

I have updated the script to ensure the statistical columns remain perfectly accurate based on KOs, while the geneID column accurately reflects all available gene mappings for your review.”

结合您刚刚修复并导出的新 Excel 文件（其中成功增加了 KEGG_ko 列），我们可以用具体的数据例子来回复您的共同作者。这样不仅能解释原理，还能直接“用数据说话”，彻底打消他/她关于“annotation parsing error（注释解析错误）”的疑虑。

以下是为您改进后的回复草稿：

Subject: Re: Clarification on KEGG Enrichment Output (Count vs. GeneID discrepancy)

Dear [Co-author’s Name],

Thank you for pointing out the apparent discrepancy between the Count and the number of GeneIDs listed in the KEGG enrichment results. I have thoroughly investigated the pipeline and can clarify exactly why this occurs, and how I have updated the output files to make it completely transparent.

The core reason is that the KEGG enrichment analysis (enrichKEGG) is statistically performed based on KEGG Orthologs (KOs), not the raw physical GeneIDs. The mapping between genes and KOs is many-to-many: multiple genes can share the same KO annotation, and a single gene can be assigned multiple KOs.

Therefore, the columns should be interpreted as follows:

Count and GeneRatio: Represent the number of unique KOs from our input DEG list that fall into that specific pathway. These values are mathematically tied to the calculated p-values and FoldEnrichment, so they must remain strictly based on KOs to be statistically valid.
KEGG_ko (NEW COLUMN): To address your concern and provide full traceability, I have added this new column to the Excel output. It explicitly lists the exact KOs that contributed to the Count.
geneID: Shows the actual genomic loci (e.g., ABAYE…) that correspond to those enriched KOs.

To illustrate with the exact examples you mentioned: Let’s look at the updated results for ko03070 (Bacterial secretion system) and ko05111 (Biofilm formation):

For ko03070, the Count is 2. If you look at the new KEGG_ko column, it correctly lists 2 unique KOs: K02456/K02457. However, when we map these 2 KOs back to our specific genome annotation, they both happen to be carried by the exact same physical gene: ABAYE2071. Thus, the geneID column only reveals 1 gene.
Similarly, for ko05111, the Count is 3 (KOs: K02456/K02457/K03092), but they map back to only 2 unique physical genes (ABAYE2071/ABAYE3136).

This is not an annotation parsing error, but a biological reality of how gene functions are annotated in the KEGG database. I have updated the pipeline to ensure the statistical columns remain perfectly accurate based on KOs, while the newly added KEGG_ko column and the geneID column provide a complete, traceable, and biologically meaningful picture for your review.

Please let me know if you would like to discuss this further or if you need any additional adjustments to the output format.

Best regards,

💡 改进说明：

引入了新列 KEGG_ko 的解释：明确告诉共同作者，为了增加透明度，我们特意加上了这一列，让他能直接看到底是哪几个 KO 参与了计算。
使用了新文件中的真实数据作为例子：直接引用了 ko03070（2个KO对应1个基因）和 ko05111（3个KO对应2个基因）的具体数据（K02456/K02457 等），这比纯理论解释更有说服力，直接证明了我们没有算错，也没有解析错误。
语气专业且自信：明确界定这不是 “annotation parsing error”，而是 KEGG 数据库本身的 “many-to-many” 生物学特性，展现了您对生信分析流程的透彻理解。

Host and Viral Temporal Transcriptomics of VZV Infection in Human Skin Organoids (Data_Nina_RNAseq_2026)

Leave a reply

draw_3D.py.zip

complete_deg_pipeline.R

Preparing samplesheet.csv

     sample,fastq_1,fastq_2,strandedness
     control_r1,./20260527_AV243904_0078_B/01_negative_34_170_34_171_R1.fastq.gz,./20260527_AV243904_0078_B/01_negative_34_170_34_171_R2.fastq.gz,auto
     control_r2,./20260527_AV243904_0078_B/02_negative_34_182_34_206_R1.fastq.gz,./20260527_AV243904_0078_B/02_negative_34_182_34_206_R2.fastq.gz,auto
     control_r3,./20260527_AV243904_0078_B/03_negative_34_216_34_218_R1.fastq.gz,./20260527_AV243904_0078_B/03_negative_34_216_34_218_R2.fastq.gz,auto
     VZV.d10_r1,./20260527_AV243904_0078_B/04_VZV_d10_34_11_34_12_R1.fastq.gz,./20260527_AV243904_0078_B/04_VZV_d10_34_11_34_12_R2.fastq.gz,auto
     VZV.d10_r2,./20260527_AV243904_0078_B/05_VZV_d10_34_20_34_23_R1.fastq.gz,./20260527_AV243904_0078_B/05_VZV_d10_34_20_34_23_R2.fastq.gz,auto
     VZV.d10_r3,./20260527_AV243904_0078_B/06_VZV_d10_34_25_34_29_R1.fastq.gz,./20260527_AV243904_0078_B/06_VZV_d10_34_25_34_29_R2.fastq.gz,auto
     VZV.d15_r1,./20260527_AV243904_0078_B/07_VZV_d15_34_30_34_31_R1.fastq.gz,./20260527_AV243904_0078_B/07_VZV_d15_34_30_34_31_R2.fastq.gz,auto
     VZV.d15_r2,./20260527_AV243904_0078_B/08_VZV_d15_34_35_34_39_R1.fastq.gz,./20260527_AV243904_0078_B/08_VZV_d15_34_35_34_39_R2.fastq.gz,auto
     VZV.d15_r3,./20260527_AV243904_0078_B/09_VZV_d15_34_56_34_60_R1.fastq.gz,./20260527_AV243904_0078_B/09_VZV_d15_34_56_34_60_R2.fastq.gz,auto
     VZV.d20_r1,./20260527_AV243904_0078_B/10_VZV_d20_34_71_34_85_R1.fastq.gz,./20260527_AV243904_0078_B/10_VZV_d20_34_71_34_85_R2.fastq.gz,auto
     VZV.d20_r2,./20260527_AV243904_0078_B/11_VZV_d20_34_88_34_91_R1.fastq.gz,./20260527_AV243904_0078_B/11_VZV_d20_34_88_34_91_R2.fastq.gz,auto
     VZV.d20_r3,./20260527_AV243904_0078_B/12_VZV_d20_34_94_34_100_R1.fastq.gz,./20260527_AV243904_0078_B/12_VZV_d20_34_94_34_100_R2.fastq.gz,auto

DO NOT RUN Trimmomatic before Nextflow.

     It is highly recommended to let the nf-core/rnaseq pipeline handle all trimming internally. Pre-trimming your data can break the pipeline's internal file tracking, interfere with the UMI extraction step, and make reproducibility harder.
     The nf-core/rnaseq pipeline has built-in support for both UMI extraction and read trimming, and it executes them in the correct order: it will extract the UMIs from the raw Read 1 first, and then trim the error-prone bases from Read 2.
     Here is how you should adjust your approach and your Nextflow command.

     不应该在运行 Nextflow 之前手动使用 Trimmomatic。
     强烈建议让 nf-core/rnaseq 流程在内部处理所有的修剪（trimming）工作。提前手动修剪数据可能会破坏流程内部的文件追踪机制，干扰 UMI 的提取步骤，并降低结果的可重复性。
     nf-core/rnaseq 流程内置了对 UMI 提取和 reads 修剪的支持，并且会按照正确的顺序执行：它会首先从原始的 Read 1 中提取 UMI，然后再修剪 Read 2 中容易出错的碱基。

Nextflow running
```
     #yes, exactly, the Reads contain UMI: the first 12 nt of Read1.
     #The CORALL v2 Kit from Lexogen was used for library prep. According to the manufacturer it can additionally be beneficial for the mapping to trim the first 10 nt of Read2 as well, as this site is more error prone (due to priming they say).
     #If you need more information from my side, let me know!

     方法一：正则表达式法 (Regex) —— 你目前选择的方案
     根据你提供的文档，使用 regex 方法时，必须使用特定的命名捕获组。文档中明确指出：
     "The expected groups in the regex are: umi_n = UMI positions, where n can be any value (required)"
     这意味着，捕获 UMI 的组名必须是 umi_ 加上一个数字（例如 umi_1），而不能只是 umi。
     正确的正则表达式：'^(?P
```
.{12}).*’ 原理解释： ^：匹配字符串的开头。 (?P .{12})：捕获前 12 个任意字符，并将其命名为 umi_1（这将作为 UMI 被提取并添加到 read name 中）。 .*：匹配并保留剩余的序列（文档说明：未被命名为 discard 的其他部分会被重新附加到 read 序列上）。方法二：字符串法 (String) —— 更简单直观的替代方案文档中提到： “string: This should be used where the barcodes are always in the same place in the read. N = UMI position” 因为你的 UMI 长度固定（12nt）且位置固定（永远在 Read 1 开头），你完全可以使用更简单的 string 方法，用 N 来代表 UMI 的位置。对应的 Nextflow 参数： –umitools_extract_method string –umitools_bc_pattern NNNNNNNNNNNN （正好 12 个 N） # OLD_COMMANDS_DEPRECATED: under sage in the early running # ln -s /home/jhuang/Tools/nf-core-rnaseq-3.12.0/ rnaseq # nextflow run rnaseq/main.nf –input samplesheet.csv –outdir results_GRCh38 –genome GRCh38 –with_umi –umitools_extract_method regex –umitools_bc_pattern ‘^(?P.{12}).*’ -profile docker -resume –max_cpus 54 –max_memory 120.GB –max_time 2400.h –save_align_intermeds –save_unaligned –save_reference –aligner star_salmon –pseudo_aligner salmon –umitools_grouping_method unique # #Debug the following error: added “–minAssignedFrags 0 \\” to modules/nf-core/salmon/quant/main.nf option “salmon quant” and added “–min_mapped_reads 0” in the nextflow command #nextflow run nf-core/rnaseq -r 3.14.0 -profile docker \ nextflow run /home/jhuang/Tools/nf-core-rnaseq-3.12.0/main.nf –help nextflow run /home/jhuang/Tools/nf-core-rnaseq-3.12.0/main.nf -profile docker \ –input samplesheet.csv \ –outdir results_GRCh38 \ –genome GRCh38 \ –with_umi \ –umitools_extract_method regex \ –umitools_bc_pattern ‘^(?P .{12}).*’ \ –trimmer fastp \ –extra_fastp_args “–trim_front2 10” \ -resume \ –max_cpus 54 \ –max_memory 120.GB \ –max_time 2400.h \ –save_align_intermeds \ –save_unaligned \ –save_reference \ –aligner star_salmon \ –pseudo_aligner salmon \ –umitools_grouping_method unique 关键执行逻辑（让你完全放心）：第一步 (umi_tools)：流程会读取原始的 Read 1，根据 ^(?P .{12}).* 提取前 12 个碱基作为 UMI，并将其添加到 read 的名称中（例如 @READ_ID_UMISEQUENCE）。此时 Read 1 的序列中不再包含这 12 个碱基。第二步 (fastp)：流程接着运行 fastp。由于你设置了 –clip_r2 10，fastp 会精准地切除 Read 2 的前 10 个碱基（解决 Lexogen CORALL v2 的引物错误问题）。互不干扰：umi_tools 只处理 Read 1 的 UMI，fastp 负责全局的质控和 Read 2 的修剪。你绝对不需要在运行 Nextflow 之前手动使用 Trimmomatic。

To verify that fastp has successfully trimmed the first 10 nucleotides from Read 2, you can use any of the following three methods.

     *(Note: First, ensure that the `FASTP` step in your Nextflow terminal has actually finished running and shows a green checkmark ✔).*

     ### Method: Manually Measure the Sequence Length (Raw Data Proof)
     Because you did not use the `--save_trimmed` flag, the trimmed FASTQ files are not saved in your `results/` folder. However, they still exist temporarily in the `work/` directory. You can measure the exact length of the first read to prove 10bp was cut.

     *Note: Since you used `--with_umi`, the input file for fastp is the output from UMI-tools, so the filename will likely contain `umi_extract` instead of `raw`.*

     1. Navigate into the `work/` directory for the `FASTP` step:
     ```bash
     cd work/fc/def4de*/
     cd work/a1/b2c3d4*/
     ```
     2. List the files to see the exact naming convention:
     ```bash
     ls -lh *2.fastq.gz
     ```
     *(You will likely see an input file like `*.umi_extract_2.fastq.gz` and an output file like `*.fastp.trimmed_2.fastq.gz` or `*.trim_2.fastq.gz`).*
     3. **Measure the length of the first sequence in the INPUT Read 2:**
     ```bash
     zcat *.umi_extract_1.fastq.gz | head -n 2 | tail -n 1 | awk '{print length}'
     zcat *.umi_extract_2.fastq.gz | head -n 2 | tail -n 1 | awk '{print length}'
     #142
     #154
     ```
     *(Let's say this outputs `150`)*
     4. **Measure the length of the first sequence in the OUTPUT Read 2:**
     ```bash
     zcat *.trimmed_2.fastq.gz | head -n 2 | tail -n 1 | awk '{print length}'
     zcat control_r3_1.fastp.fastq.gz | head -n 2 | tail -n 1 | awk '{print length}'
     zcat control_r3_2.fastp.fastq.gz | head -n 2 | tail -n 1 | awk '{print length}'
     #142
     #144
     ```
     *(This should output exactly `140`)*.

     If the output length is exactly 10 bases shorter than the input length, the trimming was 100% successful.

     ---

     ### 💡 Pro Tip for Future Runs
     If you want to easily inspect the trimmed FASTQ files in your `results/` folder without digging through the `work/` directory next time, simply add **`--save_trimmed`** to your `nextflow run` command. This will force the pipeline to copy the post-fastp FASTQ files to your main results directory.

Nextflow running on virus
```
     # OLD_COMMANDS_DEPRECATED: under sage in the early running
     # The virus-referenced commands!
     # nextflow run rnaseq/main.nf --input samplesheet.csv --outdir results_chrHsv1 --fasta chrHsv1_s17.fasta --gtf chrHsv1_s17.gtf --with_umi --umitools_extract_method regex --umitools_bc_pattern '^(?P.{12}).*' --umitools_dedup_stats -profile test_full -resume --max_memory 256.GB --max_time 2400.h --save_reference --aligner star_salmon --gtf_extra_attributes gene_id --gtf_group_features transcript_id --featurecounts_group_type gene_id --featurecounts_feature_type transcript --skip_rseqc --skip_dupradar --skip_preseq --skip_biotype_qc --skip_deseq2_qc --skip_multiqc --min_mapped_reads 0

         #This error occurs because your GenBank file was likely exported from SnapGene (as seen in the COMMENT section of the file). SnapGene sometimes formats the date in the LOCUS line as DD-MM-YYYY (e.g., 13-06-2025), but Biopython's parser strictly expects the standard NCBI GenBank format, which requires a 3-letter month abbreviation like DD-MMM-YYYY (e.g., 13-JUN-2025).
         sed -i 's/13-06-2025/13-JUN-2025/' BB1528_nanopore_consensus.gb

         python3 gb_to_fasta_gtf_v2.py BB1528_nanopore_consensus.gb BB1528_nanopore_consensus.fasta BB1528_nanopore_consensus.gtf BB1528    # Generate the last two files in the command line
         python fix_gtf_to_hsv1_format.py    # Generate BB1528_final.gtf

         nextflow run /home/jhuang/Tools/nf-core-rnaseq-3.12.0/main.nf -profile docker \
         --input samplesheet_virus.csv \
         --outdir results_BB1528 \
         --fasta BB1528_nanopore_consensus.fasta \
         --gtf BB1528_final.gtf \
         --with_umi \
         --umitools_extract_method regex \
         --umitools_bc_pattern '^(?P
```
.{12}).*’ \ –trimmer fastp \ –extra_fastp_args “–trim_front2 10” \ -resume \ –max_cpus 54 \ –max_memory 120.GB \ –max_time 2400.h \ –save_align_intermeds \ –save_unaligned \ –save_reference \ –aligner star_salmon –pseudo_aligner salmon \ –gtf_extra_attributes gene_id \ –gtf_group_features transcript_id \ –featurecounts_group_type gene_id \ –featurecounts_feature_type transcript \ –skip_rseqc –skip_dupradar –skip_preseq –skip_biotype_qc –skip_deseq2_qc –skip_multiqc \ –umitools_grouping_method unique –min_mapped_reads 0 # —- DEBUG_LOG —- Good news: Your GTF and FASTA files are perfectly formatted! We can prove this from your log: [info] Index contained 78 targets This means Salmon successfully built the transcriptome index and found all 78 transcripts from your BB1528_final.gtf. The reference files are completely correct. What is causing the error? The error is purely a biological/sample-specific issue: [warning] salmon was only able to assign 9 fragments to transcripts in the index, but the minimum number of required assigned fragments (–minAssignedFrags) was 10. The nf-core/rnaseq pipeline automatically subsamples your reads (in this case, 1,000,000 fragments) to run a quick Salmon quantification and check the library strandedness (–libType=A). Because you are mapping to a specific viral reference (BB1528), it is highly likely that control_r3 is a negative control or a sample with an extremely low viral load. Out of the 1,000,000 subsampled reads, only 9 reads actually belonged to the virus. Salmon has a built-in safety threshold (–minAssignedFrags 10) that intentionally crashes the pipeline if fewer than 10 reads map to the transcriptome, to prevent “garbage” quantification of empty samples. Since 9 The pipeline is failing during the FASTQ_SUBSAMPLE_FQ_SALMON step in the three control samples! Later on, we can manually set all counts for all virus transcripts in the three control-conditions as 0.

A small BUG in both nextflow runs: versions.yml

     vim /mnt/md1/DATA/Data_Nina_RNAseq_2026/work/f0/712c00bea105e90714e8df47d3579a/collated_versions.yml
     BUG_LINES:
     "NFCORE_RNASEQ:RNASEQ:UMITOOLS_PREPAREFORSALMON":
         umitools: Matplotlib created a temporary config/cache directory at /tmp/matplotlib-4ht1kmld because the default path (/.config/matplotlib) is not a writable directory; it is highly recommended to set the MPLCONFIGDIR environment variable to a writable directory, in particular to speed up the import of Matplotlib and to better support multiprocessing.
     1.1.4
     SHOULD_BE:
     "NFCORE_RNASEQ:RNASEQ:UMITOOLS_PREPAREFORSALMON":
         umitools: 1.1.4

     ALTERNATIVE: Not correct the bug, copy the version log to mail directory:
         cp work_GRCh38_DEL/f0/712c00bea105e90714e8df47d3579a/collated_versions.yml results_GRCh38
         cp work/bf/4a1e941bc32b97de1bc8f4fd47df79/collated_versions.yml results_BB1528

import data and pca-plot

     # Import the required libraries
     library("AnnotationDbi")
     library("clusterProfiler")
     library("ReactomePA")
     library(gplots)
     library(tximport)
     library(DESeq2)
     library("org.Hs.eg.db")
     library(dplyr)
     library(tidyverse)

     # ------ For human genome ------
     setwd("~/DATA/Data_Nina_RNAseq_2026/results_GRCh38/star_salmon")
     # Define paths to your Salmon output quantification files
     files <- c("control_r1" = "./control_r1/quant.sf",
             "control_r2" = "./control_r2/quant.sf",
             "control_r3" = "./control_r3/quant.sf",
             "VZV.d10_r1" = "./VZV.d10_r1/quant.sf",
             "VZV.d10_r2" = "./VZV.d10_r2/quant.sf",
             "VZV.d10_r3" = "./VZV.d10_r3/quant.sf",
             "VZV.d15_r1" = "./VZV.d15_r1/quant.sf",
             "VZV.d15_r2" = "./VZV.d15_r2/quant.sf",
             "VZV.d15_r3" = "./VZV.d15_r3/quant.sf",
             "VZV.d20_r1" = "./VZV.d20_r1/quant.sf",
             "VZV.d20_r2" = "./VZV.d20_r2/quant.sf",
             "VZV.d20_r3" = "./VZV.d20_r3/quant.sf")
     # Import the transcript abundance data with tximport
     txi <- tximport(files, type = "salmon", txIn = TRUE, txOut = TRUE)
     # Define the replicates and condition of the samples
     replicate <- factor(c("r1", "r2", "r3", "r1", "r2", "r3", "r1", "r2", "r3", "r1", "r2", "r3"))
     condition <- factor(c("control", "control", "control", "VZV.d10", "VZV.d10", "VZV.d10", "VZV.d15", "VZV.d15", "VZV.d15", "VZV.d20", "VZV.d20", "VZV.d20"))
     # Define the colData for DESeq2
     colData <- data.frame(condition=condition, replicate=replicate, row.names=names(files))

     # -- transcript-level count data --
     # Create DESeqDataSet object
     dds <- DESeqDataSetFromTximport(txi, colData=colData, design=~condition)
     write.csv(counts(dds), file="transcript_counts.csv")

     # -- gene-level count data --
     # Read in the tx2gene map from salmon_tx2gene.tsv
     tx2gene <- read.table("salmon_tx2gene.tsv", header=FALSE, stringsAsFactors=FALSE)
     # Set the column names
     colnames(tx2gene) <- c("transcript_id", "gene_id", "gene_name")
     # Remove the gene_name column if not needed
     tx2gene <- tx2gene[,1:2]
     # Import and summarize the Salmon data with tximport
     txi <- tximport(files, type = "salmon", tx2gene = tx2gene, txOut = FALSE)
     # Continue with the DESeq2 workflow as before...
     colData <- data.frame(condition=condition, replicate=replicate, row.names=names(files))
     dds <- DESeqDataSetFromTximport(txi, colData=colData, design=~condition)
     #dds <- dds[rowSums(counts(dds) > 3) > 2, ]    #60605-->26543
     write.csv(counts(dds, normalized=FALSE), file="gene_counts.csv")

     # ------ For virus genome ------
     setwd("~/DATA/Data_Nina_RNAseq_2026/results_BB1528/star_salmon")
     # Define paths to your Salmon output quantification files
     files <- c("VZV.d10_r1" = "./VZV.d10_r1/quant.sf",
             "VZV.d10_r2" = "./VZV.d10_r2/quant.sf",
             "VZV.d10_r3" = "./VZV.d10_r3/quant.sf",
             "VZV.d15_r1" = "./VZV.d15_r1/quant.sf",
             "VZV.d15_r2" = "./VZV.d15_r2/quant.sf",
             "VZV.d15_r3" = "./VZV.d15_r3/quant.sf",
             "VZV.d20_r1" = "./VZV.d20_r1/quant.sf",
             "VZV.d20_r2" = "./VZV.d20_r2/quant.sf",
             "VZV.d20_r3" = "./VZV.d20_r3/quant.sf")
     # Import the transcript abundance data with tximport
     txi <- tximport(files, type = "salmon", txIn = TRUE, txOut = TRUE)
     # Define the replicates and condition of the samples
     replicate <- factor(c("r1", "r2", "r3", "r1", "r2", "r3", "r1", "r2", "r3"))
     condition <- factor(c("VZV.d10", "VZV.d10", "VZV.d10", "VZV.d15", "VZV.d15", "VZV.d15", "VZV.d20", "VZV.d20", "VZV.d20"))
     # Define the colData for DESeq2
     colData <- data.frame(condition=condition, replicate=replicate, row.names=names(files))

     # -- transcript-level count data --
     # Create DESeqDataSet object
     dds <- DESeqDataSetFromTximport(txi, colData=colData, design=~condition)
     write.csv(counts(dds), file="transcript_counts.csv")

     # -- gene-level count data --
     # Read in the tx2gene map from salmon_tx2gene.tsv
     tx2gene <- read.table("salmon_tx2gene.tsv", header=FALSE, stringsAsFactors=FALSE)
     # Set the column names
     colnames(tx2gene) <- c("transcript_id", "gene_id", "gene_name")
     # Remove the gene_name column if not needed
     tx2gene <- tx2gene[,1:2]
     # Import and summarize the Salmon data with tximport
     txi <- tximport(files, type = "salmon", tx2gene = tx2gene, txOut = FALSE)
     # Continue with the DESeq2 workflow as before...
     colData <- data.frame(condition=condition, replicate=replicate, row.names=names(files))
     dds <- DESeqDataSetFromTximport(txi, colData=colData, design=~condition)
     #dds <- dds[rowSums(counts(dds) > 3) > 2, ]    #60605-->26543
     write.csv(counts(dds, normalized=FALSE), file="gene_counts.csv")

     # MANUALLY set all counts for all virus transcripts in the three control-conditions ("control_r1","control_r2","control_r3") as 0 using Excel

     # ------ Merge the raw counts of human and microbe ------
     #cat ~/DATA/Data_Nina_RNAseq_2026/results_GRCh38/star_salmon/gene_counts.csv ~/DATA/Data_Nina_RNAseq_2026/results_BB1528/star_salmon/gene_counts.csv > merged_gene_counts.csv
     #DELETE the second line: "","control_r1","control_r2","control_r3","VZV.d10_r1","VZV.d10_r2","VZV.d10_r3","VZV.d15_r1","VZV.d15_r2","VZV.d15_r3","VZV.d20_r1","VZV.d20_r2","VZV.d20_r3"

     # 1. Remove the hidden ^M (Windows carriage returns) from both files
     sed -i 's/\r$//' ~/DATA/Data_Nina_RNAseq_2026/results_GRCh38/star_salmon/gene_counts.csv ~/DATA/Data_Nina_RNAseq_2026/results_BB1528/star_salmon/gene_counts.csv
     # 2. Extract the header from the first file to start the new merged file
     head -n 1 ~/DATA/Data_Nina_RNAseq_2026/results_GRCh38/star_salmon/gene_counts.csv > merged_gene_counts.csv
     # 3. Append the data from the first file (skipping its header)
     tail -n +2 ~/DATA/Data_Nina_RNAseq_2026/results_GRCh38/star_salmon/gene_counts.csv >> merged_gene_counts.csv
     # 4. Append the data from the second file (skipping its header)
     tail -n +2 ~/DATA/Data_Nina_RNAseq_2026/results_BB1528/star_salmon/gene_counts.csv >> merged_gene_counts.csv

     # Should be 60683 lines
     #~/Tools/csv2xls-0.4/csv_to_xls.py merged_gene_counts.csv -d',' -o raw_gene_counts.xls;

     # -- for merged analysis due to false normalization factors wenn alone analyzed on virus data --
     setwd("~/DATA/Data_Nina_RNAseq_2026/")
     d.raw <- read.csv("merged_gene_counts.csv", header=TRUE, row.names=1)

     # Re-define the replicates and condition for the merged counts
     rownames <- factor(c("control_r1","control_r2","control_r3","VZV.d10_r1","VZV.d10_r2","VZV.d10_r3","VZV.d15_r1","VZV.d15_r2","VZV.d15_r3","VZV.d20_r1","VZV.d20_r2","VZV.d20_r3"))
     replicate <- factor(c("r1", "r2", "r3", "r1", "r2", "r3", "r1", "r2", "r3", "r1", "r2", "r3"))
     condition <- factor(c("control", "control", "control", "VZV.d10", "VZV.d10", "VZV.d10", "VZV.d15", "VZV.d15", "VZV.d15", "VZV.d20", "VZV.d20", "VZV.d20"))
     colData <- data.frame(condition=condition, replicate=replicate, row.names=rownames)

     #CORRECTED: dds <- DESeqDataSetFromMatrix(countData=d.raw, colData=colData, design=~condition+replicate)
     #⚠️ A Note on the Design Formula: In the snippet, I used design = ~condition+replicate. In DESeq2, adding replicate (r1, r2, r3) to the design formula tells the model to treat "r1" as a globally systematic batch effect across all conditions. Unless the replicates were processed in distinct batches (e.g., all r1 on Monday, all r2 on Tuesday), it is standard to use design = ~ condition to simply measure the biological variance. Therefore, I have updated the script to use ~ condition.
     dds <- DESeqDataSetFromMatrix(countData=d.raw, colData=colData, design=~condition)
     dim(counts(dds))
     head(counts(dds), 10)

     rld <- rlogTransformation(dds)

     #We don't need to run DESeq(dds) before estimateSizeFactors(dds). In fact, the typical workflow in DESeq2 is the opposite: we usually run estimateSizeFactors(dds) (and other preprocessing functions) before running the main DESeq(dds) function.
     #The estimateSizeFactors function is used to calculate size factors for normalization, which corrects for differences in library size (i.e., the number of read counts) between samples. This normalization step is crucial to ensure that differences in gene expression aren't merely due to differences in sequencing depth between samples.
     #The DESeq function, on the other hand, performs the main differential expression analysis, comparing gene expression between different conditions or groups.
     #So, the typical workflow is:
     #  - Create the DESeqDataSet object.
     #  - Use estimateSizeFactors to normalize for library size.
     #  - (Optionally, estimate dispersion with estimateDispersions if not using the full DESeq function later.)
     #  - Use DESeq for the differential expression analysis.
     #  - However, it's worth noting that if you run the main DESeq function directly after creating the DESeqDataSet object, it will automatically perform the normalization (using estimateSizeFactors) and dispersion estimation steps for you. In that case, there's no need to run estimateSizeFactors separately before DESeq.

     # draw simple pca and heatmap
     library(gplots)
     library("RColorBrewer")
     #mat <- assay(rld)
     #mm <- model.matrix(~condition, colData(rld))
     #mat <- limma::removeBatchEffect(mat, batch=rld$batch, design=mm)
     #assay(rld) <- mat
     # -- pca --
     png("pca.png", 1200, 800)
     plotPCA(rld, intgroup=c("condition"))
     dev.off()
     # -- heatmap --
     png("heatmap.png", 1200, 800)
     distsRL <- dist(t(assay(rld)))
     mat <- as.matrix(distsRL)
     hc <- hclust(distsRL)
     hmcol <- colorRampPalette(brewer.pal(9,"GnBu"))(100)
     heatmap.2(mat, Rowv=as.dendrogram(hc),symm=TRUE, trace="none",col = rev(hmcol), margin=c(13, 13))
     dev.off()

draw 3D PCA plots.

     library(gplots)
     library("RColorBrewer")

     library(ggplot2)
     data <- plotPCA(rld, intgroup=c("condition", "replicate"), returnData=TRUE)
     write.csv(data, file="plotPCA_data.csv")
     #calculate all PCs including PC3 with the following codes
     library(genefilter)
     ntop <- 500
     rv <- rowVars(assay(rld))
     select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))]
     mat <- t( assay(rld)[select, ] )
     pc <- prcomp(mat)
     pc$x[,1:3]
     #df_pc <- data.frame(pc$x[,1:3])
     df_pc <- data.frame(pc$x)
     identical(rownames(data), rownames(df_pc)) #-->TRUE

     data$PC1 <- NULL
     data$PC2 <- NULL
     merged_df <- merge(data, df_pc, by = "row.names")
     #merged_df <- merged_df[, -1]
     row.names(merged_df) <- merged_df$Row.names
     merged_df$Row.names <- NULL  # remove the "name" column
     merged_df$name <- NULL
     merged_df <- merged_df[, c("PC1","PC2","PC3","PC4","PC5","PC6","PC7","PC8","PC9","PC10","PC11","PC12","group","condition","replicate")]
     write.csv(merged_df, file="merged_df_12PCs.csv")
     summary(pc)
     #                           PC1     PC2     PC3
     #Proportion of Variance  0.5969  0.1732 0.09444

     python3 draw_3D.py
     # Edit on the generated SVG-figure.
     #/usr/bin/convert PCA_3D.png -crop 2900x1600+250+700 PCA_3D_cropped.png

(Optional) estimate size factors and dispersion values.

     #Size Factors: These are used to normalize the read counts across different samples. The size factor for a sample accounts for differences in sequencing depth (i.e., the total number of reads) and other technical biases between samples. After normalization with size factors, the counts should be comparable across samples. Size factors are usually calculated in a way that they reflect the median or mean ratio of gene expression levels between samples, assuming that most genes are not differentially expressed.
     #Dispersion: This refers to the variability or spread of gene expression measurements. In RNA-seq data analysis, each gene has its own dispersion value, which reflects how much the counts for that gene vary between different samples, more than what would be expected just due to the Poisson variation inherent in counting. Dispersion is important for accurately modeling the data and for detecting differentially expressed genes.
     #So in summary, size factors are specific to samples (used to make counts comparable across samples), and dispersion values are specific to genes (reflecting variability in gene expression).

     sizeFactors(dds)
     #NULL
     # Estimate size factors
     dds <- estimateSizeFactors(dds)
     # Estimate dispersions
     dds <- estimateDispersions(dds)
     #> sizeFactors(dds)
     #control_r1 control_r2 control_r3 VZV.d10_r1 VZV.d10_r2 VZV.d10_r3 VZV.d15_r1
     #1.2072182  1.1009757  0.8610196  1.1965615  0.8575169  1.1714809  0.8344160
     #VZV.d15_r2 VZV.d15_r3 VZV.d20_r1 VZV.d20_r2 VZV.d20_r3
     #0.5278583  1.1768657  1.1006353  1.1288974  1.2630530

     # If alone with virus data, the following BUG occured:
     #Still NULL --> BUG --> using manual calculation method for sizeFactor calculation!
                         HeLa_TO_r1                      HeLa_TO_r2
                         0.9978755                       1.1092227
     data.frame(genes = rownames(dds), dispersions = dispersions(dds))

     #Given the raw counts, the control_r1 and control_r2 samples seem to have a much lower sequencing depth (total read count) than the other samples. Therefore, when normalization methods are applied, the normalization factors for these control samples will be relatively high, boosting the normalized counts.

     1/0.9978755=1.002129023
     1/1.1092227=

     #bamCoverage --bam ../markDuplicates/${sample}Aligned.sortedByCoord.out.bam -o ${sample}_norm.bw --binSize 10 --scaleFactor  --effectiveGenomeSize 2864785220
     bamCoverage --bam ../markDuplicates/HeLa_TO_r1Aligned.sortedByCoord.out.markDups.bam -o HeLa_TO_r1.bw --binSize 10 --scaleFactor 1.002129023     --effectiveGenomeSize 2864785220
     bamCoverage --bam ../markDuplicates/HeLa_TO_r2Aligned.sortedByCoord.out.markDups.bam -o HeLa_TO_r2.bw --binSize 10 --scaleFactor  0.901532217        --effectiveGenomeSize 2864785220

     raw_counts <- counts(dds)
     normalized_counts <- counts(dds, normalized=TRUE)
     #write.table(raw_counts, file="raw_counts.txt", sep="\t", quote=F, col.names=NA)
     #write.table(normalized_counts, file="normalized_counts.txt", sep="\t", quote=F, col.names=NA)

     #convert bam to bigwig using deepTools by feeding inverse of DESeq’s size Factor
     estimSf <- function (cds){
         # Get the count matrix
         cts <- counts(cds)

         # Compute the geometric mean
         geomMean <- function(x) prod(x)^(1/length(x))

         # Compute the geometric mean over the line
         gm.mean  <-  apply(cts, 1, geomMean)

         # Zero values are set to NA (avoid subsequentcdsdivision by 0)
         gm.mean[gm.mean == 0] <- NA

         # Divide each line by its corresponding geometric mean
         # sweep(x, MARGIN, STATS, FUN = "-", check.margin = TRUE, ...)
         # MARGIN: 1 or 2 (line or columns)
         # STATS: a vector of length nrow(x) or ncol(x), depending on MARGIN
         # FUN: the function to be applied
         cts <- sweep(cts, 1, gm.mean, FUN="/")

         # Compute the median over the columns
         med <- apply(cts, 2, median, na.rm=TRUE)

         # Return the scaling factor
         return(med)
     }
     #https://dputhier.github.io/ASG/practicals/rnaseq_diff_Snf2/rnaseq_diff_Snf2.html
     #http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#data-transformations-and-visualization
     #https://hbctraining.github.io/DGE_workshop/lessons/02_DGE_count_normalization.html
     #https://hbctraining.github.io/DGE_workshop/lessons/04_DGE_DESeq2_analysis.html
     #https://genviz.org/module-04-expression/0004/02/01/DifferentialExpression/
     #DESeq2’s median of ratios [1]
     #EdgeR’s trimmed mean of M values (TMM) [2]
     #http://www.nathalievialaneix.eu/doc/html/TP1_normalization.html  #very good website!
     test_normcount <- sweep(raw_counts, 2, sizeFactors(dds), "/")
     sum(test_normcount != normalized_counts)

Merged Host-Virus Counts -> DESeq2 -> DEGs of 6 Comparisons

     Rscript complete_deg_pipeline.R

     # Key changes made:
     #     1. Virus gene identification: Added detection of virus genes using pattern matching ("^ORF|cat_RNA|redF_RNA|repE_RNA|sopA_RNA|sopB_RNA|UL45_RNA")
     #     2. Updated significance cutoff: Changed from |log2FC| >= 1 to |log2FC| >= 2
     #     3. Color scheme:
     #        * Green for virus RNAs (all virus genes)
     #        * Dark green for significantly regulated virus genes
     #        * Red for significantly up-regulated host genes
     #        * Blue for significantly down-regulated host genes
     #        * Gray for non-significant host genes
     #     4. Enhanced legend: The legend now clearly distinguishes between virus RNAs and host genes with their regulation status
     #     5. Summary statistics: Added detailed output showing counts of virus vs host genes and their regulation status

     #  ================================================================================
     # 📊 FINAL SUMMARY OF ALL 6 COMPARISONS
     #                 name total  up down virus sig_total pct_sig
     # 03_VZV.d20_vs_control 17594 193  220    77       413     2.3
     # 02_VZV.d15_vs_control 17594 139   94    77       233     1.3
     # 05_VZV.d20_vs_VZV.d10 17594  35  114    77       149     0.8
     # 01_VZV.d10_vs_control 17594  89   20    77       109     0.6
     # 06_VZV.d20_vs_VZV.d15 17594  28   38    77        66     0.4
     # 04_VZV.d15_vs_VZV.d10 17594   6   13    77        19     0.1

Protected: 自闭症和强迫症

Enter your password to view comments.

GOOD_QUESTIONS

Leave a reply

Validation of Mutant-Specific Effects The additional methodological details on mutant construction and genome validation are appreciated. However, the study still lacks functional complementation or equivalent validation to demonstrate that the observed phenotypes are directly attributable to loss of the targeted efflux pumps rather than secondary regulatory consequences. Reliance on knockout mutants alone limits confidence in causal interpretation, particularly for broad transcriptomic and stress-response phenotypes. 突变体特异性效应的验证感谢作者补充了关于突变体构建和基因组验证的额外方法学细节。然而，本研究仍缺乏功能回补或等效的验证手段，来证明所观察到的表型是直接由靶向外排泵的缺失引起的，而非继发性的调控后果。仅依赖敲除突变体限制了因果推断的可靠性，尤其是对于广泛的转录组与应激反应表型而言。

93只恒生指数成分股名单

Leave a reply

根据维基百科和恒生指数公司的最新信息，从2026年6月8日起，恒生指数成分股已增加至93只[[2]]。以下是完整的93只恒生指数成分股名单，按行业分类：

金融业（10只成份股）

汇丰控股 (0005)
友邦保险 (1299)
建设银行 (0939)
工商银行 (1398)
香港交易所 (0388)
中国平安 (2318)
中国银行 (3988)
中国人寿 (2628)
招商银行 (3968)
中银香港 (2388)

非必需性消费及必需性消费（30只成份股）

阿里巴巴-SW (9988)
美团-W (3690)
比亚迪股份 (1211)
京东集团-SW (9618)
百度集团-SW (9888)
创科实业 (0669)
快手-W (1024)
吉利汽车 (0175)
泡泡玛特 (9992)
安踏体育 (2020)
携程集团-S (9961)
农夫山泉 (9633)
理想汽车-W (2015)
万洲国际 (0288)
银河娱乐 (0027)
港铁公司 (0066)
美的集团 (0300)
蒙牛乳业 (2319)
海尔智家 (6690)
华润啤酒 (0291)
李宁 (2331)
申洲国际 (2313)
金沙中国 (1928)
老铺黄金 (6181)
新东方-S (9901)
海底捞 (6862)
康师傅控股 (0322)
周大福 (1929)
恒安国际 (1044)
百威亚太 (1876)

资讯科技业（6只成份股）

腾讯控股 (0700)
小米集团-W (1810)
中芯国际 (0981)
网易-S (9999)
联想集团 (0992)
比亚迪电子 (0285)

能源业、原材料业、工业及综合事业（18只成份股）

中国海洋石油 (0883)
中国石油股份 (0857)
紫金矿业 (2899)
长和 (0001)
中国神华 (1088)
宁德时代 (3000)
中国宏桥 (1378)
中国石油化工股份 (0386)
中信股份 (0267)
中通快递 (2057)
洛阳钼业 (3993)
极兔速递 (1519) 新增
舜宇光学科技 (2382)
中国铝业 (2600) 新增
京东物流 (2618)
信义玻璃 (0868)
东方海外国际 (0316)
信义光能 (0968)

电讯业及公用事业（9只成份股）

中国移动 (0941)
中电控股 (0002)
电能实业 (0006)
香港中华煤气 (0003)
中国电信 (0728)
中国联通 (0762)
新奥能源 (2688)
华润电力 (0836)
长江基建集团 (1038)

地产建筑业（10只成份股）

新鸿基地产 (0016)
华润置地 (1109)
领展房产基金 (0823)
长实集团 (1113)
中国海外发展 (0688)
恒基地产 (0012)
九龙仓置业 (1997)
华润万象生活 (1209)
龙湖集团 (0960)
恒隆地产 (0101)

医疗保健业（10只成份股）

百济神州 (6160) 新增
信达生物 (1801)
药明生物 (2269)
药明康德 (2359)
石药集团 (1093)
中国生物制药 (1177)
翰森制药 (3692)
京东健康 (6618)
阿里健康 (0241)
国药控股 (1099)

最新变动： 在2026年5月22日的季度检讨中，恒生指数公司宣布新增3只蓝筹股：极兔速递(1519)、中国铝业(2600) 和 百济神州(6160)，令恒指成分股数目由90只增加至93只，并于2026年6月8日起生效[[3]][[7]]。

在港股市场中，股票名称后面的 -W、-S 甚至 -SW、-B 等字母后缀，是港交所为了区分公司的特殊股权结构或上市类型而设置的“身份标签”。

这能帮助投资者一眼看出这家公司的治理结构或上市背景。具体含义如下：

1. `-W`：同股不同权 (Weighted Voting Rights)

全称：Weighted Voting Rights (WVR)。
含义：指这家公司采用了“同股不同权”的架构。在这种架构下，公司的股票分为普通股份和特别表决权股份。普通投资者买到的股票（1股只有1票投票权），而创始人或管理层持有的特别股份（1股拥有多票投票权，比如1股抵10票）。
目的：确保创始团队在公司经过多轮融资、自身股权被稀释的情况下，依然能牢牢掌握公司的控制权和经营决策权（这在互联网和科技公司中非常常见）。
代表公司：美团-W (3690)、小米集团-W (1810)。

2. `-S`：第二上市 (Secondary Listing)

全称：Secondary Listing。
含义：指这家公司已经在海外（如美国纽交所、纳斯达克）完成了“主要上市”，然后再次来到港交所发行股票并挂牌交易。
目的：这类公司通常是为了拓宽融资渠道、吸引亚洲及大中华区的投资者、提升股票流动性，或者为了规避单一海外市场的监管风险（近年来很多中概股回港都属于此类）。
代表公司：百度集团-S (9888)。
(注：如果一家第二上市的公司后来申请转为“双重主要上市”，它的 -S 标记通常会被取消。)

3. 组合后缀：`-SW`

含义：如果您看到像 阿里巴巴-SW (9988) 或 京东集团-SW (9618) 这样的名字，说明这家公司既是第二上市（-S），又采用了同股不同权架构（-W）。

补充：其他常见后缀

除了 -W 和 -S，您在港股名单中可能还会看到：

-B：代表 Biotech（生物科技公司）。
- 含义：指这家公司是未通过港交所常规盈利测试、但具备核心研发能力的生物科技公司。因为还在研发阶段，可能尚未实现盈利，带有 -B 是交易所特意提醒投资者注意其未盈利的风险。
- 代表公司：歌礼制药-B (1672)。

总结： 这些字母后缀只是交易所的提示性标记，相当于公司的“特殊身份证”，它们不会影响股票的正常交易，也不直接代表公司基本面的好坏。普通投资者在买卖时，只需知道它们代表了特殊的股权或上市背景即可。

不是的，比亚迪并非只在香港上市。它是一家典型的 “A+H”两地上市的公司，同时在中国大陆的A股市场和香港的H股市场挂牌交易。

具体情况如下：

1. A股（中国大陆）

上市地点：深圳证券交易所（深交所）
股票简称：比亚迪
股票代码：002594
上市时间：2011年6月（比亚迪最早是在A股上市的，当时发行的是中小板，后来中小板与主板合并）。

2. H股（中国香港）

上市地点：香港联合交易所（港交所）
股票简称：比亚迪股份
股票代码：01211（在港股软件中常显示为 1211）
上市时间：2002年7月（比亚迪最早其实是先在港股上市的，比A股早了将近9年）。

为什么会产生“只在港股上市”的错觉？

在您之前看到的“恒生指数成分股”名单中，列出的是 “比亚迪股份 (1211)”。因为恒生指数是港股的指数，所以名单里只能收录它在港股上市的代码和名称。这可能会让人误以为它只有港股。

A股和H股是什么关系？

同一家公司：002594（A股）和 01211（H股）代表的是完全同一家公司（比亚迪股份有限公司），享有同等的股东权益（同股同权）。
价格不同：因为A股和H股是在两个完全独立的市场（大陆市场和香港市场）进行交易，面对的投资者群体、资金流动性、市场情绪都不同，所以两者的股价（折算成同一种货币后）通常是不一样的。这就产生了所谓的“A/H股溢价”。

补充：在其他地方上市了吗？

除了中国大陆（A股）和中国香港（H股）之外，比亚迪没有在美国（如纳斯达克、纽交所）或欧洲（如伦敦交易所）等其他海外资本市场上市。

总结：如果您在大陆的股票软件里搜索“比亚迪”（002594），或者在港股软件里搜索“比亚迪股份”（01211），买到的都是同一家公司的股票，只是交易的币种（人民币 vs 港币）和所在的市场不同而已。

这个问题需要分两个层面来回答：

如果按大众认知和实际市场主导地位来看，美国主要有 2个核心股票交易所；但如果按美国证券交易委员会（SEC）的官方注册来看，美国目前共有 16个 全国性的合法证券交易所。

为了让您清晰了解，我们可以将它们分为“两大巨头”、“其他专业交易所”和“场外市场”三个梯队：

第一梯队：两大绝对巨头（占据 99% 以上的市值和交易量）

这是普通投资者最熟悉、也是全球最著名的两家交易所。绝大多数知名的美国上市公司都在这两家挂牌：

纽约证券交易所 (NYSE – New York Stock Exchange)
- 特点：历史最悠久（1792年成立），位于华尔街11号。采用“专家撮合”制度（虽然现在很多也电子化了），门槛极高。
- 代表公司：主要是传统行业巨头、超级蓝筹股、大型金融机构。比如：伯克希尔·哈撒韦、沃尔玛、可口可乐、摩根大通、强生等。
纳斯达克 (Nasdaq)
- 特点：成立于1971年，是全球第一个全电子化的股票交易市场。没有实体的交易大厅，门槛相对纽交所略低，对科技和创新企业非常友好。
- 代表公司：全球顶尖的科技巨头、成长型公司。比如：苹果、微软、英伟达、特斯拉、亚马逊、Meta（脸书）、谷歌等。

第二梯队：其他全国性证券交易所（共 14 个左右）

除了纽交所和纳斯达克，美国SEC还注册了十几个其他的证券交易所。这些交易所普通散户平时很少直接感知，它们主要处理机构订单、提供流动性或主打特定理念。 它们包括：

纽交所旗下的其他交易所：
- NYSE Arca： originally 太平洋证券交易所，现在是ETF（交易型开放式指数基金） 交易的绝对主力。
- NYSE American：原美国证券交易所（AMEX），现在主要面向中小型市值公司和初创企业。
- NYSE Chicago / NYSE National：主要处理部分电子订单，提供交易通道。
芝加哥期权交易所 (Cboe) 旗下的股票交易所：
- Cboe 本身是期权交易巨头，但它旗下也拥有 4 个股票交易所（Cboe BZX, BYX, EDGA, EDGX）。很多券商（如Robinhood、盈透）会把散户的订单路由到这些交易所去撮合。
特色/独立交易所：
- IEX (Investors Exchange)：因畅销书《闪电侠》而闻名。它的特色是故意设置一个微小的“减速带”（350微秒的延迟），用来防止高频交易机构利用速度优势“插队”收割普通投资者。
- MEMX (Members Exchange)：由华尔街多家大型券商（如嘉信理财、摩根士丹利等）联合出资成立的交易所，目的是为了打破纽交所和纳斯达克的垄断，降低交易手续费。
- LTSE (Long-Term Stock Exchange)：长期证券交易所，由硅谷著名风投家提出，旨在奖励长期持有的股东（比如持股时间越长，投票权越大），对抗华尔街的短期逐利行为。

补充：场外交易市场 (OTC – Over-The-Counter)

除了上述正规的“交易所”，美国还有一个庞大的场外交易市场（严格意义上它不是交易所，而是一个交易商网络）。

OTC Markets：大家常听说的 “粉单市场” (Pink Sheets) 就在这里。
特点：这里上市的公司通常达不到纽交所或纳斯达克的财务标准。里面充斥着仙股（几毛钱甚至几分钱的股票）、退市公司、破产重组公司，以及一些不想承担高昂合规成本的微型企业。风险极高。

总结

如果您问的是 “美国最主要的股票交易场所”，答案是 2个（纽交所和纳斯达克）。
如果您问的是 “美国官方认可的证券交易所数量”，答案是 16个。美国之所以有这么多交易所，是因为美国有一个“全国市场系统（NMS）”，允许各家交易所通过降低手续费、提供更快的网速或更公平的机制来互相竞争，从而让整体的交易成本保持在极低的水平。

以下是截至2026年6月最新的纳斯达克100指数（NASDAQ-100）完整成分股列表及权重[[1]]。

重要说明：虽然指数名为“100”，但由于Alphabet（谷歌母公司） 同时发行了A类（GOOGL）和C类（GOOG）两种股票，这两类股票均被单独计入指数，因此下表实际包含101个条目，代表的是100家公司。权重数据基于各公司市值在指数总市值中的占比动态计算[[1]]。

纳斯达克100指数完整成分股及权重表

排名	公司中文名称	公司英文名称	股票代码	权重
1	英伟达	Nvidia Corp	NVDA	12.74%
2	苹果	Apple Inc.	AAPL	11.11%
3	微软	Microsoft Corp	MSFT	7.19%
4	亚马逊	Amazon.Com Inc	AMZN	6.55%
5	谷歌A类股	Alphabet Inc. Class A	GOOGL	5.83%
6	谷歌C类股	Alphabet Inc. Class C	GOOG	5.44%
7	博通	Broadcom Inc.	AVGO	4.89%
8	特斯拉	Tesla, Inc.	TSLA	3.80%
9	Meta (Facebook母公司)	Meta Platforms, Inc.	META	3.70%
10	美光科技	Micron Technology, Inc.	MU	3.12%
11	沃尔玛	Walmart Inc.	WMT	2.39%
12	超威半导体 (AMD)	Advanced Micro Devices	AMD	2.19%
13	阿斯麦 (ASML)	ASML Holding NV	ASML	1.88%
14	英特尔	Intel Corp	INTC	1.67%
15	泛林集团	Lam Research Corp	LRCX	1.25%
16	应用材料	Applied Materials Inc	AMAT	1.25%
17	思科	Cisco Systems, Inc.	CSCO	1.19%
18	Arm控股	Arm Holdings plc	ARM	1.18%
19	开市客 (Costco)	Costco Wholesale Corp	COST	1.09%
20	奈飞 (Netflix)	NetFlix Inc	NFLX	0.83%
21	科磊 (KLA)	KLA Corporation	KLAC	0.83%
22	帕兰泰尔 (Palantir)	Palantir Technologies Inc.	PLTR	0.80%
23	闪迪 (SanDisk)	Sandisk Corporation	SNDK	0.77%
24	德州仪器	Texas Instruments Incorporated	TXN	0.71%
25	美满电子 (Marvell)	Marvell Technology, Inc.	MRVL	0.68%
26	西部数据	Western Digital Corp.	WDC	0.66%
27	希捷科技	Seagate Technology Holdings PLC	STX	0.64%
28	林德气体	Linde plc	LIN	0.61%
29	高通	Qualcomm Inc	QCOM	0.59%
30	派拓网络	Palo Alto Networks, Inc.	PANW	0.59%
31	亚德诺半导体 (ADI)	Analog Devices, Inc.	ADI	0.52%
32	T-Mobile	T-Mobile US, Inc.	TMUS	0.50%
33	百事可乐	PepsiCo, Inc.	PEP	0.49%
34	安进	Amgen Inc	AMGN	0.47%
35	CrowdStrike	CrowdStrike Holdings, Inc.	CRWD	0.44%
36	AppLovin	Applovin Corporation	APP	0.41%
37	吉利德科学	Gilead Sciences Inc	GILD	0.40%
38	霍尼韦尔	Honeywell International, Inc.	HON	0.37%
39	直观外科	Intuitive Surgical Inc.	ISRG	0.37%
40	Shopify	Shopify Inc.	SHOP	0.36%
41	Booking Holdings	Booking Holdings Inc.	BKNG	0.34%
42	福泰制药	Vertex Pharmaceuticals Inc	VRTX	0.29%
43	星巴克	Starbucks Corp	SBUX	0.29%
44	拼多多	PDD Holdings Inc.	PDD	0.29%
45	楷登电子	Cadence Design Systems	CDNS	0.28%
46	飞塔 (Fortinet)	Fortinet, Inc.	FTNT	0.27%
47	万豪国际	Marriott International	MAR	0.27%
48	星座能源	Constellation Energy Corporation	CEG	0.24%
49	怪物饮料	Monster Beverage Corporation	MNST	0.23%
50	新思科技	Synopsys Inc	SNPS	0.23%
51	自动数据处理公司 (ADP)	Automatic Data Processing	ADP	0.22%
52	CSX运输	CSX Corporation	CSX	0.22%
53	爱彼迎	Airbnb, Inc.	ABNB	0.21%
54	美卡多	Mercado Libre, Inc	MELI	0.21%
55	康卡斯特	Comcast Corp	CMCSA	0.21%
56	Datadog	Datadog, Inc.	DDOG	0.20%
57	亿滋国际	Mondelez International, Inc.	MDLZ	0.20%
58	奥多比 (Adobe)	Adobe Inc.	ADBE	0.20%
59	恩智浦半导体	NXP Semiconductors N.V.	NXPI	0.20%
60	罗斯百货	Ross Stores Inc	ROST	0.19%
61	单片电源系统 (MPS)	Monolithic Power Systems, Inc.	MPWR	0.19%
62	奥莱利汽车	O’Reilly Automotive, Inc.	ORLY	0.19%
63	Intuit	Intuit Inc	INTU	0.19%
64	DoorDash	DoorDash, Inc.	DASH	0.18%
65	美国电力公司 (AEP)	American Electric Power Company, Inc.	AEP	0.18%
66	Lumentum	Lumentum Holdings Inc.	LITE	0.18%
67	Cintas	Cintas Corp	CTAS	0.17%
68	华纳兄弟探索	Warner Bros. Discovery, Inc.	WBD	0.17%
69	再生元制药	Regeneron Pharmaceuticals Inc	REGN	0.16%
70	帕卡 (Paccar)	Paccar Inc	PCAR	0.16%
71	贝克休斯	Baker Hughes Company	BKR	0.15%
72	微芯科技	Microchip Technology Inc	MCHP	0.13%
73	响尾蛇能源	Diamondback Energy, Inc.	FANG	0.13%
74	艺电 (EA)	Electronic Arts Inc	EA	0.13%
75	快扣 (Fastenal)	Fastenal Co	FAST	0.13%
76	卓越能源	Xcel Energy, Inc.	XEL	0.13%
77	法罗里奥 (Ferrovial)	Ferrovial N.V.	FER	0.13%
78	爱克斯龙 (Exelon)	Exelon Corporation	EXC	0.12%
79	老道明货运	Old Dominion Freight Line	ODFL	0.12%
80	爱德士实验室 (IDEXX)	Idexx Laboratories Inc	IDXX	0.11%
81	可口可乐欧洲太平洋	Coca-Cola Europacific Partners plc	CCEP	0.11%
82	Take-Two Interactive	Take-Two Interactive Software Inc	TTWO	0.11%
83	Keurig Dr Pepper	Keurig Dr Pepper Inc.	KDP	0.11%
84	Strategy Inc (原MicroStrategy)	Strategy Inc	MSTR	0.11%
85	欧特克 (Autodesk)	Autodesk Inc	ADSK	0.10%
86	Alnylam制药	Alnylam Pharmaceuticals, Inc.	ALNY	0.10%
87	贝宝 (PayPal)	PayPal Holdings, Inc.	PYPL	0.10%
88	Paychex	Paychex Inc	PAYX	0.09%
89	汤森路透	Thomson Reuters Corporation	TRI	0.09%
90	Axon Enterprise	Axon Enterprise, Inc.	AXON	0.09%
91	柔佛科技 (Roper)	Roper Technologies, Inc.	ROP	0.08%
92	Workday	Workday, Inc.	WDAY	0.08%
93	GE医疗	GE HealthCare Technologies Inc.	GEHC	0.07%
94	德康医疗 (DexCom)	DexCom, Inc.	DXCM	0.07%
95	Copart	Copart Inc	CPRT	0.07%
96	卡夫亨氏	The Kraft Heinz Company	KHC	0.07%
97	Verisk Analytics	Verisk Analytics, Inc.	VRSK	0.06%
98	高知特 (Cognizant)	Cognizant Technology Solutions	CTSH	0.06%
99	Insmed	Insmed, Inc.	INSM	0.05%
100	Zscaler	Zscaler, Inc.	ZS	0.05%
101	特许通讯	Charter Comm Inc Del CL A New	CHTR	0.04%

注：由于股票市值每日随市场波动，上述权重数据为实时动态估算值，仅供参考。

关于Betano送你80欧元优惠的“流水要求”（Umsatzbedingungen）

Leave a reply

关于Betano送你80欧元、需要充值80欧元的这个优惠，这里的关键信息都整理好了：

💰 这80欧元能提现吗？

可以，但有条件。

这80欧元属于奖金（Bonusguthaben），不能直接提现。你需要先满足 “流水要求”（Umsatzbedingungen） 后才能提现。

具体要求是：你需要用 “充值金额+奖金总额” （即80+80=160欧元）下注，且总下注额达到这个数的5倍，也就是 800欧元。

并且，这些投注最低赔率必须是1.65。只有满足这些条件后，奖金及其产生的盈利才会转为可提现的余额。

⏳ 有使用期限吗？

有，而且有两个不同的截止日期需要注意：

完成流水要求的期限：你有90天的时间来完成上述800欧元的流水要求。如果90天后没完成，奖金和相关的盈利可能会被收回。
激活后投注的期限：奖金发放后，也需要在一定时间内用于下注，通常是90天内。

📝 其他重要规则

这是一项新用户优惠：通常只针对新注册用户。
最低充值额：要获得这个100%的奖金，最低充值额是10欧元。
注意支付方式：并非所有支付方式都符合条件，例如使用 Skrill 充值可能无法获得此奖金。
某些投注类型不算：不是所有投注都算在流水里，比如“双 chance”等特殊玩法可能不被计入。
可能还有额外赠礼：这个优惠有时还附带一张20欧元的免费投注券（Freebet）。它通常有单独的7-14天有效期，赢了只给你净利润，不包括本金。

💎 总结

简单来说，Betano的这笔奖金可以提现，但必须在90天内，用不低于1.65的赔率，完成总计800欧元的投注额。

建议你在参与前，务必登录Betano官网，仔细阅读最新的 “奖金条款”（Bonusbedingungen），因为具体规则可能会有调整。

How to Convert GenBank Files to GTF/FASTA for the nf-core/rnaseq Pipeline (Data_Nina_RNAseq_2026)

Leave a reply

Abstract: When running the nf-core/rnaseq pipeline on custom or non-model viral genomes (like VZV or HSV), standard GenBank (.gb) files often cause indexing errors in STAR or quantification failures in Salmon due to missing GTF attributes. In this post, I share a robust Python (Biopython) workflow to convert GenBank files into perfectly formatted FASTA and GTF files. We specifically address how to inject required features like exon, gene_type, and properly formatted transcript_id and gene_id attributes to ensure seamless execution of the star_salmon alignment and quantification steps.

You are absolutely right! The duplicate exon was not the root cause. The real issue is that Salmon and the nf-core/rnaseq pipeline strictly expect specific attributes (like gene_type "protein_coding") to be present in every corresponding record to correctly build the transcript-to-gene mapping and index.

By comparing your original BB1528_nanopore_consensus.gtf with the HSV1 structure, we can see the missing “additional IDs/attributes”:

gene_type "protein_coding" is missing in all records.
The transcript_id format should ideally end with _RNA (e.g., UL33_RNA) instead of _tx.
CDS records should have an orf_id.
gene records should not have gene_name, and exon records should not have exon_number (to match the clean HSV1 structure).

Instead of re-converting from the GenBank file, we can write a quick Python script to directly fix your existing BB1528_nanopore_consensus.gtf and reformat its attributes to perfectly match the HSV1 structure.

Step 1: Save and Run the GTF Formatting Script

Save the following code as fix_gtf_to_hsv1_format.py in your working directory:

import sys

def fix_gtf(input_gtf, output_gtf):
    seen_exons = set()

    with open(input_gtf, 'r') as fin, open(output_gtf, 'w') as fout:
        for line in fin:
            if line.startswith('#'):
                continue
            parts = line.strip().split('\t')
            if len(parts) < 9:
                continue

            seqname, source, feature, start, end, score, strand, frame, attributes = parts

            # Parse existing attributes into a dictionary
            attrs_dict = {}
            for attr in attributes.split(';'):
                attr = attr.strip()
                if not attr:
                    continue
                if ' ' in attr:
                    key, val = attr.split(' ', 1)
                    attrs_dict[key] = val.strip('"')

            gene_id = attrs_dict.get('gene_id', 'unknown_gene')
            # Standardize transcript_id to gene_id_RNA to match HSV1 format
            new_transcript_id = f"{gene_id}_RNA"

            # Rebuild attributes strictly matching HSV1 structure
            new_attrs = ""
            if feature == "gene":
                new_attrs = f'gene_id "{gene_id}"; gene_type "protein_coding"'
            elif feature == "transcript":
                new_attrs = f'gene_id "{gene_id}"; transcript_id "{new_transcript_id}"; gene_type "protein_coding"'
            elif feature == "exon":
                # Avoid duplicate exons
                exon_key = (seqname, start, end, strand, gene_id, new_transcript_id)
                if exon_key in seen_exons:
                    continue
                seen_exons.add(exon_key)
                new_attrs = f'gene_id "{gene_id}"; transcript_id "{new_transcript_id}"; gene_type "protein_coding"'
            elif feature == "CDS":
                orf_id = f"{gene_id}_CDS"
                new_attrs = f'gene_id "{gene_id}"; transcript_id "{new_transcript_id}"; orf_id "{orf_id}"; gene_type "protein_coding"'
            else:
                # For any other features, keep them but ensure gene_type is present
                if 'transcript_id' in attrs_dict:
                    attrs_dict['transcript_id'] = new_transcript_id
                if 'gene_type' not in attrs_dict:
                    attrs_dict['gene_type'] = 'protein_coding'
                new_attrs = '; '.join([f'{k} "{v}"' for k, v in attrs_dict.items()])

            # Write the fixed line
            fout.write(f"{seqname}\t{source}\t{feature}\t{start}\t{end}\t{score}\t{strand}\t{frame}\t{new_attrs}\n")

if __name__ == "__main__":
    fix_gtf("BB1528_nanopore_consensus.gtf", "BB1528_final.gtf")
    print("✅ GTF successfully reformatted to match HSV1 structure -> BB1528_final.gtf")

Run it in your terminal:

python3 fix_gtf_to_hsv1_format.py

Step 2: Verify the Output

If you check the new BB1528_final.gtf, the structure will now perfectly mirror your HSV1 example:

BB1528  GenBank gene    88  561 .   +   .   gene_id "ORF_S_L_transcribed"; gene_type "protein_coding"
BB1528  GenBank transcript  88  561 .   +   .   gene_id "ORF_S_L_transcribed"; transcript_id "ORF_S_L_transcribed_RNA"; gene_type "protein_coding"
BB1528  GenBank exon    88  561 .   +   .   gene_id "ORF_S_L_transcribed"; transcript_id "ORF_S_L_transcribed_RNA"; gene_type "protein_coding"
BB1528  GenBank CDS 88  561 .   +   0   gene_id "ORF_S_L_transcribed"; transcript_id "ORF_S_L_transcribed_RNA"; orf_id "ORF_S_L_transcribed_CDS"; gene_type "protein_coding"

Step 3: Re-run Nextflow with `-resume`

Now, update your Nextflow command to point to the newly formatted BB1528_final.gtf.

Important: Keep the -resume flag! Because the GTF file’s content has changed, Nextflow will automatically detect the new file hash and re-run only the genome preparation steps (like SALMON_INDEX and GTF2BED) without re-processing your FASTQ files.

nextflow run rnaseq/main.nf \
    --input samplesheet.csv \
    --outdir results_BB1528 \
    --fasta BB1528_nanopore_consensus.fasta \
    --gtf BB1528_final.gtf \
    --with_umi \
    --umitools_extract_method regex \
    --umitools_bc_pattern '^(?P
<UMI>.{12}).*' \
    --umitools_dedup_stats \
    -profile test_full \
    -resume \
    --max_memory 256.GB \
    --max_time 2400.h \
    --save_reference \
    --aligner star_salmon \
    --gtf_extra_attributes gene_id \
    --gtf_group_features transcript_id \
    --featurecounts_group_type gene_id \
    --featurecounts_feature_type transcript \
    --skip_rseqc \
    --skip_dupradar \
    --skip_preseq \
    --skip_biotype_qc \
    --skip_deseq2_qc \
    --skip_multiqc \
    --min_mapped_reads 0

This should completely resolve the SALMON_QUANT indexing error!

Good news: Your GTF and FASTA files are perfectly formatted!

We can prove this from your log: [info] Index contained 78 targets This means Salmon successfully built the transcriptome index and found all 78 transcripts from your BB1528_final.gtf. The reference files are completely correct.

What is causing the error?

The error is purely a biological/sample-specific issue:

[warning] salmon was only able to assign 9 fragments to transcripts in the index, 
but the minimum number of required assigned fragments (--minAssignedFrags) was 10.

The nf-core/rnaseq pipeline automatically subsamples your reads (in this case, 1,000,000 fragments) to run a quick Salmon quantification and check the library strandedness (--libType=A).

Because you are mapping to a specific viral reference (BB1528), it is highly likely that control_r3 is a negative control or a sample with an extremely low viral load. Out of the 1,000,000 subsampled reads, only 9 reads actually belonged to the virus.

Salmon has a built-in safety threshold (--minAssignedFrags 10) that intentionally crashes the pipeline if fewer than 10 reads map to the transcriptome, to prevent “garbage” quantification of empty samples. Since 9 < 10, the pipeline stops.

How to fix it

You can easily bypass this by creating a custom Nextflow configuration file that tells Salmon to accept samples with even just 1 mapped fragment.

Step 1: Create a config file In your working directory, create a file named salmon_fix.config and add the following code:

process {
    withName: 'SALMON_QUANT' {
        // Append --minAssignedFrags 1 to any existing Salmon arguments
        ext.args = { (task.ext.args ?: '') + " --minAssignedFrags 1" }
    }
}

Step 2: Run the pipeline with the config file Add -c salmon_fix.config to your Nextflow command. This will apply the fix to both the subsampled QC step and the main quantification step, ensuring your negative controls don’t crash the pipeline.

nextflow run /home/jhuang/Tools/nf-core-rnaseq-3.12.0/main.nf -profile docker \
    -c salmon_fix.config \
    --input samplesheet.csv \
    --outdir results_BB1528 \
    --fasta BB1528_nanopore_consensus.fasta \
    --gtf BB1528_final.gtf \
    --with_umi \
    --umitools_extract_method regex \
    --umitools_bc_pattern '^(?P
<umi_1>.{12}).*' \
    --trimmer fastp \
    --extra_fastp_args "--trim_front2 10" \
    -resume \
    --max_cpus 54 \
    --max_memory 120.GB \
    --max_time 2400.h \
    --save_align_intermeds \
    --save_unaligned \
    --save_reference \
    --aligner star_salmon \
    --gtf_extra_attributes gene_id \
    --gtf_group_features transcript_id \
    --featurecounts_group_type gene_id \
    --featurecounts_feature_type transcript \
    --skip_rseqc --skip_dupradar --skip_preseq --skip_biotype_qc --skip_deseq2_qc --skip_multiqc \
    --min_mapped_reads 0

Summary of what will happen now:

The pipeline will resume from where it left off (thanks to -resume).
Salmon will process control_r3, see that it only has 9 viral reads, but instead of crashing, it will output a warning and successfully generate an output file with 0 or near-0 counts.
The rest of your infected samples (which likely have plenty of viral reads) will quantify normally.
At the end of the pipeline, you will get a complete count matrix where your negative controls correctly show 0 expression, which is exactly what you want for downstream differential expression analysis!

Yes, all detected configuration files are being loaded and merged by Nextflow. Your custom settings will definitely be applied. What does this warning actually mean? This is not an error; it is just a best-practice reminder built into the nf-core pipeline framework. The developers of nf-core strictly enforce a rule: Pipeline parameters (the arguments starting with — like –input, –fasta, –max_memory) should never be defined inside a .config file. They must be passed directly via the command line (CLI) or a dedicated -params-file. The warning is simply the pipeline saying:

"Hey, I see you are using custom .config files. Just a reminder: do not put your --parameters inside them! Use the command line for parameters, and use .config files only for Nextflow settings (like process resources, arguments, or executor settings)."

Is your setup correct? Yes, absolutely. The salmon_fix.config file we created in the previous step only contains a process block to modify Salmon’s internal arguments (ext.args). It does not contain any params. Therefore, you are using the config file exactly as intended. Why did it trigger? Nextflow automatically detects any file ending in .config in your working directory (e.g., nextflow.config, salmon_fix.config, or profile configs). Because it found more than one, it triggered this informational warning. What should you do? You can safely ignore this warning. The pipeline will continue running normally, and your fix to allow Salmon to quantify samples with fewer than 10 assigned fragments (–minAssignedFrags 1) will be successfully applied. Did the pipeline successfully pass the SALMON_QUANT step and finish running this time? If it didn’t crash at that specific step, it means the config file worked perfectly!

Based on the error log and the reference files you provided, here is the diagnosis of the problem and how to fix it.

🔍 The Root Cause

The pipeline is failing during the FASTQ_SUBSAMPLE_FQ_SALMON step. nf-core/rnaseq uses this step to subsample 1 million reads and run Salmon to automatically detect your library’s strandedness (e.g., forward, reverse, unstranded) and verify the mapping rate.

The critical error from Salmon is:

salmon was only able to assign 8 fragments to transcripts in the index, but the minimum number of required assigned fragments (--minAssignedFrags) was 10. This could be indicative of a mismatch between the reference and sample...

Why is this happening?

Your Reference is Virus-Only: The reference you provided (BB1528_nanopore_consensus.fasta) is a Varicella-Zoster Virus (VZV) genome (~127 kb).
Your Controls are Likely Host-Only: You have samples named control_r1, control_r2, and control_r3. These are almost certainly uninfected host cells (e.g., human or monkey cell lines).
The Mismatch: Because host RNA-seq reads do not match the VZV viral genome, the mapping rate for your control samples is essentially 0% (only 8 random reads mapped out of 1,000,000). Salmon aborts because it cannot confidently determine the library type with fewer than 10 mapped reads.

🛠️ How to Fix It

Depending on your experimental goals, choose one of the following solutions:

Solution 1: Add the Host Genome to your Reference (Recommended)

If you want to analyze host gene expression or simply want the pipeline to run smoothly on uninfected controls, you must include the host genome (e.g., Human GRCh38 or Macaque) alongside the VZV genome.

Action: Concatenate the Host FASTA and VZV FASTA into a single combined_genome.fa. Do the same for the Host GTF and VZV GTF (combined_annotation.gtf).
Why: This allows the host reads in your control samples to map properly, allowing Salmon to easily detect the library strandedness and pass the QC threshold.

Solution 2: Remove the Control Samples from the Samplesheet

If you are only interested in viral transcription and intentionally built a virus-only reference, the uninfected controls will yield zero useful data anyway.

Action: Open your samplesheet.csv and delete the rows for control_r1, control_r2, and control_r3.
Why: The pipeline will only run on the VZV.dXX infected samples, which will have high mapping rates to the VZV reference, allowing Salmon to successfully detect the library type.

Solution 3: Manually Specify Strandedness

If you know the strandedness of your library preparation kit (e.g., Illumina TruSeq Stranded mRNA is usually reverse), you can tell the pipeline to skip the auto-detection step.

Action: Add a strandedness column to your samplesheet.csv and fill it with reverse, forward, or unstranded for all samples. Alternatively, pass --strandedness reverse (or appropriate value) in your Nextflow run command.
Note: While this bypasses the library-type inference failure, Salmon might still throw a warning about the low mapping rate for the controls, but it usually prevents the hard crash.

Solution 4: Verify Sample Identities

If your control samples are actually supposed to be VZV-infected (e.g., a different viral strain or a positive control), then an 8-read mapping rate indicates a severe issue.

Action: Check your FASTQ files. Ensure the control FASTQ files aren’t corrupted, empty, or accidentally swapped with a completely different organism’s sequencing data. Ensure the viral strain in the control isn’t highly divergent from the Oka parental strain reference you are using.

💡 Summary Recommendation

For viral RNA-seq experiments, Solution 1 is the standard best practice. Creating a “chimeric” reference (Host + Virus) ensures that mapping rates are accurate, host-response pathways can be evaluated, and nf-core/rnaseq’s internal QC tools function exactly as intended.

Comprehensive Reproducible Pipeline for Longitudinal Nasal Microbiome and S. epidermidis Epidome Analysis (DATA_B/Data_Luise_Epidome_longitudinal_nose)

Leave a reply

readme

phyloseq-rmd

microbiotaprocess-r

1. Study Design & Metadata

This pipeline analyzes longitudinal nasal swab samples from patients undergoing neurosurgery.

Cohort A (Aneurysm): 15 patients (45 samples)
Cohort H (Hypophysis): 20 patients (60 samples)
Total: 108 patient samples + 3 extraction controls.
Timepoints:
- .1 = Admission (Baseline)
- .2 = Surgery (Intraoperative/Immediate post-op)
- .3 = Discharge (Recovery)

Targeted Sequencing Approaches:

16S rRNA Gene Amplicon Sequencing: For overall nasal microbiome profiling (processed via QIIME1 open-reference picking against SILVA 132).
Epidome Method (S. epidermidis): Targeted amplicon sequencing of the g216 and yycH genes for high-resolution Sequence Type (ST) tracking (processed via DADA2 and custom Epidome Python/R scripts).

2. R Environment & Dependencies

To reproduce this analysis, ensure you are using R (v4.1+) and have the following packages installed.

# Core Microbiome & Phyloseq
install.packages(c("phyloseq", "microbiome", "picante", "vegan", "ape"))
# Bioconductor packages
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("DESeq2", "microbiome", "phyloseq"))

# Advanced Processing & Visualization
install.packages(c("MicrobiotaProcess", "microeco", "ggplot2", "ggpubr", "ggrepel", "aplot"))
install.packages(c("ggh4x", "gghalves", "ggalluvial", "RColorBrewer", "heatmaply", "gplots"))

# Statistical Testing & Data Wrangling
install.packages(c("dplyr", "tidyr", "rstatix", "RVAideMemoire", "openxlsx", "knitr", "kableExtra", "tibble", "reshape2", "coin", "BSDA"))

3. Workflow Part I: Overall Microbiome (16S) Analysis

Scripts involved: MicrobiotaProcess.R (Primary engine for diversity/ordination) & Phyloseq.Rmd (Reporting, DESeq2, and specific statistical tests).

3.1 Data Import & Normalization Strategy

The pipeline splits the raw OTU table (table_even9893.biom) into distinct objects based on the downstream statistical requirements:

ps_filt: Raw counts, low-depth samples removed (min 1000 reads). Used for DESeq2 and Hellinger transformation.
ps_rarefied: Rarefied to an even depth of 9,893 reads (set.seed(9242)). Used for Alpha diversity and UniFrac/Jaccard metrics.
ps_abund_rel: Relative abundance, filtered to keep only taxa with mean relative abundance > 0.1%. Used for clean taxonomic composition plots.
hellinger (via MicrobiotaProcess): Hellinger-transformed ps_filt. Used for Bray-Curtis distances and PCoA.

3.2 Alpha Diversity (Longitudinal & Cross-Sectional)

Metrics: Observed OTUs, Chao1, ACE, Shannon, Simpson, Pielou.
Longitudinal Analysis (Within Cohorts A and H):
- Overall temporal shift: Friedman test (non-parametric repeated measures) accounting for PatientID pairing.
- Post-hoc pairwise: Paired t-tests and Paired Wilcoxon signed-rank tests with Benjamini-Hochberg (BH) correction.
Cross-Sectional Analysis (Cohort A vs. Cohort H):
- All timepoints pooled: Independent t-test / Wilcoxon rank-sum test.
- Surgery timepoint ONLY: To prevent dilution of the surgical effect by baseline/discharge samples, a specific filter (SampleType == "surgery") is applied before running the independent t-test/Wilcoxon.

3.3 Beta Diversity & Ordination

Distance Metric: Bray-Curtis dissimilarity calculated on Hellinger-transformed non-rarefied counts.
Ordination: Principal Coordinates Analysis (PCoA).
Statistical Testing:
- Overall: PERMANOVA (vegan::adonis2) with 9,999 permutations.
- Post-hoc Pairwise: All 15 possible pairwise group comparisons (e.g., A1 vs A2, A1 vs H1, H2 vs H3) are calculated using adonis2. P-values are adjusted using BH (FDR) and Bonferroni corrections.
Outputs: PCoA plots (colored by group, sized by Shannon, alpha by Observed), exported as PNG/PDF/SVG. Pairwise PERMANOVA results exported to Bray_pairwise_PERMANOVA.xlsx.

3.4 Taxonomic Composition

Plots: Stacked barplots and Heatmaps at the Class level (Top 20 taxa).
Grouping: Plots generated for individual samples, grouped by Cohort (A vs H), and grouped by Timepoint (Admission, Surgery, Discharge).

3.5 Differential Abundance Analysis (DESeq2)

Method: Negative Binomial Generalized Linear Model (Wald test) using raw, non-rarefied integer counts (ps_filt).
Comparisons: Six independent pairwise comparisons per cohort:
1. Admission vs. Surgery
2. Admission vs. Discharge
3. Surgery vs. Discharge (Repeated for both Cohort A and Cohort H)
Threshold: Adjusted P-value (BH) < 0.05.
Outputs: Volcano-style plots (log2FoldChange by Family/Class) and Excel tables of significant DEGs for each comparison.

4. Workflow Part II: S. epidermidis Epidome Analysis

Script involved: Phyloseq.Rmd (Sections: ST summary, Binary Prevalence Analysis). Input Data: count_table_seq31_seq37_ST.txt (ASVs classified into STs via g216 and yycH BLAST against the Epidome DB). Normalized to median depth (56,191).

4.1 ST Composition & Alpha Diversity

Composition: Stacked barplots of relative ST abundance across the 108 samples.
Alpha Diversity: Shannon index and Observed STs compared between Cohorts (t-test) and across Timepoints (t-test).
goeBURST: Clonal complex (CC) visualization based on PubMLST allele profiles.

4.2 Continuous Abundance Testing (Specific STs)

Method: Wilcoxon signed-rank test (paired, Admission .1 vs. Discharge .3).
Targets: Analyzed individually for clinically relevant STs (ST2, ST5, ST8, ST23, ST59, ST60, ST73, ST100, ST130, ST215) and in combined clinical groups (e.g., ST5+ST87+ST130+ST210+ST384).
Handling Ties: Evaluated using exact calculations via the coin package and the Sign Test (BSDA) to ensure robustness against zero-inflation.

4.3 🌟 Novel Binary Prevalence Analysis (Cochran’s Q & McNemar)

Because continuous abundance of specific high-risk STs (ST2, ST5, ST23) showed no significant shifts, we converted the data to a binary presence/absence matrix to track strain acquisition and clearance.

Data Transformation: Abundance > 0 becomes 1 (Present); becomes (Absent).
Overall Temporal Shift: Cochran’s Q test (RVAideMemoire::cochran.qtest) to evaluate if the prevalence proportion changes across the 3 timepoints (Admission, Surgery, Discharge) for a given ST.
Pairwise Transitions: Robust McNemar tests for paired binary data (Admission vs Surgery, Surgery vs Discharge, Admission vs Discharge).
- Robustness: Custom R function handles zero-variance (no discordant pairs) by assigning $p = 1.0$ instead of crashing.
Correction: Benjamini-Hochberg (BH) FDR correction applied to the 3 pairwise p-values.
Visualization: Line charts with scatter points showing the Prevalence (%) trend across the 3 timepoints, annotated with $n/N$ (number of positive patients / total patients).
Export: The binary matrix is exported to ST_binary_matrix.xlsx.

5. Execution Instructions (How to Reproduce)

Step 1: Prepare the Directory

Ensure your working directory contains the following raw files:

table_even9893.biom and ../clustering/rep_set.tre (16S data)
../map3_corrected.txt (Metadata)
../count_table_seq31_seq37_ST.txt (Epidome ST counts)
adiv_even.txt, adiv_even_A.txt, adiv_even_H.txt (Pre-calculated 16S alpha metrics from QIIME)
alpha_diversity_metrics_samples_ST.csv (Pre-calculated ST alpha metrics)

Step 2: Run the MicrobiotaProcess Pipeline (Beta/Alpha/Composition)

Open your terminal or RStudio and run the MicrobiotaProcess.R script. This will generate the figures_MP/ and figures_All_Combined/ directories.

Rscript MicrobiotaProcess.R

Outputs: PCoA plots, rarefaction curves, Bray-Curtis distance matrices, pairwise PERMANOVA Excel files, and Class-level heatmaps.

Step 3: Render the Comprehensive HTML Report

Open RStudio, load Phyloseq.Rmd, and render it. This integrates the outputs from Step 2, runs DESeq2, performs the Epidome statistical tests, and generates the final interactive report.

# In R Console:
rmarkdown::render('Phyloseq.Rmd', output_file = 'Phyloseq.html')

Outputs: Phyloseq.html (the master report), figures/ directory (DESeq2 plots, ST prevalence plots, goeBURST), and various .xlsx files for ST binary matrices and DEG tables.

6. Summary of Key Output Files

File / Directory	Description
`Phyloseq.html`	Master Report. Contains all 16S and Epidome statistics, tables, and embedded figures.
`figures_MP/`	High-res PCoA, rarefaction, alpha diversity, and composition heatmaps from `MicrobiotaProcess`.
`figures_All_Combined/`	Pairwise PERMANOVA results (`Bray_pairwise_PERMANOVA.xlsx`) and combined PCoA plots.
`figures/`	DESeq2 volcano plots, ST prevalence trend lines, and Alpha diversity boxplots.
`DEGs_*.xlsx`	Excel tables of differentially abundant OTUs for all 6 pairwise comparisons per cohort.
`ST_binary_matrix.xlsx`	The 0/1 presence/absence matrix for ST2, ST5, and ST23 used for the McNemar analysis.

Pipeline maintained by Jiabin Huang. For questions regarding the Epidome Python scripts or upstream QIIME1/DADA2 processing, refer to the respective bash/R scripts in the epidome/scripts/ directory.

Prerequisites

Step 1: Load Data and Create Metadata

Step 2: Reproduce Figures from wmgx_report.pdf

1. Ordination (PCoA)

2. Heatmap (Top 25 Species)

3. Stacked Barplot (Top 15 Species)

Step 3: Differential Abundance Analysis (Pre vs. Post)

Step 4: Is the difference between time points significant?

How to interpret the adonis2 output:

1. Why Enrichment is Based on KOs

2. Why the Co-author Saw a “Discrepancy”

3. The Correct Solution

The Corrected R Function

💡 How to Reply to Your Co-author

💡 改进说明：

金融业（10只成份股）

非必需性消费及必需性消费（30只成份股）

资讯科技业（6只成份股）

能源业、原材料业、工业及综合事业（18只成份股）

电讯业及公用事业（9只成份股）

地产建筑业（10只成份股）

医疗保健业（10只成份股）

1. -W：同股不同权 (Weighted Voting Rights)

2. -S：第二上市 (Secondary Listing)

3. 组合后缀：-SW

补充：其他常见后缀

1. A股（中国大陆）

2. H股（中国香港）

为什么会产生“只在港股上市”的错觉？

A股和H股是什么关系？

补充：在其他地方上市了吗？

第一梯队：两大绝对巨头（占据 99% 以上的市值和交易量）

第二梯队：其他全国性证券交易所（共 14 个左右）

补充：场外交易市场 (OTC – Over-The-Counter)

总结

纳斯达克100指数完整成分股及权重表

💰 这80欧元能提现吗？

⏳ 有使用期限吗？

📝 其他重要规则

💎 总结

Step 1: Save and Run the GTF Formatting Script

Step 2: Verify the Output

Step 3: Re-run Nextflow with -resume

What is causing the error?

How to fix it

Summary of what will happen now:

🔍 The Root Cause

🛠️ How to Fix It

Solution 1: Add the Host Genome to your Reference (Recommended)

Solution 2: Remove the Control Samples from the Samplesheet

Solution 3: Manually Specify Strandedness

Solution 4: Verify Sample Identities

💡 Summary Recommendation

1. Study Design & Metadata

2. R Environment & Dependencies

3. Workflow Part I: Overall Microbiome (16S) Analysis

3.1 Data Import & Normalization Strategy

3.2 Alpha Diversity (Longitudinal & Cross-Sectional)

3.3 Beta Diversity & Ordination

3.4 Taxonomic Composition

3.5 Differential Abundance Analysis (DESeq2)

4. Workflow Part II: S. epidermidis Epidome Analysis

4.1 ST Composition & Alpha Diversity

4.2 Continuous Abundance Testing (Specific STs)

4.3 🌟 Novel Binary Prevalence Analysis (Cochran’s Q & McNemar)

5. Execution Instructions (How to Reproduce)

Step 1: Prepare the Directory

Step 2: Run the MicrobiotaProcess Pipeline (Beta/Alpha/Composition)

Step 3: Render the Comprehensive HTML Report

6. Summary of Key Output Files

Step 2: Reproduce Figures from `wmgx_report.pdf`

How to interpret the `adonis2` output:

1. `-W`：同股不同权 (Weighted Voting Rights)

2. `-S`：第二上市 (Secondary Listing)

3. 组合后缀：`-SW`

Step 3: Re-run Nextflow with `-resume`