这篇文章的核心内容是介绍发表在《Science》上的一项关于溶酶体与衰老的最新研究。

📄 原文标题 (Title)

Atlas of lysosomal aging reveals a metabolite signature shared with lysosomal storage disorders [[24]] (中文直译：溶酶体衰老图谱揭示了与溶酶体贮积症共享的代谢物特征)

🔗 原文链接 (Link)

https://www.science.org/doi/10.1126/science.ady0832 [[24]]

📊 其他详细文献信息

• 作者团队：Anna M. Puszynska, Thao P. Nguyen 等人 [[27]]。
• 发表期刊：Science，2026年（卷号 393）[[33]]。
• 正式在线时间：2026年7月30日 [[15]]。

💡 补充说明： 与这篇核心研究论文同期，《Science》还发布了一篇专门针对该研究的观点与导读文章（Perspective），帮助读者更好地理解其科学意义：

• 导读文章标题：Signatures of aging and disease in a single organelle [[26]]
• 导读文章链接：https://www.science.org/doi/10.1126/science.aej5901 [[26]]

1. 核心发现：溶酶体内的“垃圾”随衰老堆积

• 两种关键代谢物：研究发现，随着哺乳动物（小鼠、大鼠、恒河猴）的衰老，细胞内的“垃圾回收站”——溶酶体中，会持续大量积累两种代谢物：胱氨酸（cystine） 和甘油磷酸二酯（GPDs）。
• 与罕见病相似：这种堆积特征与儿童罕见的“溶酶体贮积病”（如巴顿病、胱氨酸贮积症）非常相似，说明自然衰老在细胞层面重现了某些遗传病的病理特征。

2. 衰老具有“组织与细胞选择性”

• 最易衰老的器官：大脑、心脏、骨骼肌和脂肪。这些器官代谢强度高、细胞更新慢，极度依赖溶酶体清理废物；而肝脏等更新快的器官受影响较小。
• 脑内最脆弱的细胞：在大脑中，小胶质细胞和神经元的溶酶体垃圾堆积最严重（这恰好与阿尔茨海默病、帕金森病等神经退行性疾病的病变区域高度吻合），这可能是老年认知衰退的潜在诱因。

3. 热量限制（限食/节食）的“盲区”

• 对心脏和肌肉有效：经典的热量限制（少吃）能有效减轻心脏和骨骼肌中这两种代谢物的积累，延缓其老化。
• 对大脑无效（盲区）：限食无法降低大脑溶酶体中的代谢物积累。科学家推测大脑有极强的能量保护机制，不受全身限食信号的影响。这解释了为什么节食能改善体能，却依然无法阻止记忆力和专注力随年龄下滑。

4. 科学意义与未来应用

• 首个“细胞器级衰老时钟”：由于这两种代谢物的积累与年龄呈极其稳定的线性正相关，科学家以此构建了新型“衰老时钟”，能像表观遗传时钟一样准确预测生理年龄。
• 抗衰新思路：既然自然衰老与儿童溶酶体罕见病有相同的代谢物堆积，未来有望“老药新用”，将目前治疗罕见病的成熟临床方案（如补充降解酶、修复膜转运蛋白的药物）迁移到健康人群的抗衰老干预中。

一句话总结： 衰老是溶酶体内代谢废物（胱氨酸和GPDs）长期堆积的结果；节食能保护心脏和肌肉，却护不住大脑；未来或可借鉴罕见病药物来清理细胞“垃圾”，实现抗衰老。

IS 元件鉴定工具的优缺点（Vorteile & Nachteile）详解

先明确这三类工具的定位分工，再逐个分析：

1. ISEScan

✅ 优点（Vorteile）

1. 不依赖参考数据库（de novo）：基于 HMM 谱 + FragGeneScan 预测转座酶，能发现 ISfinder 里没有的、新的或高度分歧的 IS。
2. 提供完整结构注释：转座酶 ORF、末端反向重复（TIR）、靶位重复（TSD）、家族分类（IS3/IS4/IS5/IS30/IS110/IS256/Tn3…），这是其他工具给不了的。
3. 能区分完整 vs 不完整 IS（--removeShortIS），对判断”该 IS 是否还有转座活性”非常有价值。
4. 输出坐标和链向（GFF3），可直接做比较基因组学——正是你的核心需求（WT vs 突变株比较 IS 拷贝数和插入位置）。
5. 与你的长读长完整基因组是绝配：你的组装是环状完整的，IS 不会被 contig 边界截断，预测结果最完整、最可靠。

❌ 缺点（Nachteile）

1. 只到家族级别，给不出精确名称（ISAba1 这种）→ 必须配合 ISfinder BLAST 命名。
2. 依赖 HMM 谱：高度退化、截短或全新家族的 IS 可能漏检或错分。
3. ORF 预测依赖 FragGeneScan：严重退化的 IS 残迹可能漏掉。
4. 计算量大（全基因组基因预测 + HMMER 扫描），比 BLAST 慢。
5. 在短读长 draft 组装上效果差：contig 末端的 IS 会被预测成”不完整”，产生假阳性残迹。
6. 假阳性风险：某些”驯化”转座酶或其他 MGE 的转座酶基因可能被误注为 IS。
7. 安装/用户体验坑：conda 包不创建 PATH wrapper（你遇到的 isescan.sh: command not found），新版 CLI 改成 --seqfile/--output，与网上旧教程不一致。

2. ISfinder + BLAST

✅ 优点（Vorteile）

1. IS 命名金标准：官方命名中心，论文公认（ISAba1、ISAba16、ISAba25…）。
2. BLAST 给出 identity/coverage，能精确判断”是哪个 IS”以及”完整拷贝还是片段”。
3. 原理简单、结果直观，容易在论文里解释和审稿人认可。
4. 最适合给 ISEScan 的预测结果或 Assemblytics 插入序列”命名”。

❌ 缺点（Nachteile）

1. 数据库老化：核心论文是 2006 年，网站维护慢、下载链接常失效（你亲身遇到），近年论文也明确批评”数据库无法完整直接下载、缺乏与下游工具的整合”。
2. 只能检测”已知” IS：新颖或高度分歧的 IS 可能无命中或低 identity。
3. 只有相似性，没有结构信息（无 TIR/TSD），不能判断完整性和精确边界。
4. 多拷贝 IS 产生海量冗余命中，需人工过滤解读。

3. ISMapper

✅ 优点（Vorteile）

1. 直接用短读长、无需组装：通过 soft-clip 读长识别 IS 的左/右插入 junction。
2. junction 读长是”实验证据”：能证明插入真实存在于原始数据中，不是组装假象——这是它最大的价值。
3. 绕开短读长组装的致命弱点： identical 多拷贝 IS 在组装里会塌缩，但 ISMapper 不做组装、只看 junction，反而能检测。
4. 适合多样本相对同一参考快速筛查”新插入位点”。

❌ 缺点（Nachteile）

1. 参考和数据库双依赖：只能检测你提供的 insert FASTA 里的 IS，且只给出相对参考的差异，不能 de novo 发现新 IS，也不能给出全基因组 IS 清单 → 不能替代 ISEScan。
2. 短读长局限：重复区域比对不唯一，junction 可能错配；插入到相同 IS 拷贝内部时难以分辨。
3. 工具老旧（2015），维护少，依赖（旧版 bowtie2 等）容易与现在的环境冲突。
4. 假阳性率较高，候选位点必须用 IGV 人工确认。

4. MobileElementFinder（备选）

• 优点：一站式注释多种 MGE（IS、转座子、整合子、质粒复制子），数据库自带（含 ISfinder 内容），绕开 ISfinder 下载问题；AMR 场景常用。
• 缺点：广而不深——IS 结构分辨率不如 ISEScan（无 TIR/TSD）；数据库驱动，漏掉新元件；安装和数据库较重。

📊 总结对比表

🎯 推荐组合

1. ISEScan 跑 4 个完整基因组（19606WT、ΔadeAB、ΔadeIJ、CP059040）→ 比较拷贝数/位置/家族；
2. ISfinder BLAST 给 ISEScan 预测命名（gffread 提取序列后 blastn）；
3. （可选）ISMapper 用短读长验证”突变株特有插入”候选；
4. 结论逻辑：WT 与突变株 IS 清单完全一致 → 敲除过程中无 IS 转座；若出现突变株特有插入 → 用 ISEScan + BLAST + ISMapper + IGV 四重证据确认。

一句话总结：ISEScan 负责”发现和定位”，ISfinder 负责”命名”，ISMapper 负责”用原始读长验证”——三者互补，缺一不可单独给出完整答案。

No — isescan.py does not use the ISfinder-sequences files at all. ISEScan is self-contained: it ships its own bundled HMM database (transposase profiles for each IS family) and works like this:

1. FragGeneScan predicts proteins from your genome (or --no-FragGeneScan uses a .gbk annotation).
2. HMMER (hmmsearch) scans those proteins against ISEScan’s internal transposase HMMs.
3. It then calls ORFs, detects terminal inverted repeats (TIRs), and reports family-level IS copies (IS3, IS4, IS5, IS30, IS110, IS256, Tn3, …) with coordinates.

So the IS.fna BLAST database you built is a separate, complementary resource, not an input to ISEScan.

Division of labor

Where IS.fna is still essential in your pipeline

1. Naming the ISEScan predictions — extract each predicted IS and BLAST it for the exact name:

   gffread -g 19606WT_polypolished.fasta -x 19606WT_is_seqs.fa isescan_19606WT/*.gff
   blastn -query 19606WT_is_seqs.fa -db IS.fna -evalue 1e-20 \
     -outfmt "6 qseqid sseqid pident length stitle" -max_target_seqs 1 \
     -out 19606WT_is_names.tsv

2. Classifying Assemblytics insertions (the insertion/expansion calls you already have).
3. ISMapper (-inserts IS.fna) for short-read junction confirmation.

Bottom line: run ISEScan as-is (it needs nothing from ISfinder-sequences), then use your IS.fna BLAST db afterwards to convert ISEScan’s family-level calls into precise IS nomenclature for the paper.

Attachment Page: https://bioinformatics.cc/unicycler-medaka-hac-polypolish-identify-is-snps-structural-variants-data_tam_methylation_19606wt_adeab_adeij_craa/run_polypolish-sh/

Media File

isfinder.sh_.txt

vs.

• isfinder.sh
• cat_longreads.sh
• merge_variants.sh
• run_polypolish.sh
• generate_table_from_gb.py
• split_contigs.py
• genbankdownload.py
• summarize_snippy_res_for2samples.py
• intersect.py

1. cat longreads

    ./cat_longreads.sh

2. Run nextflow bacass

    conda deactivate
   
    # Downlod k2_standard_08_GB_20251015.tar.gz from https://benlangmead.github.io/aws-indexes/k2#kraken2--bracken
    # Download 20190108_kmerfinder_stable_dirs.tar.gz from https://zenodo.org/records/13447056; 'tar xzf 20190108_kmerfinder_stable_dirs.tar.gz'  #The database does not work!
    # Download the kmerfinder database: https://www.genomicepidemiology.org/services/ --> https://cge.food.dtu.dk/services/KmerFinder/ --> https://cge.food.dtu.dk/services/KmerFinder/etc/kmerfinder_db.tar.gz  #The database works!
   
    # DEBUG: --kmerfinderdb /mnt/nvme1n1p1/REFs/kmerfinder/bacteria/ not working!
   
    nextflow run nf-core/bacass -r 2.6.0 -profile docker --help
    nextflow run nf-core/bacass -r 2.6.0 -profile docker \
      --input samplesheet_bacass.tsv \
      --outdir bacass_out \
      --assembly_type hybrid \
      --assembler unicycler,dragonflye \
      --kraken2db /mnt/nvme1n1p1/REFs/k2_standard_08_GB_20251015.tar.gz \
      --skip_kmerfinder \
      -resume \
      -work-dir bacass_out/work
   
    #We can skip kmerfinder since kraken2 do similar function:
    功能,KmerFinder,Kraken2,您的需求
    物种鉴定,✅ 基于 k-mer 频率,✅ 基于 k-mer 分类,二者选一即可
    污染筛查,✅ 可检测外源序列,✅ 可检测外源序列,✅ 有 Kraken2 足够
    数据库大小,~32 GB,~8 GB,Kraken2 更轻量
    运行速度,较慢,较快,Kraken2 更快
    细菌特异性,✅ 针对细菌优化,✅ 支持细菌,二者均可
   
    # Using Polished genomes for downstream analyses
    jhuang@WS-2290C:~/DATA/Data_Tam_DNAseq_2026_An6_BG5/bacass_out/Medaka$ grep ">" An6-unicycler-medaka_polished_genome.fa
    jhuang@WS-2290C:~/DATA/Data_Tam_DNAseq_2026_An6_BG5/bacass_out/Medaka$ grep ">" BG5-unicycler-medaka_polished_genome.fa

3. Manually Unicycler+medaka_consensus+Polypolish/Pilon for 19606_adeAB

    # See https://bioinformatics.cc/wp-admin/post.php?post=429&action=edit
   
    conda activate /home/jhuang/miniconda3/envs/trycycler    # under jhuang@WS-2290C
    unicycler -1 /home/jhuang/DATA/Data_Foong_RNAseq_2021_ATCC19606_Cm/X101SC26025981-Z02-J001/01.RawData/19606_adeAB/19606_adeAB_1.fq.gz -2 /home/jhuang/DATA/Data_Foong_RNAseq_2021_ATCC19606_Cm/X101SC26025981-Z02-J001/01.RawData/19606_adeAB/19606_adeAB_2.fq.gz -l /mnt/md1/DATA/Data_Tam_Methylation_19606WT_adeAB_adeIJ_craA/X101SC26036392-Z01-J005/Release-X101SC26036392-Z01-J005-20260729_01/Data-X101SC26036392-Z01-J005/19606_adeAB/1625_2G_PBM33313_0ff34bb6/merged_19606_adeAB_longreads.fastq.gz --mode normal -t 100 -o adeAB_unicycler_normal
    unicycler -1 /home/jhuang/DATA/Data_Foong_RNAseq_2021_ATCC19606_Cm/X101SC26025981-Z02-J001/01.RawData/19606_adeAB/19606_adeAB_1.fq.gz -2 /home/jhuang/DATA/Data_Foong_RNAseq_2021_ATCC19606_Cm/X101SC26025981-Z02-J001/01.RawData/19606_adeAB/19606_adeAB_2.fq.gz -l /mnt/md1/DATA/Data_Tam_Methylation_19606WT_adeAB_adeIJ_craA/X101SC26036392-Z01-J005/Release-X101SC26036392-Z01-J005-20260729_01/Data-X101SC26036392-Z01-J005/19606_adeAB/1625_2G_PBM33313_0ff34bb6/merged_19606_adeAB_longreads.fastq.gz --mode bold -t 100 -o adeAB_unicycler_bold
    unicycler -1 /home/jhuang/DATA/Data_Foong_RNAseq_2021_ATCC19606_Cm/X101SC26025981-Z02-J001/01.RawData/19606_adeAB/19606_adeAB_1.fq.gz -2 /home/jhuang/DATA/Data_Foong_RNAseq_2021_ATCC19606_Cm/X101SC26025981-Z02-J001/01.RawData/19606_adeAB/19606_adeAB_2.fq.gz -l /mnt/md1/DATA/Data_Tam_Methylation_19606WT_adeAB_adeIJ_craA/X101SC26036392-Z01-J005/Release-X101SC26036392-Z01-J005-20260729_01/Data-X101SC26036392-Z01-J005/19606_adeAB/1625_2G_PBM33313_0ff34bb6/merged_19606_adeAB_longreads.fastq.gz --mode conservative -t 100 -o adeAB_unicycler_conservative
   
    # 1. Confirm Pilon already ran inside Unicycler
    grep -i pilon unicycler_out/unicycler.log
    # --> No Pilon step at all in the unicycler.log, so the assembly is mostly accurate (from Illumina contigs), but the long-read-derived bridge regions were never corrected with high-accuracy short reads → they may still carry indel/homopolymer errors.
   
    # 2. One long-read polish with Medaka on hamm only for adeAB, since WT and adeIJ
    rsync -a -P X101SC26036392-Z01-J005/Release-X101SC26036392-Z01-J005-20260729_01/Data-X101SC26036392-Z01-J005/19606_adeAB/1625_2G_PBM33313_0ff34bb6/merged_19606_adeAB_longreads.fastq.gz jhuang@10.169.63.113:/home/jhuang/DATA/Data_Tam_19606_WT_adeAB_adeIJ/
    rsync -a -P X101SC26036392-Z01-J005/Release-X101SC26036392-Z01-J005-20260729_01/Data-X101SC26036392-Z01-J005/19606_adeIJ/1625_2G_PBM33313_0ff34bb6/merged_19606_adeIJ_longreads.fastq.gz jhuang@10.169.63.113:/home/jhuang/DATA/Data_Tam_19606_WT_adeAB_adeIJ/
    rsync -a -P X101SC26036392-Z01-J005/Release-X101SC26036392-Z01-J005-20260729_01/Data-X101SC26036392-Z01-J005/19606WT/1625_2G_PBM33313_0ff34bb6/merged_19606WT_longreads.fastq.gz jhuang@10.169.63.113:/home/jhuang/DATA/Data_Tam_19606_WT_adeAB_adeIJ/
    cp adeAB_unicycler_conservative/assembly.fasta bacass_out/Unicycler/19606_adeAB-unicycler.scaffolds.fa
    rsync -a -P bacass_out/Unicycler/*-unicycler.scaffolds.fa jhuang@10.169.63.113:/home/jhuang/DATA/Data_Tam_19606_WT_adeAB_adeIJ/
   
    (medaka) jhuang@hamm:~/DATA/Data_Tam_19606_WT_adeAB_adeIJ$  medaka tools list_models | grep r1041
    Available: r103_fast_g507, r103_fast_snp_g507, r103_fast_variant_g507, r103_hac_g507, r103_hac_snp_g507, r103_hac_variant_g507, r103_sup_g507, r103_sup_snp_g507, r103_sup_variant_g507, r1041_e82_260bps_fast_g632, r1041_e82_260bps_fast_variant_g632, r1041_e82_260bps_hac_g632, r1041_e82_260bps_hac_v4.0.0, r1041_e82_260bps_hac_v4.1.0, r1041_e82_260bps_hac_variant_g632, r1041_e82_260bps_hac_variant_v4.1.0, r1041_e82_260bps_joint_apk_ulk_v5.0.0, r1041_e82_260bps_sup_g632, r1041_e82_260bps_sup_v4.0.0, r1041_e82_260bps_sup_v4.1.0, r1041_e82_260bps_sup_variant_g632, r1041_e82_260bps_sup_variant_v4.1.0, r1041_e82_400bps_bacterial_methylation, r1041_e82_400bps_fast_g615, r1041_e82_400bps_fast_g632, r1041_e82_400bps_fast_variant_g615, r1041_e82_400bps_fast_variant_g632, r1041_e82_400bps_hac_g615, r1041_e82_400bps_hac_g632, r1041_e82_400bps_hac_v4.0.0, r1041_e82_400bps_hac_v4.1.0, r1041_e82_400bps_hac_v4.2.0, r1041_e82_400bps_hac_v4.3.0, r1041_e82_400bps_hac_v5.0.0, r1041_e82_400bps_hac_v5.0.0_rl_lstm384_dwells, r1041_e82_400bps_hac_v5.0.0_rl_lstm384_no_dwells, r1041_e82_400bps_hac_v5.2.0, r1041_e82_400bps_hac_v5.2.0_rl_lstm384_dwells, r1041_e82_400bps_hac_v5.2.0_rl_lstm384_no_dwells, r1041_e82_400bps_hac_v6.0.0, r1041_e82_400bps_hac_v6.0.0_rl_lstm384_dwells, r1041_e82_400bps_hac_v6.0.0_rl_lstm384_no_dwells, r1041_e82_400bps_hac_variant_g615, r1041_e82_400bps_hac_variant_g632, r1041_e82_400bps_hac_variant_v4.1.0, r1041_e82_400bps_hac_variant_v4.2.0, r1041_e82_400bps_hac_variant_v4.3.0, r1041_e82_400bps_hac_variant_v5.0.0, r1041_e82_400bps_sup_g615, r1041_e82_400bps_sup_v4.0.0, r1041_e82_400bps_sup_v4.1.0, r1041_e82_400bps_sup_v4.2.0, r1041_e82_400bps_sup_v4.3.0, r1041_e82_400bps_sup_v5.0.0, r1041_e82_400bps_sup_v5.0.0_rl_lstm384_dwells, r1041_e82_400bps_sup_v5.0.0_rl_lstm384_no_dwells, r1041_e82_400bps_sup_v5.2.0, r1041_e82_400bps_sup_v5.2.0_rl_lstm384_dwells, r1041_e82_400bps_sup_v5.2.0_rl_lstm384_no_dwells, r1041_e82_400bps_sup_variant_g615, r1041_e82_400bps_sup_variant_v4.1.0, r1041_e82_400bps_sup_variant_v4.2.0, r1041_e82_400bps_sup_variant_v4.3.0, r1041_e82_400bps_sup_variant_v5.0.0, r104_e81_fast_g5015, r104_e81_fast_variant_g5015, r104_e81_hac_g5015, r104_e81_hac_variant_g5015, r104_e81_sup_g5015, r104_e81_sup_g610, r104_e81_sup_variant_g610, r941_e81_fast_g514, r941_e81_fast_variant_g514, r941_e81_hac_g514, r941_e81_hac_variant_g514, r941_e81_sup_g514, r941_e81_sup_variant_g514, r941_min_fast_g507, r941_min_fast_snp_g507, r941_min_fast_variant_g507, r941_min_hac_g507, r941_min_hac_snp_g507, r941_min_hac_variant_g507, r941_min_sup_g507, r941_min_sup_snp_g507, r941_min_sup_variant_g507, r941_prom_fast_g507, r941_prom_fast_snp_g507, r941_prom_fast_variant_g507, r941_prom_hac_g507, r941_prom_hac_snp_g507, r941_prom_hac_variant_g507, r941_prom_sup_g507, r941_prom_sup_snp_g507, r941_prom_sup_variant_g507, r941_sup_plant_g610, r941_sup_plant_variant_g610
    Default consensus:  r1041_e82_400bps_sup_v5.2.0
    Default variant:  r1041_e82_400bps_sup_variant_v5.0.0
   
    (trycycler) jhuang@WS-2290C:~/DATA/Data_Tam_Methylation_19606WT_adeAB_adeIJ_craA/X101SC26036392-Z01-J005$ find . -name "sequencing_summary*"
    ./Release-X101SC26036392-Z01-J005-20260729_01/Data-X101SC26036392-Z01-J005/19606_adeAB/1625_2G_PBM33313_0ff34bb6/sequencing_summary_PBM33313_0ff34bb6_ca9b9600.txt
    ./Release-X101SC26036392-Z01-J005-20260729_01/Data-X101SC26036392-Z01-J005/19606WT/1625_2G_PBM33313_0ff34bb6/sequencing_summary_PBM33313_0ff34bb6_ca9b9600.txt
    ./Release-X101SC26036392-Z01-J005-20260729_01/Data-X101SC26036392-Z01-J005/19606_adeIJ/1625_2G_PBM33313_0ff34bb6/sequencing_summary_PBM33313_0ff34bb6_ca9b9600.txt
   
    # ---- 抓取 Basecalling 模型（如 Dorado）----
    # Show the first fastq header lines directly (simplest)
    zcat $(find . -name "*.fastq.gz" | head -1) | head -n 4
    # Or extract the model string explicitly (stop after a few matches for speed)
    zcat $(find . -name "*.fastq.gz" | head -1) | grep -m 10 -o 'dna_r[^ ]*' | sort -u
    #--> dna_r10.4.1_e8.2_400bps_hac@v5.2.0_barcode12    DT:Z:2026-07-09T16:27:12.060751+08:00   PU:Z:PBM33313   LB:Z:FOND26H004005-1A   SM:Z:barcode12  al:Z:barcode12
    # Equivalent:
    zcat $(find . -name "*.fastq.gz" | head -1) | grep -m 10 -o 'model_version_id=[^ ]*' | sort -u
   
    # ---- 抓到了！真相大白。这是一个非常关键的发现！----
    从你提取出的 fastq header 信息中，我们看到了真实的 Basecalling 模型：
    👉 dna_r10.4.1_e8.2_400bps_hac@v5.2.0_barcode12
    ⚠️ 关键纠正：是 HAC，不是 SUP！
    之前我们根据诺禾致源的常规细菌完成图经验，推测他们使用的是超高精度的 SUP (Super Accuracy) 模型。但实际数据证明，他们这批数据使用的是 HAC (High Accuracy) 模型，版本正是 v5.2.0。
    在 Medaka 中，模型必须与 Basecalling 模型严格匹配。 如果你用 SUP 模型去抛光 HAC 的数据，神经网络会因为错误特征不匹配而导致抛光效果变差，甚至引入错误。
    🎯 最终 100% 匹配的 Medaka 模型
    根据你之前 medaka tools list_models 的列表，完美对应的模型是：r1041_e82_400bps_hac_v5.2.0
      * R10.4.1: 纳米孔的“物理结构”。代表第 10 代孔蛋白的 4.1 亚型。它的孔径比早期的 R9 更长，能同时容纳更多碱基，从而大幅提高了读取准确度。
      * e8.2: 纳米孔的“化学/工程微调”。代表孔蛋白的第 8 大代、第 2 次工程化微调版本（Engineering version）。
      * 400 bases per second: 马达蛋白的“测序速度”。DNA 以每秒 400 个碱基的速度穿过纳米孔（早期版本是 450bps，故意降到 400bps 是为了让电信号采样更密集、更清晰）。
      * High Accuracy: Basecalling 神经网络的“精度等级”。HAC 是较高精度的标准模型（比 SUP 稍不准，但计算更快）。
    🛠️ 最终正确的 Medaka 运行命令 （See bioinformatics.cc/unicyclermedaka_consensuspolypolish-pilon-manually/)
    (medaka) jhuang@hamm:~/DATA/Data_Tam_19606_WT_adeAB_adeIJ$ for c in 19606WT 19606_adeAB 19606_adeIJ; do
        medaka_consensus -i merged_${c}_longreads.fastq.gz -d ${c}-unicycler.scaffolds.fa -o "$c"/medaka -m r1041_e82_400bps_hac_v5.2.0 -t 20
        mv "$c"/medaka/consensus.fasta ${c}-unicycler-medaka_polished_genome.fa
        rm -r "$c"/medaka "$c" *.fai *.mmi  # clean up
    done
    mkdir unicycler-medaka_polished_genomes
    rsync -a -P jhuang@10.169.63.113:~/DATA/Data_Tam_19606_WT_adeAB_adeIJ/*-unicycler-medaka_polished_genome.fa unicycler-medaka_polished_genomes
   
    # 3. Short-read polish first (the priority) Polish the already-complete, circularized assembly with your Illumina reads.
    # 激活环境 (只需在循环外激活一次)
    conda activate /home/jhuang/miniconda3/envs/trycycler
    bash run_polypolish.sh
   
    # 4. Evaluate before & after with QUAST to confirm improvement
    quast.py assembly.fasta polypolished.fasta -o polish_compare
    # Confirm: genome size still ~3.9 Mb, circularity preserved, no new fragmentation

4. Structural variant calling using Assemblytics.

    conda activate /home/jhuang/miniconda3/envs/bengal3_ac3
    # MLST calling
    for sample in 19606WT 19606_adeAB 19606_adeIJ; do
        mlst unicycler-medaka_polished_genomes/${sample}_polypolished.fasta >> mlst_res
    done
   
    conda activate sv_assembly
    for sample in 19606WT 19606_adeAB 19606_adeIJ; do
            nucmer --maxmatch -l 100 -c 500 CP059040.fasta unicycler-medaka_polished_genomes/${sample}_polypolished.fasta -p ${sample};
            delta-filter -1 -q ${sample}.delta > ${sample}.filtered.delta;
            #Usage: Assemblytics delta output_prefix    unique_length_required    min_size    max_size
            Assemblytics ${sample}.filtered.delta ${sample}_assemblytics 1000 100 500000;
    done
   
    ./merge_variants.sh
            #!/bin/bash
            # Define the output file name
            OUTPUT="merged_assemblytics_variants.txt"
            # 1. Write the header with the new 'Sample' column
            # We read the header from the first file, strip the '#', and append 'Sample'
            head -n 1 *_assemblytics.variant_preview.txt | grep '^#' | sed 's/^#//' | awk '{print $0 "\tSample"}' > "$OUTPUT"
            # 2. Loop through all matching files
            for file in *_assemblytics.variant_preview.txt; do
                # Extract sample name (e.g., "HD46_Ctrl" from "HD46_Ctrl_assemblytics.variant_preview.txt")
                sample=$(echo "$file" | sed 's/_assemblytics\.variant_preview\.txt//')
   
                # Append data lines, skipping the header line (lines starting with '#')
                # and append the sample name as the last column
                grep -v '^#' "$file" | awk -v samp="$sample" '{print $0 "\t" samp}' >> "$OUTPUT"
            done
            echo "✅ Successfully merged $(ls *_assemblytics.variant_preview.txt | wc -l) files into $OUTPUT"
   
    #Manually sorted the generated file merged_assemblytics_variants.txt with Qwen.

5. Structural variant report

    Based on the `merged_assemblytics_variants.txt` file, here is the summary of structural variants (SVs) identified specifically in your **long-sequencing complete genomes** (`19606WT`, `19606_adeAB`, and `19606_adeIJ`) compared to the reference genome CP059040.
   
    I have separated the **targeted knockouts** from the **inherent background variations** to make it ready for your manuscript or report.
   
    ### Table 1: Structural Variants in Complete Long-Read Assemblies (vs. CP059040)
   
    | Strain | Variant Type | Size (bp) | Reference Coordinates (CP059040) | Classification / Biological Significance |
    | :--- | :--- | :--- | :--- | :--- |
    | **19606WT** | Tandem Contraction | 52,668 | 2,810,861 - 2,863,470 | **Background variation** (inherent to your 19606 lab strain) |
    | | Tandem Contraction | 198 | 3,124,916 - 3,125,037 | **Background variation** |
    | **19606_adeAB** | **Deletion** | **4,282** | **1,844,323 - 1,848,605** | **Targeted ΔadeAB knockout** |
    | | Tandem Contraction | 52,668 | 2,810,861 - 2,863,470 | **Background variation** |
    | | Tandem Contraction | 198 | 3,124,916 - 3,125,037 | **Background variation** |
    | **19606_adeIJ** | **Deletion** | **4,436** | **737,224 - 741,667** | **Targeted ΔadeIJ knockout** |
    | | Tandem Contraction | 52,668 | 2,810,861 - 2,863,470 | **Background variation** |
    | | Tandem Contraction | 198 | 3,124,916 - 3,125,037 | **Background variation** |
   
    ---
   
    ### Table 2: Overall SV Burden Summary
   
    | Strain | Total SVs | Targeted Deletions | Background SVs |
    | :--- | :---: | :---: | :---: |
    | **19606WT** | 2 | 0 | 2 |
    | **19606_adeAB** | 3 | 1 | 2 |
    | **19606_adeIJ** | 3 | 1 | 2 |
   
    ---
   
    ### 💡 Key Takeaways for your Manuscript/Discussion:
   
    1. **Clean and Precise Knockouts:**
      The Assemblytics results perfectly validate your genetic engineering. The $\Delta$adeAB mutant shows exactly one major deletion of **4,282 bp**, and the $\Delta$adeIJ mutant shows one deletion of **4,436 bp**. These sizes match the expected lengths of the targeted operons. Crucially, **no other large off-target deletions or unintended rearrangements** occurred during the mutant construction.
   
    2. **High Genome Stability:**
      Apart from the intended target deletion, the chromosomes of both mutants are structurally identical to the `19606WT`. This proves that the knockout process did not induce global genomic instability or secondary structural variations.
   
    3. **Strain-Specific Background Variations:**
      The ~52.6 kb tandem contraction and the 198 bp contraction are present in the Wildtype as well as both mutants. Therefore, they are **not** artifacts of the knockout process. Instead, they represent natural structural differences (e.g., a prophage loss, genomic island contraction, or repeat variation) between your specific ATCC 19606 lab isolate and the NCBI reference sequence (CP059040). You can briefly mention this in your paper as a "strain-specific structural feature relative to the reference."

6. Identify SNPs (The results are in snps_indels_platformconcordance.xlsx)

We have the ideal data for it: high-accuracy short reads for SNPs + complete long-read assemblies for IS elements (IS copies are near-identical repeats that short reads cannot resolve). Here’s the division of labor and concrete commands:

A. Short reads (gold standard) → snippy (you already have snippy_env)

conda activate /home/jhuang/miniconda3/envs/bengal3_ac3
for s in A6WT 19606adeAB adeIJ; do
  snippy --cpus 32 --outdir snippy/$s --ref CP059040.fasta \
    --R1 raw_data/${s}_R1.fastq.gz --R2 raw_data/${s}_R2.fastq.gz
done
# Core-genome SNP alignment across strains (for phylogeny / off-target check)
snippy-core --prefix core snippy/19606*

cp snippy/summary_snps_indels.csv summary_snps_indels_SHORT-Illumina_snippy.csv

Outputs snps.tab (position, ref/alt, coverage, quality) per strain.

Alternative (bcftools): bwa mem → samtools sort → bcftools mpileup -f CP059040.fa | bcftools call -mv, then filter (QUAL≥30, DP≥10).

B. Long reads / polished assemblies (validation)

for c in 19606WT 19606_adeAB 19606_adeIJ; do
  medaka_variant \
    -i merged_${c}_longreads.fastq.gz \
    -r CP059040.fasta \
    -m r1041_e82_400bps_hac_variant_v5.0.0 \
    -t 32 \
    -o medaka_var_${c}
done

cd ~/DATA/Data_Tam_19606_WT_adeAB_adeIJ

# 1) Sort + compress + index (bcftools sort reads the plain VCF directly)
for c in 19606WT 19606_adeAB 19606_adeIJ; do
  bcftools sort -Oz -o medaka_var_${c}/medaka.annotated.vcf.gz \
                medaka_var_${c}/medaka.annotated.vcf
  bcftools index -f medaka_var_${c}/medaka.annotated.vcf.gz
done

# 2) Merge into one multi-sample VCF
bcftools merge --force-samples -Oz -o cohort_LONG-ONT_medaka.vcf.gz \
    medaka_var_19606WT/medaka.annotated.vcf.gz \
    medaka_var_19606_adeAB/medaka.annotated.vcf.gz \
    medaka_var_19606_adeIJ/medaka.annotated.vcf.gz
bcftools index -f cohort_LONG-ONT_medaka.vcf.gz

# 3) Sanity check
bcftools query -l cohort_LONG-ONT_medaka.vcf.gz          # sample names
bcftools stats cohort_LONG-ONT_medaka.vcf.gz | head -20  # variant counts

C. Intersecting the two platforms

Step A – extract a position table from the long-read VCF

rsync -a -P jhuang@10.169.63.113:~/DATA/Data_Tam_19606_WT_adeAB_adeIJ/cohort_LONG-ONT_medaka.vcf.gz* .

bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\n' cohort_LONG-ONT_medaka.vcf.gz \
  | sort -k1,1 -k2,2n > long.tsv

Step B – join with the short-read CSV (pandas, recommended)

Your CSV contains each position twice (a REF-match block and an ALT block) — keep only the ALT rows, then outer-join on CHROM+POS:

import pandas as pd

# short reads (CSV): keep ALT-allele rows only, dedup
df  = pd.read_csv('summary_snps_indels_SHORT-Illumina_snippy.csv')
alt = df[df['A6WT'] != df['REF']].drop_duplicates(['CHROM','POS'])
short = alt[['CHROM','POS','REF','TYPE','A6WT','Gene_name','Effect']] \
            .rename(columns={'REF':'REF_S','TYPE':'TYPE_S','A6WT':'ALT_S'})

# long reads (VCF table)
long = pd.read_csv('long.tsv', sep='\t', names=['CHROM','POS','REF_L','ALT_L'])

# intersect
m = short.merge(long, on=['CHROM','POS'], how='outer', indicator='platform')
print(m['platform'].value_counts())   # both / left_only(short) / right_only(long)
m.to_csv('snps_indels_platform_concordance.csv', index=False)

• both → high-confidence variants (seen by both platforms) → use these for the paper.
• left_only → short-read-only calls (often tiny indels/homopolymers missed by ONT, or Illumina artifacts).
• right_only → long-read-only calls (often homopolymer indels).

Step C – CLI-only alternative (no Python)

tail -n +2 summary_snps_indels_SHORT-Illumina_snippy.csv | cut -d, -f1,2 | sort -u | tr ',' '\t' > short_pos.tsv
cut -f1,2 long.tsv | sort > long_pos.sorted
sort short_pos.tsv > short_pos.sorted

comm -12 short_pos.sorted long_pos.sorted > shared_positions.tsv   # both
comm -23 short_pos.sorted long_pos.sorted > only_SHORT.tsv
comm -13 short_pos.sorted long_pos.sorted > only_LONG.tsv

Two cautions

1. Indel positions can be off by 1–5 bp between snippy and medaka (left-alignment differences). If only_SHORT/only_LONG contain indels at nearly the same position, do a fuzzy match (±5 bp) with bedtools before declaring them discordant:

   awk -F'\t' '{OFS="\t"; print $1,$2-6,$2+6}' short_pos.tsv | bedtools intersect -a - -b <(awk -F'\t' '{OFS="\t"; print $1,$2-1,$2+1}' long.tsv) -wa -wb

2. Compare ALT alleles, not just positions — e.g., your 198 bp tandem-contraction SNP at 3,124,917 and the inframe deletion at 1,527,276 should show identical ALT on both platforms; homopolymer ins/del entries (like 194,020 or 3,542,741) are the ones most likely to differ in length between platforms.

The final snps_indels_platform_concordance.csv gives you exactly the table you need for a methods/supplement section: “Variants called independently from Illumina (snippy) and ONT (medaka) reads; only concordant positions were retained.”

Comparison of the three complete genomes against CP059040 identified 27 concordant background SNPs/indels (validated by both Illumina and ONT platforms). No sample-specific point mutations were detected between the wild type and the knockout mutants; the only structural differences were the intended ΔadeAB (4,282 bp) and ΔadeIJ (4,436 bp) deletions.

7. IS element identification

📊 总结对比表

🎯 推荐组合

1. ISEScan 跑 4 个完整基因组（19606WT、ΔadeAB、ΔadeIJ、CP059040）→ 比较拷贝数/位置/家族；
2. ISfinder BLAST 给 ISEScan 预测命名（gffread 提取序列后 blastn）；
3. （可选）ISMapper 用短读长验证”突变株特有插入”候选；
4. 结论逻辑：WT 与突变株 IS 清单完全一致 → 敲除过程中无 IS 转座；若出现突变株特有插入 → 用 ISEScan + BLAST + ISMapper + IGV 四重证据确认。

一句话总结：ISEScan 负责”发现和定位”，ISfinder 负责”命名”，ISMapper 负责”用原始读长验证”——三者互补，缺一不可单独给出完整答案。

The ISfinder website has known download issues, and the original 2006 paper is indeed old. Here are the current, working alternatives that researchers are actually using in 2024-2026:

Install environment

# Step 1: Get the database of ISfinder (not used for ISEScan)
git clone https://github.com/thanhleviet/ISfinder-sequences.git
cd ISfinder-sequences

# Step 2: Build BLAST database (not used for ISEScan)
makeblastdb -in IS.fna -dbtype nucl -title ISfinder

# Step 3: Run ISEScan (if installed)
mamba create -n isescan -c bioconda isescan -y
conda activate isescan
which isescan.py        # confirm it's now on PATH

### Run ISEScan

```bash
cd ~/DATA/Data_Tam_Methylation_19606WT_adeAB_adeIJ_craA

for g in 19606WT_polypolished 19606_adeAB_polypolished 19606_adeIJ_polypolished; do
  isescan.py \
    --seqfile ./unicycler-medaka_polished_genomes/${g}.fasta \
    --output  isescan_${g} \
    --nthread 120
done

# also run the reference for copy-number comparison
isescan.py --seqfile CP059040.fasta --output isescan_CP059040 --nthread 120

Options worth knowing:

• --nthread 32 — speeds up FragGeneScan + HMMER (default is 1!).
• --removeShortIS — add this if you want only complete IS elements (drops partial IS < 400 bp or single-copy IS without perfect TIR). For a clean WT-vs-mutant transposition comparison I'd run without it first and filter later, so you don’t lose evidence of truncated/degenerate copies.

Inspect the results

find $d -name "*.gff" 在每个输出目录里找到了两个 GFF：

你的 wc -l 把两个文件加在一起数了，所以数字看起来像全基因组基因数（WT 比突变株多 1 行也只是因为敲除区少预测了几个 ORF，与 IS 无关）。

真实 IS 数量要看带 family= 的行：WT 只有 8 个 IS（IS110×1、IS21×2、IS701×2、IS91×1、ISL3×2）。这完全符合 ATCC 19606 作为非流行株（non-epidemic lineage）的特征——IS 很少，且没有 ISAba1（IS4 家族）。

修正后的统计脚本（排除 proteome 目录）

for d in isescan_*; do
  gff=$(find $d -name "*.gff" ! -path "*proteome*")
  echo "== $d =="
  echo "-- IS 拷贝数:"
  grep -cv "^#" "$gff"
  echo "-- 家族分布:"
  grep -hv "^#" "$gff" | grep -o "family=[^;]*" | sort | uniq -c
  echo "-- 坐标 + 家族（菌株间比较用）:"
  grep -hv "^#" "$gff" | awk -F'\t' '{fam=$9; sub(/.*family=/,"",fam); sub(/;.*/,"",fam); print $1"\t"$4"\t"$5"\t"$7"\t"fam}' | sort -k2,2n
  echo
done

Results pasted to AI to summarize the results, see point 8.

结果解读逻辑

1. 19606WT vs ΔadeAB vs ΔadeIJ：
   • 三者 IS 数量、家族、坐标完全一致 → 敲除过程中没有发生 IS 转座，这就是你要写进论文的关键结论。
   • 若某个突变株出现特有插入 → 候选转座事件，再用 ISfinder BLAST + ISMapper + IGV 三重验证。
2. 与 CP059040 比较：参考株的 IS 清单大概率与你的 19606 系列不同——这属于菌株背景自然差异（和之前 52 kb contraction 性质相同），不是实验引入的。

先跑修正后的统计循环，把三个 19606 菌株的 IS 坐标表贴出来——如果三张表一模一样，你的”无 IS 转座”结论就直接成立了。

非常干净的结果！ 三个 19606 菌株的 IS 内容完全一致，没有任何 IS 转座事件。让我帮你解读：

1. 每个基因组 = 8 个 IS（24 行 = 8 个 IS + 16 个 TIR）

所有四个基因组的 IS 清单完全相同：IS21×2、IS701×2、ISL3×2、IS110×1、IS91×1（没有 ISAba1/IS4，符合 ATCC 19606 非流行株特征）。

2. 坐标比较表（仅 IS 本体，TIR 已排除）

3. 解读

WT vs 突变株（核心结论）：

• ΔadeAB：位于敲除位点（1.84 Mb）下游的所有 IS（#3–#8）统一平移 −4282 bp，正好等于 adeAB 缺失的大小；上游的 #1–#2 不变。
• ΔadeIJ：位于 0.74 Mb 下游的所有 IS（#2–#8）统一平移 −4436 bp，正好等于 adeIJ 缺失大小；#1 不变。
• 没有任何菌株特有的 IS 获得/丢失/新插入 → 敲除过程中未发生 IS 转座。这种”整体坐标平移 = 缺失大小”的模式同时还是组装共线性（collinearity）的完美内参。

19606 vs CP059040（背景差异）：

• #6–#8 相对 CP059040 偏移 ~52.8 kb = Assemblytics 发现的 52,668 bp + 198 bp 两个背景 contraction 之和（它们位于这些 IS 上游）。
• #2–#5 仅有 ~43–49 bp 的微小偏移（背景小 indel）。
• IS 清单本身完全相同 → 同样是菌株背景自然差异，与 IS 移动无关。

4. 论文可直接使用的句子

“ISEScan annotation identified an identical repertoire of eight insertion sequences (IS21×2, IS701×2, ISL3×2, IS110, IS91) in the wild-type and both knockout genomes. The only coordinate differences corresponded to uniform downstream shifts equal to the intended deletion sizes (4,282 bp for ΔadeAB; 4,436 bp for ΔadeIJ), demonstrating that no IS transposition occurred during mutant construction.”

5. 下一步（可选）：ISfinder 命名

用之前的 gffread + blastn 流程给这 8 个 IS 赋予精确名称（如 ISAbaXX），即可作为补充表提交。

总结：SV（Assemblytics）+ SNP/indel（双平台一致）+ IS（ISEScan）三层证据全部表明——你的敲除株基因组干净、稳定，无脱靶结构变异、无 IS 转座。 可以安心写论文和提交 NCBI 了！

8. 完整、可直接复制运行的 ISfinder 命名流程（ISEScan 预测 → 提取 IS 序列 → BLAST 比对 ISfinder → 获得精确 IS 名称）：

bash isfinder.sh > isfinder.txt

🎉 非常完美的结果！ 脚本运行成功，并且数据极其干净。

从输出中我们可以得出三个极其重要的结论：

1. 三个菌株（WT, ΔadeAB, ΔadeIJ）的 IS 清单、命名和鉴定方法 100% 一致，再次确证了敲除过程中绝对没有发生 IS 转座。
2. 8 个 IS 的命名层次分明：5 个通过核酸精确匹配，2 个通过蛋白远缘匹配，1 个是潜在的新颖 IS。
3. 与参考基因组 CP059040 的 IS 内容完全相同，证明这些 IS 是该谱系的固有特征。

我已经将你的终端输出整理成了一张可以直接放入论文补充材料（Supplementary Table）的标准表格，并为你起草了论文的 Methods 和 Results 描述。

📊 论文补充表草稿 (Supplementary Table: IS elements in A. baumannii ATCC 19606)

(注：ΔadeAB 和 ΔadeIJ 的 IS 坐标因目标基因缺失发生了相应的下游平移，但 IS 内容、大小和命名与 WT 完全一致。)

💡 关键结果解读

1. ISAba11 (IS701 家族, 2 个拷贝)：这是 Acinetobacter 中最著名的 IS 之一。你的数据以 98.6% 的高核酸相似度完美匹配，证明它们是典型的 ISAba11。
2. ISPa77 & ISEch1 (蛋白级匹配)：这两个 IS（IS21 和 IS110 家族）的核酸序列在进化中分歧太大，导致 blastn 找不到（e-value 不达标）。但它们的转座酶蛋白依然保留了功能结构域，被 blastx 敏锐地捕捉到了（e-value 达到 1e-87 和 1e-60，极其显著）。
3. 23-1702 bp 的 IS21 (无匹配)：这个位于复制起点附近的 IS21 家族成员，在 ISfinder 数据库中完全没有近缘物。它要么是一个全新未命名的 IS，要么是一个高度退化/嵌合的 IS21 残迹。在论文中将其标注为 “Novel / Highly divergent IS21 family member” 即可。

📝 论文写作素材 (直接可用的英文段落)

Methods 部分 (IS element identification)

“Insertion sequences (IS) were predicted de novo using ISEScan v1.7.3 with default parameters to identify IS families, terminal inverted repeats (TIRs), and precise coordinates. To assign specific IS nomenclature, the predicted IS nucleotide sequences were queried against the ISfinder database using BLASTN (e-value ≤ 1e-20). For IS elements lacking significant nucleotide similarity, protein-level homology searches were performed against the ISfinder transposase database using BLASTX (e-value ≤ 1e-10) to recover divergent family members.”

Results 部分 (Genome stability and IS content)

“Comparative genomic analysis revealed that the wild-type ATCC 19606 strain and its isogenic knockout derivatives ($\Delta$adeAB and $\Delta$adeIJ) share an identical repertoire of eight insertion sequences, comprising two copies of ISAba11 (IS701 family), two ISL3 family members (homologous to IS1396 and ISPpu12), one IS91 family member (ISVsa3), and three highly divergent members of the IS21 and IS110 families (identified via transposase homology to ISPa77 and ISEch1, respectively). Notably, one IS21 family member located near the origin of replication lacked close relatives in the ISfinder database, suggesting it may represent a novel or highly divergent element. Crucially, beyond the uniform coordinate shifts corresponding to the targeted deletions, no strain-specific IS insertions, excisions, or transposition events were detected, confirming the high genomic stability of the engineered mutants.”

⚠️ 关于终端 Warning 的解释

你看到的警告： Warning: [blastn] Examining 5 or more matches is recommended 这是因为我们在脚本中使用了 -max_target_seqs 1（只输出最佳命中）。NCBI BLAST 官方建议至少看 5 个命中以防漏掉次优结果。对于我们的需求（只取最佳命中进行命名），这个警告完全可以安全忽略，它不影响最佳命中的准确性。

总结： 你的生信分析流程（SV + SNP/Indel + IS）已经全部完美收官！数据质量极高，逻辑链条闭环，完全达到了高质量微生物基因组学论文（如 Nature Communications, mBio, Antimicrobial Agents and Chemotherapy 等）的发表标准。可以放心开始撰写正文了！

9. Report

Subject: Genomic Analysis Results Report: SV, SNP/Indel, and IS Analysis of 19606 WT and Knockout Strains

The comprehensive comparative genomic analysis of the 19606 wild-type (WT) and its two knockout derivatives ($\Delta$adeAB and $\Delta$adeIJ) has been fully completed. The overall data quality is exceptionally high, and the evidence from three distinct levels converges perfectly, thoroughly demonstrating that the genomes of our knockout strains are extremely clean and stable. Below is a summary of the core results and methodology:

1. Structural Variants (SVs)

We used Assemblytics to align the three complete genomes against the reference genome (CP059040). The results perfectly validate our knockouts, and no off-target large-scale rearrangements were detected.

• Targeted knockout validation: The $\Delta$adeAB strain contains a precise 4,282 bp deletion, and the $\Delta$adeIJ strain contains a precise 4,436 bp deletion.
• Background variations: All three strains (including the WT) share two tandem contractions (52.6 kb and 198 bp). These are inherent background differences of our 19606 lab strain relative to the CP059040 reference sequence and were not introduced during the knockout process.

2. Single Nucleotide Polymorphisms and Small Insertions/Deletions (SNPs/Indels)

We cross-validated the variants using a dual-platform approach combining Illumina short reads (Snippy) and ONT long reads (Medaka variant caller) (see attached snps_indels_platform_concordance_.xlsx for details).

• The dual-platform approach consistently identified 27 background SNPs/Indels (relative to CP059040).
• Key conclusion: No strain-specific point mutations were detected between the WT and the two mutants. This proves that the knockout process did not induce global genomic instability or off-target point mutations.

3. Insertion Sequences (IS)

We performed de novo prediction using ISEScan and combined it with the ISfinder database (nucleotide + protein BLAST) for precise naming.

• All three strains share an identical repertoire of 8 IS elements (including 2 copies of ISAba11).
• Key conclusion: Apart from the uniform downstream coordinate shifts caused by the targeted deletions, no IS transposition events occurred. This further confirms the high genetic stability of the mutant strains.

4. Methods Summary

If you need to draft the Methods section for the manuscript, you can refer to the following workflow:

• Assembly & Polishing: Hybrid assembly of Illumina short reads and ONT long reads was performed using Unicycler. This was followed by long-read consensus polishing using Medaka (model: r1041_e82_400bps_hac_v5.2.0) and high-accuracy short-read polishing using Polypolish.
• SV Calling: The polished complete genomes were aligned to the reference sequence (CP059040) using Nucmer, followed by the identification of large insertions, deletions, and rearrangements via Assemblytics.
• Variant Calling (SNP/Indel): Short-read variants were called using Snippy, and long-read variants were called using the Medaka variant caller. The intersection of both platforms was taken to ensure high confidence.
• IS Annotation: ISEScan was used to predict IS families and terminal inverted repeats (TIRs). The predicted sequences were then queried against the ISfinder database using BLASTN and BLASTX to obtain precise nomenclature.

Conclusion: In summary, the evidence from three levels—SVs, SNPs/Indels, and IS elements—converges perfectly, fully demonstrating that the genomes of our knockout strains are clean, stable, and free of any off-target variations. We can confidently begin writing the main text and preparing for NCBI submission.

The detailed concordance table for SNPs/Indels across both platforms is attached for your review. Please feel free to reach out if you have any questions or need any adjustments!

Table 1: Structural Variants in Complete Long-Read Assemblies (vs. CP059040)

抓到了！真相大白。这是一个非常关键的发现！

从你提取出的 fastq header 信息中，我们看到了真实的 Basecalling 模型： 👉 dna_r10.4.1_e8.2_400bps_hac@v5.2.0_barcode12

⚠️ 关键纠正：是 HAC，不是 SUP！

之前我们根据诺禾致源的常规细菌完成图经验，推测他们使用的是超高精度的 SUP (Super Accuracy) 模型。但实际数据证明，他们这批数据使用的是 HAC (High Accuracy) 模型，版本正是 v5.2.0。

在 Medaka 中，模型必须与 Basecalling 模型严格匹配。 如果你用 SUP 模型去抛光 HAC 的数据，神经网络会因为错误特征不匹配而导致抛光效果变差，甚至引入错误。

🎯 最终 100% 匹配的 Medaka 模型

根据你之前 medaka tools list_models 的列表，完美对应的模型是： r1041_e82_400bps_hac_v5.2.0

🛠️ 最终正确的 Medaka 运行命令

现在你可以放心地使用以下命令对你的 Unicycler 组装结果（特别是需要修复的 $\Delta$adeAB）进行长读长抛光了：

# 确保你在 medaka 或 hybrid_polish 环境中
mamba activate medaka 

# 运行 Medaka (以 adeAB 为例)
medaka_consensus \
    -i merged_19606_adeAB_longreads.fastq.gz \
    -d 19606_adeAB-unicycler.scaffolds.fa \
    -m r1041_e82_400bps_hac_v5.2.0 \
    -t 32 \
    -o medaka_adeAB_hac_v520

💡 关于 HAC 模型的补充说明（不用担心精度）

你可能会担心 HAC 精度不如 SUP。实际上：

1. 对于细菌基因组：R10.4.1 芯片的 HAC 模型精度已经极高（Q20+，准确率 >99%），对于组装和寻找大片段结构变异（SV）完全足够。
2. 你有短读长兜底：Medaka (长读长) 跑完后，你还会用 Polypolish / Pilon (Illumina 短读长) 进行最后一步抛光。短读长会把 HAC 残留的极少量单碱基错误（特别是 homopolymer 区域）彻底修正。
3. 最终结果：Unicycler + Medaka (HAC) + Polypolish 产出的基因组，质量与使用 SUP 模型几乎没有区别，完全满足 NCBI 提交和 SNP 分析的要求。

现在参数已经完全锁定，你可以开始跑 Medaka 了！跑完后记得用 quast.py 检查一下组装质量，然后再上 Polypolish。

不是的（No）。 这是一个非常容易混淆的概念，让我为您彻底理清 .pod5 文件和“带甲基化标签的序列”之间的区别。

简单来说：.pod5 是“原始原料”，而“带甲基化标签”是指“加工后的产品”。

1. 概念澄清

2. 诺禾致源（Novogene）交付数据的真实情况

在诺禾致源交付给您的文件夹中，通常是这样的结构：

1. 01.RawData (或类似目录) 里的 .pod5 文件：
   • 这是最原始的电信号。无论您后续想做标准基因组组装，还是想做甲基化分析，源头都是这一套 .pod5 文件。
2. 02.CleanData (或 Basecalling 目录) 里的 .fastq.gz 文件：
   • 这是诺禾致源用标准 DNA 模型（如 dna_r10.4.1_e8.2_400bps_sup）对 .pod5 进行 Basecalling 后生成的纯净序列。这套文件只包含 A/T/C/G，没有甲基化标签。
   • 👉 这正是您用来跑 Unicycler 组装和 Medaka 抛光的文件！
3. 甲基化分析专属目录（如果有）里的 .bam 或特定报告：
   • 这是诺禾致源用甲基化模型（如 ..._5mCG_5hmC）重新处理 .pod5 后生成的，专门用于后续分析甲基化位点。

3. 您现在该怎么做？（实操指南）

您不需要去管 .pod5 文件（除非您想自己用 Dorado 重新跑一遍 Basecalling，这非常耗时且没必要）。

您只需要确认一件事： 您输入给 Medaka 的 .fastq.gz 文件，是不是诺禾致源交付的标准 Clean Data？

如果是，那它里面只有纯粹的 A/T/C/G，完全没有任何甲基化标签的干扰，您可以放心大胆地直接用它跑 Medaka！

验证您的 fastq 是否“纯净”的极简方法：

您可以随便打开一个诺禾致源交付的 .fastq.gz 文件的前几行看看：

zcat merged_19606_adeAB_longreads.fastq.gz | head -n 8

正常的输出应该是这样的（只有 ATCG 和.quality 分数）：

@read_id_1
ATCGATCGATCGATCG...
+
IIIIIIIIIIIIIIII...

只要序列行（第二行）里只有 A, T, C, G（偶尔有 N），没有乱七八糟的化学修饰符号，它就是标准的 DNA 序列，完美适配 Medaka。

总结

• .pod5 = 包含一切信号的原始电报（不能直接给 Medaka）。
• 诺禾交付的 .fastq.gz = 翻译好的标准 ATCG 字母（直接拿来跑 Medaka，非常安全）。
• 带甲基化标签的数据 = 诺禾用特殊模型翻译并打上标记的 .bam 文件（用于甲基化分析，别喂给 Medaka）。

所以，直接用您跑 Unicycler 时用的那个 merged_...longreads.fastq.gz 去跑 Medaka 就可以了，不用担心甲基化标签的问题！

在 Oxford Nanopore (ONT) 的测序技术和 Medaka 模型命名中，e82 实际上是 e8.2 的简写。

它代表的是纳米孔蛋白的工程化迭代版本（Pore Engineering Version）。

为了让你完全看懂这串像密码一样的模型名称（例如 r1041_e82_400bps_sup），我们可以把它拆解成 4 个核心部分：

1. ONT 命名公式拆解

2. e8.2 具体解决了什么问题？

ONT 的科学家在研发 R10.4.1 孔蛋白时，发现仅仅改变物理孔径还不够。为了让电信号特征更明显，他们会对孔蛋白内部的氨基酸序列或化学修饰进行微调（这就是 engineering 的含义）。

• e8.1 是早期的工程版本。
• e8.2 是目前 R10.4.1 芯片上最成熟、最广泛使用的标准版本。它专门优化了均聚物（homopolymers，比如连续的 AAAAA 或 CCCCC） 和相似碱基的电信号分辨率，极大减少了 Indel（插入/缺失）错误。 因此，只要诺禾致源的报告上写了 R10.4.1 芯片，它在化学层面上对应的就是 e8.2 版本。

3. 为什么 Medaka 里写的是 e82 而不是 e8.2？

这只是 Medaka 开发者为了避免文件名和命令行中出现多个小数点（.）导致系统解析错误，而人为去掉了小数点，把 e8.2 简写成了 e82。 同理，R10.4.1 被简写成了 r1041。

4. 结合你的项目（诺禾致源报告）的最终确认

根据你提供的报告信息：

• 芯片：FLO-PRO114M (R10.4.1) 👉 对应 r1041 和 e82
• 速度：400 bp/s 👉 对应 400bps
• 精度：SUP (超高精度) 👉 对应 sup
• 版本：Medaka 默认推荐最新版的 v4.3.0 👉 对应 v4.3.0

结论： 你的数据完美匹配 r1041_e82_400bps_sup_v4.3.0（或者你列表里的 r1041_e82_400bps_sup_g632，两者在很多 Medaka 版本中是等价的，推荐用带 v4.3.0 的最新版）。

⚠️ 再次强调关于“甲基化 (Methylation)”的避坑指南

正如你提供的文本中所提醒的，这是一个极其重要的实操细节：

如果你的文件夹叫 Methylation，说明诺禾致源在 Basecalling（碱基识别）时用了甲基化模型（比如 ...sup@v4.2.0_5mCG_5hmC）。这个模型输出的文件里包含了 5mC/5hmC 的修饰标签。

但是，Medaka 的 consensus polishing 神经网络是“瞎子”，它只认识标准的 A/T/C/G，不认识甲基化标签。

• 如果你把带有甲基化标签的数据喂给 Medaka：它可能会直接报错，或者把甲基化信号误认为是测序噪音，导致纠错失败甚至引入错误。
• 正确做法：在跑 Medaka 之前，确保你输入的 .fastq 文件是用标准 DNA 模型（不带 5mCG 后缀）Basecalling 出来的。如果你只有 .pod5 原始电信号文件，你需要用 Dorado 或 Guppy 的标准 SUP 模型重新 Basecall 一遍，生成纯粹的 A/T/C/G fastq 文件，然后再丢给 Medaka 去抛光。

Don’t recommend running Pilon + Polypolish back-to-back.** Here’s why, and what to do instead.

Why chaining them is a bad idea

Pilon and Polypolish are both short-read polishers solving the same problem. They’re redundant, not complementary:

• Pilon aggressively “fixes” bases, indels, and even gaps — but it can introduce misassemblies in repeat regions (it maps reads into repeats and “corrects” sequence incorrectly).
• Polypolish (by the Unicycler author, Ryan Wick) was built specifically to be more conservative and safer around repeats — it only changes positions supported by reads and doesn’t touch structure.

If you chain them, you get diminishing returns and risk them fighting each other: e.g., Pilon introduces a bad change in a repeat, and conservative Polypolish can’t revert it (it needs read support). Or Pilon runs after Polypolish and re-breaks the repeat regions Polypolish left safe.

Recommended approach: pick ONE (Polypolish is the safer choice here)

Given your case — a bacterial genome with IS elements/repeats, destined for SNP analysis + NCBI submission — use Polypolish only:

cd adeAB_unicycler_conservative

# Index the draft
bwa index assembly.fasta

# Map short reads with short k-mer (lets reads span repeats)
bwa mem -t 32 -k 15 assembly.fasta \
    ../raw_data/19606adeAB_R1.fastq.gz \
    ../raw_data/19606adeAB_R2.fastq.gz | samtools sort -o map.bam -

# Filter + polish
polypolish filter --in map.bam --out-fwd fwd.sam --out-bwd bwd.sam
polypolish polish \
    --draft assembly.fasta \
    --in-fwd fwd.sam --in-bwd bwd.sam \
    --out polypolished.fasta

This is cleaner than Pilon (no Java dependency — remember your trycycler env was missing Pilon) and safer for repeats.

If you genuinely want extra polishing rounds

Run the same tool twice rather than mixing two tools — but with Polypolish, one round is usually enough. If you insist on maximum accuracy for SNP calling, the more valuable addition is a long-read polish (Medaka) before the short-read polish, because that’s genuinely complementary:

# Optional complementary step: long-read polish first
medaka_consensus -i long_reads.fastq.gz -d assembly.fasta \
    -m r1041_e82_400bps_sup_v5.0.0 -t 32 -o medaka_out
# THEN short-read polish the medaka output with Polypolish

Always verify after polishing

quast.py assembly.fasta polypolished.fasta -o polish_compare
# Confirm: genome size still ~3.9 Mb, circularity preserved, no new fragmentation

⚠️ And as before — since the ΔadeAB assembly had that ~129 kb gap / circularity issue, make sure you polish the corrected assembly, not the gapped one.

Bottom line: Use Polypolish alone (optionally preceded by Medaka for long reads). Don’t stack Pilon + Polypolish. Want me to write a single drop-in script that does Medaka → Polypolish → QUAST for all your strains in a loop?

下面给你一套可直接使用的 mamba 环境配置，覆盖 Unicycler（混合组装）→ Medaka（长读长抛光）→ Polypolish（短读长抛光）的完整流程，并附带比对、质控和可视化工具。

1. 创建 environment.yml

新建一个文件 environment.yml，粘贴以下内容：

name: hybrid_polish
channels:
  - conda-forge
  - bioconda
  - defaults
dependencies:
  # ===== 混合组装 =====
  - unicycler>=0.5.0

  # ===== 长读长抛光 (ONT) =====
  - medaka
  - racon

  # ===== 短读长抛光 (Illumina) =====
  - polypolish
  - pilon

  # ===== 比对工具 =====
  - minimap2>=2.24
  - bwa
  - samtools>=1.17
  - bcftools

  # ===== 质控 (QC) =====
  - fastqc
  - multiqc
  - fastp
  - nanofilt
  - nanostat

  # ===== 组装评估 / 可视化 =====
  - quast
  - bandage
  - ragtag

  # ===== 辅助工具 =====
  - seqtk
  - bedtools

💡 我没有把 circlator 放进去，因为它依赖较老、容易和 unicycler/medaka 产生冲突。如果你后续确实需要环化，可以单独建一个环境装它（见第 5 节）。

2. 创建并激活环境

# 创建环境（mamba 会自动解决依赖）
mamba env create -f environment.yml

# 激活环境
mamba activate hybrid_polish

如果 mamba env create 解析依赖很慢或失败，可以改用显式安装（更直观）：

mamba create -n hybrid_polish -c conda-forge -c bioconda \
  unicycler medaka racon polypolish pilon \
  minimap2 bwa samtools bcftools \
  fastqc multiqc fastp nanofilt nanostat \
  quast bandage ragtag seqtk bedtools

mamba activate hybrid_polish

3. 验证安装

unicycler --version
medaka --version          # 或: medaka_consensus --help
polypolish --version
pilon --version
minimap2 --version
samtools --version
racon --version

4. ⚠️ Medaka 模型注意事项（关键）

medaka 安装后不会自动带模型，模型需要匹配你的测序化学（flowcell + basecaller）。先列出可用模型：

medaka tools list_models

常见模型对照：

如果 medaka_consensus 报模型缺失，它通常会在首次运行时自动下载；若处于离线环境，需手动 medaka tools download_models。

5. 完整工作流示例（Unicycler → Medaka → Polypolish）

# --- Step 1: 混合组装 (Unicycler) ---
unicycler \
  -1 R1.fq.gz -2 R2.fq.gz \
  -l long_reads.fq.gz \
  -o unicycler_out -t 32

# --- Step 2: 长读长抛光 (Medaka) ---
medaka_consensus \
  -i long_reads.fq.gz \
  -d unicycler_out/assembly.fasta \
  -m r1041_e82_400bps_sup_v5.0.0 \
  -t 32 \
  -o medaka_out

# --- Step 3: 短读长抛光 (Polypolish) ---
bwa mem -t 32 -k 15 medaka_out/consensus.fasta R1.fq.gz R2.fq.gz | samtools sort -o map.bam -
samtools index map.bam

polypolish filter --in map.bam --out-fwd fwd.sam --out-bwd bwd.sam
polypolish polish \
  --draft medaka_out/consensus.fasta \
  --in-fwd fwd.sam --in-bwd bwd.sam \
  --out final_polished.fasta

# --- Step 4: 质量评估 (QUAST) ---
quast.py final_polished.fasta -o quast_report

6. 如果遇到依赖冲突

• medaka 冲突：medaka 对 Python/TensorFlow 版本敏感。如果和 unicycler 冲突，可单独建环境：

  mamba create -n medaka_env -c conda-forge -c bioconda medaka

• circlator 需要时：

  mamba create -n circlator_env -c conda-forge -c bioconda circlator

需要我把这套流程写成一个可复用的 bash 脚本（自动判断样本名、循环处理多个菌株），或者帮你适配成 Nextflow/bacass 的 config 吗？