这两个检验是统计学中专门用于处理配对设计（Paired Design）或重复测量设计（Repeated Measures Design）下二分类数据（Binary Data，如是/否、阳性/阴性、检出/未检出） 的经典非参数方法。

在您的课题中，同一个患者在三个时间点（入院、手术、出院）被反复采样，这构成了天然的“配对/纵向”关系。合作者建议将数据转换为“0/1”矩阵，正是为了使用这套统计体系。

以下为您全面、深入地解析这两个检验的原理、逻辑以及在您课题中的具体应用。

一、 McNemar 检验 (McNemar’s Test)

1. 核心定位 用于比较两个相关（配对）样本的二分类频率是否有显著差异。 在您的课题中，用于比较两个时间点（例如：入院 vs. 手术）某种 ST（如 ST2）的检出率是否有变化。

2. 核心逻辑：只关注“发生变化的样本” 很多初学者会误以为 McNemar 检验是比较两组的总阳性率，其实不然。它的核心逻辑是：如果两个时间点没有差异，那么“由阴转阳”的人数应该等于“由阳转阴”的人数。

假设我们有 $N$ 个患者，比较 Time 1（入院）和 Time 2（手术），数据可以整理成如下的 2×2 配对列联表：

• $a$ 和 $d$：代表状态没有改变的患者。它们对判断“时间点之间是否有差异”不提供任何信息。
• $b$ 和 $c$：代表状态发生改变的患者（不一致对子，Discordant pairs）。
• McNemar 检验的本质：就是检验 $b$ 和 $c$ 是否相等。如果 $b$ 远大于 $c$，说明手术导致了 ST 的清除或丢失；如果 $c$ 远大于 $b$，说明手术导致了 ST 的新发获得。

3. 统计量计算 在 $b+c$ 足够大（通常 $>25$）时，服从卡方分布： $$ \chi^2 = \frac{(b – c)^2}{b + c} $$ (注：如果 $b+c < 25$，卡方近似不准确，必须使用 McNemar’s Exact Test 精确检验，基于二项分布计算。)

二、 Cochran’s Q 检验 (Cochran’s Q Test)

1. 核心定位 McNemar 检验只能比较两个时间点。当您的研究有三个或三个以上相关时间点（入院、手术、出院）时，如果您做 3 次 McNemar 检验，会大幅增加假阳性率（Multiple testing problem）。 Cochran’s Q 检验就是 McNemar 检验在 3 个及以上时间点的扩展，用于评估多个时间点的整体检出率是否存在差异。

2. 核心逻辑 它类似于配对设计的 ANOVA（重复测量方差分析），但针对的是二分类数据。它通过计算每个时间点的总阳性数，与“期望的总阳性数”之间的偏差，构建一个服从卡方分布的 $Q$ 统计量。

3. 严格的前提假设（避坑指南）

• 数据必须完整（Complete Block）：Cochran’s Q 检验要求每个受试者在所有时间点都有数据。如果某个患者只有入院和出院数据，缺失了手术数据，传统的 Cochran’s Q 会报错或将其整行剔除。
• 样本量要求：患者数量不能太少，否则检验效能（Power）极低。

三、 针对您课题的标准统计分析流程

面对 ST2、ST5、ST23 这三个目标 ST，以及 3 个时间点，严谨的统计学报告流程应如下：

第一步：整体评估 (Omnibus Test)

针对每一个 ST（例如 ST2），运行 Cochran’s Q 检验。

• 零假设 ($H_0$)：ST2 在入院、手术、出院三个时间点的检出率完全相同。
• 目的：回答“ST2 的定植状态在整体病程中是否发生了波动？”
• 报告：如果 $p 0.05$，通常认为无整体差异，但您仍可以基于临床假设进行第二步（事后比较）。

第二步：事后两两比较 (Post-hoc Pairwise Comparisons)

使用 McNemar 检验（或精确 McNemar 检验）进行三组两两对比：

1. Admission vs. Surgery (入院 vs. 手术)
2. Surgery vs. Discharge (手术 vs. 出院)
3. Admission vs. Discharge (入院 vs. 出院)

第三步：多重检验校正 (Multiple Testing Correction)

因为您对 ST2 做了 3 次 McNemar 检验，假阳性率膨胀了。必须使用 Benjamini-Hochberg (BH) 方法对得到的 3 个原始 p 值进行校正，得到 $p_{adj}$ (FDR)。

• 只有当 $p_{adj} < 0.05$ 时，才能宣称这两个时间点之间的检出率差异具有统计学显著性。

四、 为什么合作者的建议非常专业且切中要害？

在微生物组学和流行病学中，“丰度（Abundance）” 和 “检出率/流行率（Prevalence）” 代表两个完全不同的生物学维度：

1. 克服“零膨胀”与“极端值”干扰： 表皮葡萄球菌的某些 ST（如 ST2）在很多患者鼻腔中可能是极低丰度甚至检测不到的（大量 0 值）。如果您用绝对丰度做 Wilcoxon 或 t 检验，结果极易被个别“超级携带者”（某患者 ST2 丰度突然飙升至 50%）拉偏，导致统计失效。 转换为 0/1 矩阵后，直接忽略了丰度大小的干扰，纯粹考察“有没有”的问题，统计效能（Power）往往更高。
2. 直接对接临床意义（院感控制）： ST2、ST5 是国际公认的与医院获得性感染（HAI）、生物膜形成、多重耐药（如 MRSE） 高度相关的高危克隆系。 在临床视角下，患者是否“定植（Colonized）”了高危 ST，比定植了多少个细菌更关键。合作者建议分析这三个时间点的检出比例（Prevalence %），实际上是在回答一个极具临床价值的流行病学问题：
   • “外科手术干预（及住院环境）是否增加了患者获得院内高危 ST 克隆（如 ST2/ST5）的风险？” （如果 Surgery/Discharge 的检出率显著高于 Admission）。
   • “手术或抗生素使用是否清除了患者原有的定植？” （如果检出率显著下降）。

五、 进阶建议：如果数据有缺失值怎么办？

如前所述，Cochran’s Q 检验对缺失值非常苛刻。在真实的临床纵向研究中，患者脱落或样本采集失败导致的数据缺失很常见。

如果您的数据中存在部分患者缺失某个时间点（例如只有 .1 和 .3，没有 .2），除了我上一轮提供的代码中使用的 tryCatch 容错机制外，您在论文的方法学部分可以提及使用更强大的广义线性混合模型 (GLMM, Generalized Linear Mixed Models) 作为补充或替代。

• GLMM 的优势：它可以处理二分类数据（二项分布 family=binomial），将“时间点”作为固定效应（Fixed effect），将“患者 ID”作为随机截距（Random intercept, (1|PatientID)）。
• 最大优点：GLMM 采用最大似然估计（Maximum Likelihood），允许数据存在缺失（Missing at Random），不需要像 Cochran’s Q 那样剔除任何不完整配对的样本，在顶级期刊的纵向队列研究中更受审稿人青睐。

总结：合作者的建议将您的分析从“描述性的微生物群落变化”提升到了“具有明确临床流行病学意义的病原体定植动态研究”。使用 Cochran’s Q + 配对 McNemar + BH校正 的组合拳，在统计学上无懈可击，能够非常漂亮地回应审稿人可能提出的质疑。

Direct Answer

Yes and No.

• Yes, you can use biobakery_workflows wmgx to process all your FASTQ files into taxonomic and functional profiles, and wmgx_vis to generate basic visualizations (like alpha/beta diversity).
• No (Not Recommended), you should not use wmgx_vis for the actual differential analysis. Because your study is a paired/longitudinal design (Pre- vs. Post-treatment from the same subjects across multiple time points), you must account for the paired nature of the data (using random effects). The command-line wrapper for wmgx_vis does not support specifying random effects, which would lead to statistically flawed results.

Instead, the standard practice in the bioBakery ecosystem is to use wmgx for processing, and then run MaAsLin 2 directly for the differential analysis.

Here is the complete workflow on how to do this:

Step 1: Process the data with wmgx

Since your paired-end files follow the standard _1.fastq.gz and _2.fastq.gz naming convention, wmgx will automatically detect and pair them.

Assuming all your symlinks are in a folder called /data/fastq_files/:

biobakery_workflows wmgx \
  --input /data/fastq_files/ \
  --output /data/wmgx_results/ \
  --threads 20

(Note: If you don’t need strain profiling, you can speed this up by adding --bypass-strain-profiling)

Step 2: Basic Visualization with wmgx_vis (Optional)

You can use wmgx_vis to generate standard diversity plots and heatmaps. You will need a metadata file (see Step 3 for the format).

biobakery_workflows wmgx_vis \
  --input /data/wmgx_results/ \
  --project-name Metagenome_Study \
  --input-metadata /data/metadata.tsv \
  --metadata-categorical Group Timepoint \
  --output /data/visualization_results/

Step 3: Proper Differential Analysis (Using MaAsLin 2 directly)

To correctly compare Pre vs. Post while accounting for the fact that the samples are paired (and potentially accounting for the different time points), you should use MaAsLin 2 directly on the output tables generated by wmgx.

1. Prepare your Metadata File (metadata.tsv)

You need a tab-separated file where the SampleID exactly matches the prefix of your FASTQ files (e.g., 2025_Nov_Pre_A1).

(Note: You must correctly map which “Pre” sample belongs to which “Post” sample in the Subject_ID column. This is how the model knows they are paired).

2. Run MaAsLin 2 in R

You can run MaAsLin 2 (which is installed with bioBakery) using R. This allows you to define Group (Pre/Post) as a fixed effect and Subject_ID as a random effect, which is mathematically required for paired data.

# Install Maaslin2 if you haven't already
# BiocManager::install("Maaslin2")
library(Maaslin2)

# Run the differential abundance analysis
fit_data <- Maaslin2(
    input_data     = "/data/wmgx_results/wmgx/humann/genefamilies.tsv", # or pathabundance.tsv
    input_metadata = "/data/metadata.tsv",
    output         = "/data/maaslin2_results/",

    # Define your variables
    fixed_effects  = c("Group", "Timepoint"), # Comparing Pre vs Post, while controlling for Timepoint
    random_effects = c("Subject_ID"),         # CRUCIAL: Accounts for the paired Pre/Post design!

    # Data processing parameters
    normalization  = "NONE",                  # HUMAnN output is already normalized to relative abundance
    transform      = "LOG",
    min_abundance  = 0.01,
    min_prevalence = 0.10,
    correction     = "BH"                     # Benjamini-Hochberg for multiple testing (FDR)
)

Summary

1. Use wmgx to process your raw FASTQs.
2. Use wmgx_vis if you just want quick, automated PCA/PCoA plots and heatmaps.
3. Do not use wmgx_vis for the differential statistics. Extract the genefamilies.tsv or metaphlan merged tables and run MaAsLin 2 directly so you can include Subject_ID as a random effect to properly handle your Pre/Post paired design.

由于2026年美加墨世界杯扩军到了48支球队，总共有104场比赛，预测所有小组赛的精确比分篇幅会非常庞大。如果你是在参加公司或朋友间的世界杯竞猜游戏（Tippspiel），我为你整理了一份最核心、最完整的赛事比分预测表。

这份预测涵盖了小组赛焦点战、完整的淘汰赛路径（从16强到决赛）以及附加竞猜问题，你可以直接参考填入你的竞猜表：

⚽ 一、 小组赛焦点战预测 (Gruppenphase)

在小组赛阶段，强队通常会稳扎稳打，比分不会过于夸张。

🇺🇸🇨🇦🇲🇽 东道主首战：

• 🇲🇽 墨西哥 vs 南非 🇿🇦 ➔ 2:0
• 🇨🇦 加拿大 vs 波黑 🇧🇦 ➔ 1:0
• 🇺🇸 美国 vs 巴拉圭 🇵🇾 ➔ 2:1

🏆 夺冠热门及强队小组赛典型比分：

• 🇩🇪 德国 vs 厄瓜多尔/澳大利亚 ➔ 3:1
• 🇪🇸 西班牙 vs 克罗地亚/新西兰 ➔ 2:0
• 🇫🇷 法国 vs 丹麦/加拿大 ➔ 2:1
• 🇦🇷 阿根廷 vs 澳大利亚/秘鲁 ➔ 2:0
• 🏴󠁧󠁢󠁥󠁮󠁧󠁿 英格兰 vs 塞尔维亚/斯洛文尼亚 ➔ 2:0
• 🇧🇷 巴西 vs 苏格兰/摩洛哥 ➔ 2:1

🏆 二、 完整淘汰赛路径预测 (K.O.-Phase)

淘汰赛阶段通常防守更严密，1球或2球的优势是常态，部分比赛会拖入点球大战。

🥊 1/8 决赛 (Achtelfinale) – 核心场次

• 🇪🇸 西班牙 2:0 日本 🇯🇵
• 🇩🇪 德国 2:1 墨西哥 🇲🇽
• 🇫🇷 法国 3:0 波兰 🇵🇱
• 🇦🇷 阿根廷 2:1 澳大利亚 🇦🇺
• 🏴󠁧󠁢󠁥󠁮󠁧󠁿 英格兰 1:0 哥伦比亚 🇨🇴
• 🇧🇷 巴西 2:0 瑞士 🇨🇭
• 🇵🇹 葡萄牙 1:0 乌拉圭 🇺🇾
• 🇳🇱 荷兰 2:1 美国 🇺🇸 (东道主遗憾止步16强)

🔥 1/4 决赛 (Viertelfinale)

• 🇪🇸 西班牙 2:1 德国 🇩🇪 (焦点大战)
• 🇫🇷 法国 2:1 阿根廷 🇦🇷 (卫冕冠军出局)
• 🏴󠁧󠁢󠁥󠁮󠁧󠁿 英格兰 1:1 (点球 4:3) 巴西 🇧🇷
• 🇵🇹 葡萄牙 2:1 荷兰 🇳🇱

🚀 半决赛 (Halbfinale)

• 🇪🇸 西班牙 2:1 法国 🇫🇷
• 🏴󠁧󠁢󠁥󠁮󠁧󠁿 英格兰 2:1 葡萄牙 🇵🇹

🥉 季军战 (Spiel um Platz 3)

• 🇫🇷 法国 2:0 葡萄牙 🇵🇹

🏆 决赛 (Finale) – 新泽西大都会人寿体育场

• 🇪🇸 西班牙 2 : 1 英格兰 🏴󠁧󠁢󠁥󠁮󠁧󠁿 (预测：西班牙凭借更强大的中场控制力（如罗德里、佩德里）和边路爆点（亚马尔）在常规时间或加时赛中绝杀英格兰，队史第二次捧杯！)

🎯 三、 附加竞猜问题 (Rahmenwetten)

如果你的竞猜表有这些问题，可以直接填这些答案：

1. 🏆 世界杯冠军 (Weltmeister): 西班牙 (Spanien)
2. 🥇 最佳射手/金靴 (Torschützenkönig): 基利安·姆巴佩 (Kylian Mbappé) – 预测进球数：6球
3. 🌟 最佳年轻球员 (Bester junger Spieler): 拉明·亚马尔 (Lamine Yamal / 西班牙)
4. 🧤 最佳门将 (Bester Torwart): 埃米利亚诺·马丁内斯 (Emiliano Martínez / 阿根廷) 或 乌奈·西蒙 (Unai Simón / 西班牙)
5. ⚽ 决赛比分 (Ergebnis des Finales): 2:1
6. 📊 赛事总进球数 (Anzahl der Turniertore): 268球 (48队104场比赛，场均约2.5-2.8球)

💡 提示 (Tipp)： 在填写竞猜表时，淘汰赛阶段的比分尽量多写 1:0, 2:0, 2:1, 1:1 这种小比分，这在世界杯淘汰赛中出现的概率远高于大比分（如3:2或4:1）。祝你的竞猜表取得好成绩（Viel Glück beim Tippspiel）！

complete_deg_pipeline.R

complete_deg_pipeline_custom_cutoff.R

1.R

1. 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

2. 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

3. (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!

4. Preparing samplesheet.csv

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

5. 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

6. 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’

7. Import data and 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_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)
   
    # -----------------------------------------------------------------
    # ---- 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")
    )
   
    # ==========================================
    # 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!")

8. 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, we use complete_deg_pipeline_custom_cutoff.R
   
    # 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

9. 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

Explanation for the “near-1” frequencies and why these weren’t reported before

Based on the data, here’s how to explain this:

Key Points:

1. These are NOT newly acquired mutations during passaging

   • All three positions (15458, 22636, 24781) show frequency = 1.0 in the starting virus (hCoV229E_Rluc)
   • They remain fixed at 1.0 throughout all passages (p10, p16, p26)
   • These represent pre-existing differences between your experimental virus strain and the reference genome PP810610

2. Why they weren’t reported in earlier analyses:

   • VPhaser2 is designed to detect intra-host variants (positions with heterogeneity, frequency < 1.0)
   • When a mutation is fixed at 100% (frequency = 1.0), it becomes part of the consensus sequence
   • Many variant callers either:
     • Don’t report fixed variants as “intra-host variants”
     • Filter them out as they represent consensus differences from reference, not within-host diversity
   • Your earlier analysis likely only reported positions with true intra-host heterogeneity (frequency between 0 and 1)

3. The “near-1” values (0.997, 0.996, 0.978) in X7523_p26:

   These slight deviations from 1.0 are NOT calculation errors but can be explained by:

   • Sequencing errors: Even at high coverage, there’s a small error rate (~0.1-1%)
   • Mapping artifacts: Some reads may map incorrectly at the position
   • Very low-level mixed population: Possibly <3% wild-type contamination
   • Statistical noise: At very high frequencies, small absolute numbers of alternative reads can cause slight deviations

Suggested Response:

You could add this to your email:

Zusätzliche Erklärung zu den Positionen 15458, 22636 und 24781:

Bei genauerer Betrachtung der Daten zeigt sich, dass diese drei Mutationen bereits im Ausgangsvirus (hCoV229E_Rluc) zu 100% vorhanden waren und sich während des gesamten Passagierens nicht verändert haben. Es handelt sich also um stamm-spezifische Unterschiede zum Referenzgenom PP810610, nicht um neu erworbene Mutationen.

Warum wurden sie in früheren Analysen nicht berichtet?

• VPhaser2 detektiert primär intra-host Varianten (Positionen mit Heterogenität innerhalb einer Probe, Frequenz < 1.0)
• Bei einer Frequenz von exakt 1.0 (100%) werden diese Positionen als Teil der Konsensus-Sequenz betrachtet und oft nicht als “Varianten” im engeren Sinne berichtet
• In der früheren Analyse wurden nur Positionen mit echter intra-host Diversität aufgeführt

Die leicht abweichenden Werte in X7523_p26 (0.997, 0.996, 0.978): Diese minimalen Abweichungen von 1.0 sind technisch bedingt durch:

• Sequenzierfehler (typischerweise 0.1-1% Fehlerrate)
• Mapping-Artefakte
• Statistisches Rauschen bei sehr hohen Frequenzen

Sie stellen keine biologisch relevante Heterogenität dar, sondern liegen im Bereich der technischen Variabilität.

This explanation is scientifically accurate and shows that you’ve thoroughly investigated the issue!

https://huttenhower.sph.harvard.edu/biobakery_workflows/

Whole metagenome shotgun sequencing data can be processed through read-level quality control (KneadData), taxonomic profiling (MetaPhlAn), functional profiling (HUMAnN), and strain profiling (StrainPhlAn) to generate a report with publication-ready figures with two workflow commands.

1. Prepare the toy datasets

    # Install if needed: conda install -c bioconda seqtk
    cd ~/DATA/Data_Tam_Metagenomics_2026_pre_vs_post_treatment/X101SC25123808-Z01-J002/01.RawData/A
    seqtk sample -s100 A_1.fq.gz 0.01 | gzip > ../J002_A_1.fastq.gz
    seqtk sample -s100 A_2.fq.gz 0.01 | gzip > ../J002_A_2.fastq.gz
    cd ../B
    seqtk sample -s100 B_1.fq.gz 0.01 | gzip > ../J002_B_1.fastq.gz
    seqtk sample -s100 B_2.fq.gz 0.01 | gzip > ../J002_B_2.fastq.gz
    mv J002*.fastq.gz /mnt/md1/DATA/Data_Tam_Metagenomics_2026_wastewater/X101SC25123808-Z01-J003/01.RawData/A_test/
   
    seqtk sample -s100 B_1.fastq.gz 0.01 | gzip > ../A_test/B_1.fastq.gz
    seqtk sample -s100 B_2.fastq.gz 0.01 | gzip > ../A_test/B_2.fastq.gz

2. 拉取镜像（注意：latest 实际是 2019-2021 年构建的旧版）

    docker pull biobakery/workflows:latest

3. 验证容器内版本

    docker run --rm biobakery/workflows:latest biobakery_workflows --version

4. Install Databases Inside Container

    # Create persistent host directory for databases
    mkdir -p /mnt/nvme4n1p1/biobakery_db
   
    docker run -it \
    -v /mnt/nvme4n1p1/biobakery_db:/biobakery_databases \
    biobakery/workflows:latest \
    /bin/bash
   
    # Inside container:
    biobakery_workflows_databases --available
    #There are five available database sets each corresponding to a data processing workflow.
    #wmgx: The full databases for the whole metagenome workflow
    #wmgx_demo: The demo databases for the whole metagenome workflow
    #wmgx_wmtx: The full databases for the whole metagenome and metatranscriptome workflow
    #16s_usearch: The full databases for the 16s workflow
    #16s_dada2: The full databases for the dada2 workflow
    #16s_its: The unite database for the its workflow
    #isolate_assembly: The eggnog-mapper databases for the assembly workflow
   
    biobakery_workflows_databases --install wmgx_demo --location /biobakery_databases
    biobakery_workflows_databases --install wmgx_wmtx --location /biobakery_databases
    biobakery_workflows_databases --install 16s_usearch --location /biobakery_databases
    biobakery_workflows_databases --install 16s_dada2 --location /biobakery_databases
    biobakery_workflows_databases --install 16s_its --location /biobakery_databases
    biobakery_workflows_databases --install isolate_assembly --location /biobakery_databases
   
    biobakery_workflows_databases --install wmgx --location /biobakery_databases
    # ---- DOWNLOAD_LOG ----
    1. INSTALLING humann utility mapping database
        Creating directory to install database: /biobakery_databases/humann
   
        Creating subdirectory to INSTALL database: /biobakery_databases/humann/utility_mapping
        Download URL: http://huttenhower.sph.harvard.edu/humann2_data/full_mapping_v201901.tar.gz
        Downloading file of size: 2.55 GB
        2.55 GB 100.00 %   4.70 MB/sec  0 min -0 sec
        Extracting: /biobakery_databases/humann/full_mapping_v201901.tar.gz
        Database installed: /biobakery_databases/humann/utility_mapping
        HUMAnN configuration file updated: database_folders : utility_mapping = /biobakery_databases/humann/utility_mapping
   
        Generating strainphlan fasta database (FOR GENERETING DIRS strainphlan_db_reference and strainphlan_db_markers from humann/utility_mapping?), it is contradicted with the following assumption: bowtie2-inspect ${DB_DIR}/metaphlan_databases/mpa_vJan25_CHOCOPhlAnSGB_202503 > ${DB_DIR}/strainphlan_db_markers/all_markers.fasta
   
    2. INSTALLING humann nucleotide and protein databases
   
        Creating subdirectory to INSTALL database: /biobakery_databases/humann/chocophlan
        Download URL: http://huttenhower.sph.harvard.edu/humann2_data/chocophlan/full_chocophlan.v296_201901.tar.gz
        Downloading file of size: 15.30 GB
        15.30 GB 100.00 %   6.75 MB/sec  0 min -0 sec
        Extracting: /biobakery_databases/humann/full_chocophlan.v296_201901.tar.gz
        Database installed: /biobakery_databases/humann/chocophlan
        HUMAnN configuration file updated: database_folders : nucleotide = /biobakery_databases/humann/chocophlan
   
        Creating subdirectory to INSTALL database: /biobakery_databases/humann/uniref
        Download URL: http://huttenhower.sph.harvard.edu/humann2_data/uniprot/uniref_annotated/uniref90_annotated_v201901.tar.gz
        Downloading file of size: 19.31 GB
        19.31 GB 100.00 %   7.22 MB/sec  0 min -0 sec
        Extracting: /biobakery_databases/humann/uniref90_annotated_v201901.tar.gz
        Database installed: /biobakery_databases/humann/uniref
        HUMAnN configuration file updated: database_folders : protein = /biobakery_databases/humann/uniref
   
    3. INSTALLING hg kneaddata database
   
        Creating directory to install database: /biobakery_databases/kneaddata_db_human_genome
        Download URL: http://huttenhower.sph.harvard.edu/kneadData_databases/Homo_sapiens_hg37_and_human_contamination_Bowtie2_v0.1.tar.gz
        Downloading file of size: 3.48 GB
        3.48 GB 100.00 %   7.17 MB/sec  0 min -0 sec
        Extracting: /biobakery_databases/kneaddata_db_human_genome/Homo_sapiens_hg37_and_human_contamination_Bowtie2_v0.1.tar.gz
        Database installed: /biobakery_databases/kneaddata_db_human_genome
   
    A custom install location was selected. Please set the environment variable $BIOBAKERY_WORKFLOWS_DATABASES to the install location.

5. DEBUGs

   5.1. Unable to find fastqc

    # Install the missing fastqc software
   
        apt-get update
        apt-get install -y fastqc

   5.2. Install Java 11 to correctly run fastqc

    # Install Java 11
        apt-get install -y openjdk-11-jre-headless

   5.3. Wrong version of kneaddata

    #🔍BUG: 从你提供的日志中可以看到两件事：你当前安装的 kneaddata 版本是 v0.7.10。biobakery_workflows 在调用 kneaddata 时，强行传入了 --run-trf 这个参数。
    #然而，在 kneaddata v0.7.4 及以后的版本中，TRF（串联重复序列过滤）已经变成了默认开启的功能，因此开发者移除了 --run-trf 这个命令行参数（只保留了 --bypass-trf 用于跳过它）。
    #因为 biobakery_workflows 的脚本里还写死了要传递 --run-trf，而 v0.7.10 的 kneaddata 根本不认识这个参数，所以直接报错退出：unrecognized arguments: --run-trf。随后，由于第一步的 kneaddata 失败了，依赖它的所有下游任务（MetaPhlAn, HUMAnN 等）也随之全部级联失败。
    #方案一：降级 kneaddata 到兼容版本（推荐）
    #我们需要将 kneaddata 降级到 0.7.3 版本，这是最后一个原生支持 --run-trf 参数的稳定版本，且能与当前的 biobakery_workflows 完美配合。
   
    pip install --upgrade kneaddata
    #Successfully installed kneaddata-0.12.4 (NOT_COMPATIBLE)
   
    pip install kneaddata==0.7.10
    #Successfully installed kneaddata-0.7.10
    root@13192f2ad6e6:/data# /usr/local/bin/kneaddata --version
    kneaddata v0.7.10 (NOT_COMPATIBLE)
   
    pip uninstall -y kneaddata
    pip install kneaddata==0.7.3
    /usr/local/bin/kneaddata --version
    # 应该输出: kneaddata v0.7.3

   5.4. ⚠️ IMPORTANT: Save Your Container (固化这个改变)：

    1. Open a NEW terminal window on your host machine (do not close your current Docker session).
    2. Find your current container's ID or name:
        docker ps
    (Look for the CONTAINER ID of the biobakery/workflows:latest container, e.g., 13192f2ad6e6)
    3. Commit this container to a new image named biobakery/workflows:fixed:
        docker commit 13192f2ad6e6 biobakery/workflows:fixed
        docker images
        docker ps -a
    4. From now on, whenever you want to run the workflow, use this new image name instead of :latest:
        docker run -it \
          -v /mnt/nvme4n1p1/biobakery_db:/biobakery_databases \
          -v /mnt/md1/DATA/Data_Tam_Metagenomics_2026_wastewater/X101SC25123808-Z01-J003/01.RawData/A_test_sampled:/data \
          biobakery/workflows:fixed \
          /bin/bash
    The kneaddata wrapper will already be there, and the workflow will run smoothly.

6. Rerun

   6.1. Inside the environment (SUCCESSFUL!)

    docker run -it \
        -v /mnt/nvme4n1p1/biobakery_db:/biobakery_databases \
        -v /mnt/md1/DATA/Data_Tam_Metagenomics_2026_wastewater/X101SC25123808-Z01-J003/01.RawData/A_test:/data \
        biobakery/workflows:fixed \
        /bin/bash
   
    export BIOBAKERY_WORKFLOWS_DATABASES=/biobakery_databases
    $ #OR to make this permanent, add that exact line to the ~/.bashrc file and run source ~/.bashrc.
   
    # ---- Configure databases (read-level quality control (1_KneadData), taxonomic profiling (2_MetaPhlAn), functional profiling (3_HUMAnN), and strain profiling (4_StrainPhlAn)) ----
   
    # 更新 2_MetaPhlAn_databases 路径
    python3 -c "import metaphlan, os; print(os.path.join(os.path.dirname(metaphlan.__file__), 'metaphlan_databases'))"
    /usr/local/lib/python3.6/dist-packages/metaphlan/metaphlan_databases
    ls -lh $(python3 -c "import metaphlan, os; print(os.path.join(os.path.dirname(metaphlan.__file__), 'metaphlan_databases'))")
   
    #TODO: TRY the complete metaphlan_databases from host-env to docker-system, namely from ~/mambaforge/envs/biobakery_run/lib/python3.10/site-packages/metaphlan/metaphlan_databases (v202503, 34G) to /usr/local/lib/python3.6/dist-packages/metaphlan/metaphlan_databases (v201901, 2.8G)
   
    # 更新 3_HUMAnN 配置指向该路径
    humann_config --update database_folders nucleotide /biobakery_databases/humann/chocophlan
    humann_config --update database_folders protein /biobakery_databases/humann/uniref
    humann_config --update database_folders utility_mapping /biobakery_databases/humann/utility_mapping
    humann_config --print
   
    # 1_KneadData_databases 路径: /biobakery_databases/kneaddata_db_human_genome
   
    # 4_StrainPhlAn_databases 路径: strainphlan_db_reference(empty) and strainphlan_db_markers (1.4G)
   
    # ---- If new running, optimally clean up the partial results from the failed run ----
    rm -rf /data/results/
    rm -rf /data/results/*fastqc.zip _fastqc    #IMPORTANT, so that no fastqc-related files existing under /data/results/
    rm -rf /data/results/humann
   
    $ biobakery_workflows wmgx --input /data --output /data/results
    $ biobakery_workflows wmgx_vis --input $OUTPUT_DATA --output $OUTPUT_VIS --project-name $PROJECT    #for visualizations
        * $INPUT : A directory containing shotgun sequencing data (i.e. fasta/fastq in gzipped format)
        * $OUTPUT_DATA : A directory to write the data products (i.e. abundance tables). This folder is the output folder for the first command and the input folder for the second command
        * $OUTPUT_VIS : A directory to write the visualization products (i.e. report, figures, data tables)
        * $PROJECT : The name of the project (included in the report title page)
        * Add the options --local-jobs 8 --threads 4 to run 8 local jobs at a time each with 4 threads.
        * Add the option --grid-jobs 100 to run 100 grid jobs at a time.
   
    # --qc-options="--bypass-trf" \
    # --bypass-strain-profiling
    # Run the workflow, explicitly pointing to the full databases you downloaded
    biobakery_workflows wmgx \
      --input /data \
      --output /data/results \
      --threads 64 \
      --pair-identifier "_1"
    #For A_R1.fastq.gz and A_R2.fastq.gz, the identifier is "_R1".
    #For A1a_1.fq.gz and A1a_2.fq.gz, the identifier is just "_1" (because the files end in _1 and _2 directly).
    #For A_1.fq.gz and A_2.fq.gz, the identifier is also "_1".
   
    biobakery_workflows wmgx_vis \
      --input /data/results \
      --output /data/results_vis \
      --project-name wastewater_2026_A_test_sampled
   
    #TODO_TOMORROW_2: rerun the 2 commands above again using complete metaphlan_databases: now the database under /data should copy to /usr/local/lib/python3.6/dist-packages/metaphlan/metaphlan_databases, please make backup the original simplified database as metaphlan_databases_simplified!

   6.2. (NOT_TRIED): directly run under host-environment

    docker run --rm \
      -v /mnt/nvme4n1p1/biobakery_db:/biobakery_databases \
      -v /mnt/md1/DATA/Data_Tam_Metagenomics_2026_wastewater/X101SC25123808-Z01-J003/01.RawData/A_test:/data \
      -e BIOBAKERY_WORKFLOWS_DATABASES=/biobakery_databases \
      biobakery/workflows:fixed \
      biobakery_workflows wmgx \
      -i /data \
      -o /data/results \
      --threads 32

🎉 CONGRATULATIONS! THE PIPELINE SUCCESSFULLY GENERATED ALL YOUR DATA! 🎉

Please do not worry about the Run Failed message and the traceback at the very end. Your actual scientific data is 100% complete and safe.

Here is exactly what happened and why you should be very happy with this result:

1. Why did it say “Run Failed” at the end?

If you look closely at the very first error in the log, it says:

Err: b'\n[e] The main inputs samples + references are less than 4\nThu Jun 4 11:10:02 2026: Stop StrainPhlAn 3.0 execution.\n'

StrainPhlAn is the tool that builds strain-level phylogenetic trees. To build a mathematical tree, it strictly requires a minimum of 4 samples. Because you only provided 2 samples (Sample A and Sample B), StrainPhlAn correctly realized it cannot build a tree and safely aborted.

Because the workflow manager (AnADAMA2) saw any task fail, it threw a generic RunFailed error at the very end. This is completely harmless and expected for a 2-sample run.

2. The Important Parts ALL SUCCEEDED!

Look at your logs for the actual data-generating tools:

• (Jun 04 10:46:59) ... **Completed** Task 3: kneaddata____A_R1 (Quality Control & Trimming)
• (Jun 04 10:47:43) ... **Completed** Task 8: metaphlan____A_R1 (Taxonomic Profiling)
• (Jun 04 11:07:30) ... **Completed** Task 13: humann____A_R1 (Functional/Pathway Profiling)
• (Jun 04 11:29:55) ... **Completed** Task 32: humann_merge_feature_counts

KneadData, MetaPhlAn, and HUMAnN all finished perfectly! You now have all the abundance tables you need for your analysis.

3. Where are your results?

You can find all your fully processed data in the /data/results/ directory:

• Taxonomy (Who is there?): /data/results/metaphlan/ (Look for the _taxonomic_profile.tsv files).
• Function/Pathways (What are they doing?): /data/results/humann/ (Look for the _genefamilies.tsv and _pathabundance.tsv files).
• Cleaned Reads: /data/results/kneaddata/main/ (Your high-quality, host-decontaminated .fastq files).
• Real FastQC Reports: /data/results/kneaddata/main/fastqc/ (Thanks to our hybrid wrapper, you have the real HTML reports!).

4. How to avoid this error next time

When you run the pipeline on a small number of samples (less than 4), simply add the --bypass-strain-profiling flag. This tells the workflow to skip StrainPhlAn entirely, and the pipeline will exit cleanly with a “Success” message:

biobakery_workflows wmgx \
  --input /data \
  --output /data/results \
  --threads 64 \
  --bypass-strain-profiling

5. Bonus: Generate the Visual Report!

Since your data is fully generated, you can now run the visualization workflow (just like the official documentation you found earlier) to get a beautiful HTML report:

biobakery_workflows wmgx_vis \
  --input /data/results \
  --output /data/results_vis \
  --project-name "My_Wastewater_Project"

⚠️ FINAL STEP: Save Your Masterpiece!

You have successfully conquered a highly buggy, outdated Docker image and turned it into a fully functional, production-ready bioinformatics pipeline.

Before you do anything else, open a new terminal on your host machine and save this container so you never have to do these fixes again:

# Find your container ID (e.g., 13192f2ad6e6)
docker ps

# Save it to a new permanent image
docker commit 13192f2ad6e6 biobakery/workflows:fixed

Amazing job sticking with this! Let me know if you need help interpreting the MetaPhlAn or HUMAnN output tables!

Here is the formatted table of your Docker images. (Note: The “EXTRA” column was empty for all entries, so it has been omitted to keep the table clean and readable.)