Unicycler+medaka_consensus+Polypolish/Pilon manually

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

它代表的是纳米孔蛋白的工程化迭代版本(Pore Engineering Version)

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

1. ONT 命名公式拆解

缩写 完整含义 通俗解释(代表什么?)
r1041 R10.4.1 纳米孔的“物理结构”。代表第 10 代孔蛋白的 4.1 亚型。它的孔径比早期的 R9 更长,能同时容纳更多碱基,从而大幅提高了读取准确度。
e82 e8.2 纳米孔的“化学/工程微调”。代表孔蛋白的第 8 大代、第 2 次工程化微调版本(Engineering version)。
400bps 400 bases per second 马达蛋白的“测序速度”。DNA 以每秒 400 个碱基的速度穿过纳米孔(早期版本是 450bps,故意降到 400bps 是为了让电信号采样更密集、更清晰)。
sup Super Accuracy Basecalling 神经网络的“精度等级”。SUP 是目前最高精度的标准模型(比 HAC 更准,但计算更慢)。

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) 👉 对应 r1041e82
  • 速度: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 submissionuse 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

常见模型对照:

测序化学 推荐模型
R10.4.1 + SUP (v5) r1041_e82_400bps_sup_v5.0.0
R10.4.1 + HAC (v5) r1041_e82_400bps_hac_v5.0.0
R9.4.1 + SUP r941_prom_sup_g506
R9.4.1 + HAC r941_prom_hac_g507

如果 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 吗?

Processing Data_Patricia_Transposon_2025 v2 (Workflow for Structural Variant Calling in Nanopore Sequencing)

  1. Generate the HD46_Ctrol annotation

     mamba activate trycycler under hamm
     cd trycycler_HD46_Ctrl;
     trycycler cluster --threads 55 --assemblies assemblies/*.fasta --reads reads.fastq --out_dir trycycler;
    
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_001
     mv trycycler/cluster_001/1_contigs/J_ctg000010.fasta .
     mv trycycler/cluster_001/1_contigs/L_tig00000016.fasta .
     mv trycycler/cluster_001/1_contigs/R_tig00000001.fasta .
     mv trycycler/cluster_001/1_contigs/H_utg000001c.fasta .
    
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_002
     mv trycycler/cluster_002/1_contigs/*00000*.fasta .
     Error: unable to find a suitable common sequence
    
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_003
     mv trycycler/cluster_003/1_contigs/F_tig00000004.fasta .
     mv trycycler/cluster_003/1_contigs/L_tig00000003.fasta .
    
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_004
     mv trycycler/cluster_004/1_contigs/J_ctg000000.fasta .
     mv trycycler/cluster_004/1_contigs/P_ctg000000.fasta .
     mv trycycler/cluster_004/1_contigs/S_contig_2.fasta .
    
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_005
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_006
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_007
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_008
    
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_009
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_010
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_011
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_012
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_013
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_014
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_015
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_016
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_017
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_018
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_019
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_020
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_021
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_022
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_023
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_024
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_025
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_026
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_027
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_028
     trycycler reconcile --threads 55 --reads reads.fastq --cluster_dir trycycler/cluster_029
    
     trycycler msa --threads 55 --cluster_dir trycycler/cluster_001
     trycycler msa --threads 55 --cluster_dir trycycler/cluster_004
    
     trycycler partition --threads 55 --reads reads.fastq --cluster_dirs trycycler/cluster_001
     trycycler partition --threads 55 --reads reads.fastq --cluster_dirs trycycler/cluster_004
    
     trycycler consensus --threads 55 --cluster_dir trycycler/cluster_001
     trycycler consensus --threads 55 --cluster_dir trycycler/cluster_004
    
     #Polish --> TODO: Need to be Debugged!
     for c in trycycler/cluster_001 trycycler/cluster_004; do
         medaka_consensus -i "$c"/4_reads.fastq -d "$c"/7_final_consensus.fasta -o "$c"/medaka  -m r941_min_sup_g507 -t 12
         mv "$c"/medaka/consensus.fasta "$c"/8_medaka.fasta
         rm -r "$c"/medaka "$c"/*.fai "$c"/*.mmi  # clean up
     done
     # cat trycycler/cluster_*/8_medaka.fasta > trycycler/consensus.fasta
    
     cp trycycler/cluster_001/7_final_consensus.fasta HD46_Ctrl_chr.fasta
     cp trycycler/cluster_004/7_final_consensus.fasta HD46_Ctrl_plasmid.fasta
  2. install mambaforge https://conda-forge.org/miniforge/ (recommended)

     #download Mambaforge-24.9.2-0-Linux-x86_64.sh from website
     chmod +x Mambaforge-24.9.2-0-Linux-x86_64.sh
     ./Mambaforge-24.9.2-0-Linux-x86_64.sh
    
     To activate this environment, use:
         micromamba activate /home/jhuang/mambaforge
     Or to execute a single command in this environment, use:
         micromamba run -p /home/jhuang/mambaforge mycommand
     installation finished.
    
     Do you wish to update your shell profile to automatically initialize conda?
     This will activate conda on startup and change the command prompt when activated.
     If you'd prefer that conda's base environment not be activated on startup,
       run the following command when conda is activated:
    
     conda config --set auto_activate_base false
    
     You can undo this by running `conda init --reverse $SHELL`? [yes|no]
     [no] >>> yes
     no change     /home/jhuang/mambaforge/condabin/conda
     no change     /home/jhuang/mambaforge/bin/conda
     no change     /home/jhuang/mambaforge/bin/conda-env
     no change     /home/jhuang/mambaforge/bin/activate
     no change     /home/jhuang/mambaforge/bin/deactivate
     no change     /home/jhuang/mambaforge/etc/profile.d/conda.sh
     no change     /home/jhuang/mambaforge/etc/fish/conf.d/conda.fish
     no change     /home/jhuang/mambaforge/shell/condabin/Conda.psm1
     no change     /home/jhuang/mambaforge/shell/condabin/conda-hook.ps1
     no change     /home/jhuang/mambaforge/lib/python3.12/site-packages/xontrib/conda.xsh
     no change     /home/jhuang/mambaforge/etc/profile.d/conda.csh
     modified      /home/jhuang/.bashrc
     ==> For changes to take effect, close and re-open your current shell. <==
     no change     /home/jhuang/mambaforge/condabin/conda
     no change     /home/jhuang/mambaforge/bin/conda
     no change     /home/jhuang/mambaforge/bin/conda-env
     no change     /home/jhuang/mambaforge/bin/activate
     no change     /home/jhuang/mambaforge/bin/deactivate
     no change     /home/jhuang/mambaforge/etc/profile.d/conda.sh
     no change     /home/jhuang/mambaforge/etc/fish/conf.d/conda.fish
     no change     /home/jhuang/mambaforge/shell/condabin/Conda.psm1
     no change     /home/jhuang/mambaforge/shell/condabin/conda-hook.ps1
     no change     /home/jhuang/mambaforge/lib/python3.12/site-packages/xontrib/conda.xsh
     no change     /home/jhuang/mambaforge/etc/profile.d/conda.csh
     no change     /home/jhuang/.bashrc
     No action taken.
     WARNING conda.common.path.windows:_path_to(100): cygpath is not available, fallback to manual path conversion
     WARNING conda.common.path.windows:_path_to(100): cygpath is not available, fallback to manual path conversion
     Added mamba to /home/jhuang/.bashrc
     ==> For changes to take effect, close and re-open your current shell. <==
     Thank you for installing Mambaforge!
    
     Close your terminal window and open a new one, or run:
     #source ~/mambaforge/bin/activate
     conda --version
     mamba --version
    
     https://github.com/conda-forge/miniforge/releases
     Note
    
         * After installation, please make sure that you do not have the Anaconda default channels configured.
             conda config --show channels
             conda config --remove channels defaults
             conda config --add channels conda-forge
             conda config --show channels
             conda config --set channel_priority strict
             #conda clean --all
             conda config --remove channels biobakery
    
         * !!!!Do not install anything into the base environment as this might break your installation. See here for details.!!!!
    
     # --Deprecated method: mamba installing on conda--
     #conda install -n base --override-channels -c conda-forge mamba 'python_abi=*=*cp*'
     #    * Note that installing mamba into any other environment than base is not supported.
     #
     #conda activate base
     #conda install conda
     #conda uninstall mamba
     #conda install mamba

2: install required Tools on the mamba env

    * Sniffles2: Detect structural variants, including transposons, from long-read alignments.
    * RepeatModeler2: Identify and classify transposons de novo.
    * RepeatMasker: Annotate known transposable elements using transposon libraries.
    * SVIM: An alternative structural variant caller optimized for long-read sequencing, if needed.
    * SURVIVOR: Consolidate structural variants across samples for comparative analysis.

    mamba deactivate
    # Create a new conda environment
    mamba create -n transposon_long python=3.6 -y

    # Activate the environment
    mamba activate transposon_long

    mamba install -c bioconda sniffles
    mamba install -c bioconda repeatmodeler repeatmasker

    # configure repeatmasker database
    mamba info --envs
    cd /home/jhuang/mambaforge/envs/transposon_long/share/RepeatMasker

    #mamba install python=3.6
    mamba install -c bioconda svim
    mamba install -c bioconda survivor
  1. Test the installed tools

     # Check versions
     sniffles --version
     RepeatModeler -h
     RepeatMasker -h
     svim --help
     SURVIVOR --help
     mamba install -c conda-forge perl r
  2. Data Preparation

     Raw Signal Data: Nanopore devices generate electrical signal data as DNA passes through the nanopore.
     Basecalling: Tools like Guppy or Dorado are used to convert raw signals into nucleotide sequences (FASTQ files).
  3. Preprocessing

     Quality Filtering: Remove low-quality reads using tools like Filtlong or NanoFilt.
     Adapter Trimming: Identify and remove sequencing adapters with tools like Porechop.
  4. (Optional) Variant Calling for SNP and Indel Detection:

     Tools like Medaka, Longshot, or Nanopolish analyze the aligned reads to identify SNPs and small indels.
  5. (OFFICIAL STARTING POINT) Alignment and Structural Variant Calling: Tools such as Sniffles or SVIM detect large insertions, deletions, and other structural variants. 使用长读长测序工具如 SVIM 或 Sniffles 检测结构变异(e.g. 散在性重复序列)。

       #NOTE that the ./batch1_depth25/trycycler_WT/reads.fastq and F24A430001437_BACctmoD/BGI_result/Separate/${sample}/1.Cleandata/${sample}.filtered_reads.fq.gz are the same!
    
       # -- PREPARING the input fastq-data, merge the fastqz and move the top-directory
    
       # Under raw_data/no_sample_id/20250731_0943_MN45170_FBD12615_97f118c2/fastq_pass
       zcat ./barcode01/FBD12615_pass_barcode01_97f118c2_aa46ecf7_0.fastq.gz ./barcode01/FBD12615_pass_barcode01_97f118c2_aa46ecf7_1.fastq.gz ./barcode01/FBD12615_pass_barcode01_97f118c2_aa46ecf7_2.fastq.gz ./barcode01/FBD12615_pass_barcode01_97f118c2_aa46ecf7_3.fastq.gz ... | gzip > HD46_1.fastq.gz
       mv ./raw_data/no_sample_id/20250731_0943_MN45170_FBD12615_97f118c2/fastq_pass/HD46_1.fastq.gz ~/DATA/Data_Patricia_Transposon_2025
    
         #this are the corresponding sample names:
         #barcode 1: HD46-1
         #barcode 2: HD46-2
         #barcode 3: HD46-3
         #barcode 4: HD46-4
         mv barcode01.fastq.gz HD46_1.fastq.gz
         mv barcode02.fastq.gz HD46_2.fastq.gz
         mv barcode03.fastq.gz HD46_3.fastq.gz
         mv barcode04.fastq.gz HD46_4.fastq.gz
    
       # -- CALCULATE the coverages
         #!/bin/bash
    
         for bam in barcode*_minimap2.sorted.bam; do
             echo "Processing $bam ..."
             avg_cov=$(samtools depth -a "$bam" | awk '{sum+=$3; cnt++} END {if (cnt>0) print sum/cnt; else print 0}')
             echo -e "${bam}\t${avg_cov}" >> coverage_summary.txt
         done
    
       # ---- !!!! LOGIN the suitable environment !!!! ----
       # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #
       mamba activate transposon_long
    
       # -- TODO: AFTERNOON_DEBUG_THIS: FAILED and not_USED: Alignment and Detect structural variants in each sample using SVIM which used aligner ngmlr or mimimap2
       #mamba install -c bioconda ngmlr
       mamba install -c bioconda svim
    
       #SEARCH FOR "HD46_Ctrl_chr_plasmid.fasta" for finding the insertion-calling-commands
       # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! for all 4 options #
       # ---- Option_1: minimap2 (aligner) + SVIM (structural variant caller) --> SUCCESSFUL ----
    
       for sample in HD46_1 HD46_2 HD46_3 HD46_4 HD46_5 HD46_6 HD46_7 HD46_8 HD46_13; do
           #INS,INV,DUP:TANDEM,DUP:INT,BND
           svim reads --aligner minimap2 --nanopore minimap2+svim_${sample}    ${sample}.fastq.gz HD46_Ctrl_chr_plasmid.fasta  --cores 20 --types INS --min_sv_size 100 --sequence_allele --insertion_sequences --read_names;
       done
    
       #svim alignment svim_alignment_minmap2_1_re 1.sorted.bam CP020463_.fasta --types INS --sequence_alleles --insertion_sequences --read_names
    
       # ---- Option_2: minamap2 (aligner) + Sniffles2 (structural variant caller) --> SUCCESSFUL ----
       #Minimap2: A commonly used aligner for nanopore sequencing data.
       #    Align Long Reads to the WT Reference using Minimap2
       #sniffles -m WT.sorted.bam -v WT.vcf -s 10 -l 50 -t 60
       #  -s 20: Requires at least 20 reads to support an SV for reporting. --> 10
       #  -l 50: Reports SVs that are at least 50 base pairs long.
       #  -t 60: Uses 60 threads for faster processing.
       for sample in HD46_1 HD46_2 HD46_3 HD46_4 HD46_5 HD46_6 HD46_7 HD46_8 HD46_13; do
           #minimap2 --MD -t 60 -ax map-ont HD46_Ctrl_chr_plasmid.fasta ./batch1_depth25/trycycler_${sample}/reads.fastq | samtools sort -o ${sample}.sorted.bam
           minimap2 --MD -t 60 -ax map-ont HD46_Ctrl_chr_plasmid.fasta ${sample}.fastq.gz | samtools sort -o ${sample}_minimap2.sorted.bam
           samtools index ${sample}_minimap2.sorted.bam
           sniffles -m ${sample}_minimap2.sorted.bam -v ${sample}_minimap2+sniffles.vcf -s 10 -l 50 -t 60
           #QUAL < 20 ||
           bcftools filter -e "INFO/SVTYPE != 'INS'" ${sample}_minimap2+sniffles.vcf > ${sample}_minimap2+sniffles_filtered.vcf
       done
    
         #Estimating parameter...
         #        Max dist between aln events: 44
         #        Max diff in window: 76
         #        Min score ratio: 2
         #        Avg DEL ratio: 0.0112045
         #        Avg INS ratio: 0.0364027
         #Start parsing... CP020463
         #                # Processed reads: 10000
         #                # Processed reads: 20000
         #        Finalizing  ..
         #Start genotype calling:
         #        Reopening Bam file for parsing coverage
         #        Finalizing  ..
         #Estimating parameter...
         #        Max dist between aln events: 28
         #        Max diff in window: 89
         #        Min score ratio: 2
         #        Avg DEL ratio: 0.013754
         #        Avg INS ratio: 0.17393
         #Start parsing... CP020463
         #                # Processed reads: 10000
         #                # Processed reads: 20000
         #                # Processed reads: 30000
         #                # Processed reads: 40000
    
         # Results:
         # * barcode01_minimap2+sniffles.vcf
         # * barcode01_minimap2+sniffles_filtered.vcf
         # * barcode02_minimap2+sniffles.vcf
         # * barcode02_minimap2+sniffles_filtered.vcf
         # * barcode03_minimap2+sniffles.vcf
         # * barcode03_minimap2+sniffles_filtered.vcf
         # * barcode04_minimap2+sniffles.vcf
         # * barcode04_minimap2+sniffles_filtered.vcf
    
       #ERROR: No MD string detected! Check bam file! Otherwise generate using e.g. samtools. --> No results!
       #for sample in barcode01 barcode02 barcode03 barcode04; do
       #    sniffles -m svim_reads_minimap2_${sample}/${sample}.fastq.minimap2.coordsorted.bam -v sniffles_minimap2_${sample}.vcf -s 10 -l 50 -t 60
       #    bcftools filter -e "INFO/SVTYPE != 'INS'" sniffles_minimap2_${sample}.vcf > sniffles_minimap2_${sample}_filtered.vcf
       #done
    
       # ---- Option_3: NGMLR (aligner) + SVIM (structural variant caller) --> SUCCESSFUL ----
       for sample in HD46_1 HD46_2 HD46_3 HD46_4 HD46_5 HD46_6 HD46_7 HD46_8 HD46_13; do
           svim reads --aligner ngmlr --nanopore    ngmlr+svim_${sample}       ${sample}.fastq.gz HD46_Ctrl_chr_plasmid.fasta  --cores 10;
       done
    
       # ---- Option_4: NGMLR (aligner) + sniffles (structural variant caller) --> SUCCESSFUL ----
       for sample in HD46_1 HD46_2 HD46_3 HD46_4 HD46_5 HD46_6 HD46_7 HD46_8 HD46_13; do
           sniffles -m ngmlr+svim_${sample}/${sample}.fastq.ngmlr.coordsorted.bam -v ${sample}_ngmlr+sniffles.vcf -s 10 -l 50 -t 60
           bcftools filter -e "INFO/SVTYPE != 'INS'" ${sample}_ngmlr+sniffles.vcf > ${sample}_ngmlr+sniffles_filtered.vcf
       done
    
       #END
  6. Compare and integrate all results produced by minimap2+sniffles and ngmlr+sniffles, and check them each position in IGV!

     # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #
     mv HD46_1_minimap2+sniffles_filtered.vcf    HD46-1_minimap2+sniffles_filtered.vcf
     mv HD46_1_ngmlr+sniffles_filtered.vcf       HD46-1_ngmlr+sniffles_filtered.vcf
     mv HD46_2_minimap2+sniffles_filtered.vcf    HD46-2_minimap2+sniffles_filtered.vcf
     mv HD46_2_ngmlr+sniffles_filtered.vcf       HD46-2_ngmlr+sniffles_filtered.vcf
     mv HD46_3_minimap2+sniffles_filtered.vcf    HD46-3_minimap2+sniffles_filtered.vcf
     mv HD46_3_ngmlr+sniffles_filtered.vcf       HD46-3_ngmlr+sniffles_filtered.vcf
     mv HD46_4_minimap2+sniffles_filtered.vcf    HD46-4_minimap2+sniffles_filtered.vcf
     mv HD46_4_ngmlr+sniffles_filtered.vcf       HD46-4_ngmlr+sniffles_filtered.vcf
     mv HD46_5_minimap2+sniffles_filtered.vcf    HD46-5_minimap2+sniffles_filtered.vcf
     mv HD46_5_ngmlr+sniffles_filtered.vcf       HD46-5_ngmlr+sniffles_filtered.vcf
     mv HD46_6_minimap2+sniffles_filtered.vcf    HD46-6_minimap2+sniffles_filtered.vcf
     mv HD46_6_ngmlr+sniffles_filtered.vcf       HD46-6_ngmlr+sniffles_filtered.vcf
     mv HD46_7_minimap2+sniffles_filtered.vcf    HD46-7_minimap2+sniffles_filtered.vcf
     mv HD46_7_ngmlr+sniffles_filtered.vcf       HD46-7_ngmlr+sniffles_filtered.vcf
     mv HD46_8_minimap2+sniffles_filtered.vcf    HD46-8_minimap2+sniffles_filtered.vcf
     mv HD46_8_ngmlr+sniffles_filtered.vcf       HD46-8_ngmlr+sniffles_filtered.vcf
     mv HD46_13_minimap2+sniffles_filtered.vcf   HD46-13_minimap2+sniffles_filtered.vcf
     mv HD46_13_ngmlr+sniffles_filtered.vcf      HD46-13_ngmlr+sniffles_filtered.vcf
  7. (NOT_USED) Filtering low-complexity insertions using RepeatMasker (TODO: how to use RepeatModeler to generate own lib?)

       python vcf_to_fasta.py variants.vcf variants.fasta
       #python filter_low_complexity.py variants.fasta filtered_variants.fasta retained_variants.fasta
       #Using RepeatMasker to filter the low-complexity fasta, the used h5 lib is
       /home/jhuang/mambaforge/envs/transposon_long/share/RepeatMasker/Libraries/Dfam.h5    #1.9G
       python /home/jhuang/mambaforge/envs/transposon_long/share/RepeatMasker/famdb.py -i /home/jhuang/mambaforge/envs/transposon_long/share/RepeatMasker/Libraries/Dfam.h5 names 'bacteria' | head
       Exact Matches
       =============
       2 bacteria (blast name), Bacteria 
    (scientific name), eubacteria (genbank common name), Monera (in-part), Procaryotae (in-part), Prokaryota (in-part), Prokaryotae (in-part), prokaryote (in-part), prokaryotes (in-part) Non-exact Matches ================= 1783272 Terrabacteria group (scientific name) 91061 Bacilli (scientific name), Bacilli Ludwig et al. 2010 (authority), Bacillus/Lactobacillus/Streptococcus group (synonym), Firmibacteria (synonym), Firmibacteria Murray 1988 (authority) 1239 Bacillaeota (synonym), Bacillaeota Oren et al. 2015 (authority), Bacillota (synonym), Bacillus/Clostridium group (synonym), clostridial firmicutes (synonym), Clostridium group firmicutes (synonym), Firmacutes (synonym), firmicutes (blast name), Firmicutes (scientific name), Firmicutes corrig. Gibbons and Murray 1978 (authority), Low G+C firmicutes (synonym), low G+C Gram-positive bacteria (common name), low GC Gram+ (common name) Summary of Classes within Firmicutes: * Bacilli (includes many common pathogenic and non-pathogenic Gram-positive bacteria, taxid=91061) * Bacillus (e.g., Bacillus subtilis, Bacillus anthracis) * Staphylococcus (e.g., Staphylococcus aureus, Staphylococcus epidermidis) * Streptococcus (e.g., Streptococcus pneumoniae, Streptococcus pyogenes) * Listeria (e.g., Listeria monocytogenes) * Clostridia (includes many anaerobic species like Clostridium and Clostridioides) * Erysipelotrichia (intestinal bacteria, some pathogenic) * Tissierellia (less-studied, veterinary relevance) * Mollicutes (cell wall-less, includes Mycoplasma species) * Negativicutes (includes some Gram-negative, anaerobic species) RepeatMasker -species Bacilli -pa 4 -xsmall variants.fasta python extract_unmasked_seq.py variants.fasta.masked unmasked_variants.fasta #bcftools filter -i ‘QUAL>30 && INFO/SVLEN>100’ variants.vcf -o filtered.vcf # #bcftools view -i ‘SVTYPE=”INS”‘ variants.vcf | bcftools query -f ‘%CHROM\t%POS\t%REF\t%ALT\t%INFO\n’ > insertions.txt #mamba install -c bioconda vcf2fasta #vcf2fasta variants.vcf -o insertions.fasta #grep “SEQS” variants.vcf | awk ‘{ print $1, $2, $4, $5, $8 }’ > insertions.txt #python3 filtering_low_complexity.py # #vcftools –vcf input.vcf –recode –out filtered_output –minSVLEN 100 #bcftools filter -e ‘INFO/SEQS ~ “^(G+|C+|T+|A+){4,}”‘ variants.vcf -o filtered.vcf # — calculate the percentage of reads To calculate the percentage of reads that contain the insertion from the VCF entry, use the INFO and FORMAT fields provided in the VCF record. Step 1: Extract Relevant Information In the provided VCF entry: RE (Reads Evidence): 733 – the total number of reads supporting the insertion. GT (Genotype): 1/1 – this indicates a homozygous insertion, meaning all reads covering this region are expected to have the insertion. AF (Allele Frequency): 1 – a 100% allele frequency, indicating that every read in this sample supports the insertion. DR (Depth Reference): 0 – the number of reads supporting the reference allele. DV (Depth Variant): 733 – the number of reads supporting the variant allele (insertion). Step 2: Calculate Percentage of Reads Supporting the Insertion Using the formula: Percentage of reads with insertion=(DVDR+DV)×100 Percentage of reads with insertion=(DR+DVDV​)×100 Substitute the values: Percentage=(7330+733)×100=100% Percentage=(0+733733​)×100=100% Conclusion Based on the VCF record, 100% of the reads support the insertion, indicating that the insertion is fully present in the sample (homozygous insertion). This is consistent with the AF=1 and GT=1/1 fields. * In your VCF file generated by Sniffles, the REF=N in the results has a specific meaning: * In a standard VCF, the REF field usually contains the reference base(s) at the variant position. * For structural variants (SVs), especially insertions, there is no reference sequence replaced; the insertion occurs between reference bases. * Therefore, Sniffles uses N as a placeholder in the REF field to indicate “no reference base replaced”. * The actual inserted sequence is then stored in the ALT field.
  8. Why some records have UNRESOLVED in the FILTER field in the Excel output.

     1. Understanding the format
    
         The data appears to be structural variant (SV) calls from Sniffles, probably in a VCF-like tabular format exported to Excel:
    
             * gi|1176884116|gb|CP020463.1| → reference sequence
             * Positions: 1855752 and 2422820
             * N → insertion event
             * SVLEN=999 → size of the insertion
             * AF → allele frequency
             * GT:DR:DV → genotype, depth reference, depth variant (1/1:0:678, example values for a PASS variant)
             * FILTER → whether the variant passed filters (UNRESOLVED means it didn’t pass)
    
     2. What UNRESOLVED usually means
    
         In Sniffles:
    
         * UNRESOLVED is assigned to SVs when the tool cannot confidently resolve the exact sequence or breakpoint.
         * Reasons include:
             - Low read support (RE, DV) relative to the expected coverage
             - Ambiguous alignment at repetitive regions
             - Conflicting strand or orientation signals
             - Allele frequency inconsistent with expectations
    
     3. Examine your two records
    
         First record
    
             POS: 1855752
             SVTYPE: INS
             SVLEN: 999
             RE: 68
             AF: 1
             GT: 1/1
             FILTER: UNRESOLVED
    
         Observations:
    
         * AF = 1 → allele frequency 100%, homozygous insertion
         * RE = 68 → 68 reads support the variant, decent coverage
         * Still UNRESOLVED → likely because Sniffles could not resolve the inserted sequence precisely; sometimes long insertions in repetitive regions are hard to reconstruct fully even with good read support.
    
         Second record
    
             POS: 2422820
             SVTYPE: INS
             SVLEN: 999
             RE: 22
             AF: 0.025522
             GT: 0/0
             FILTER: UNRESOLVED
    
         Observations:
    
         * AF = 0.0255 → very low allele frequency (~2.5%)
         * RE = 22, DR = 840 → very low variant reads vs reference
         * GT = 0/0 → homozygous reference
         * Sniffles marks it UNRESOLVED because the variant is essentially noise, not confidently detected.
    
     4. Key difference between the two
         Feature First record    Second record
         Allele frequency (AF)   1 (high)    0.0255 (very low)
         Variant reads (RE)  68  22
         Genotype (GT)   1/1 0/0
         Reason for UNRESOLVED   Unresolvable inserted sequence
    
     ✅ 5. Conclusion
    
         * Sniffles marks a variant as UNRESOLVED when the SV cannot be confidently characterized.
         * Even if there is good read support (first record), complex insertions can’t always be reconstructed fully.
         * Very low allele frequency (second record) also triggers UNRESOLVED because the signal is too weak compared to background noise.
         * Essentially: “UNRESOLVED” ≠ bad data, it’s just unresolved uncertainty.
  9. (NOT_SURE_HOW_TO_USE) Polishing of assembly: Use tools like Medaka to refine variant calls by leveraging consensus sequences derived from nanopore data.

       mamba install -c bioconda medaka
       medaka-consensus -i aligned_reads.bam -r reference.fasta -o polished_output -t 4
  10. Compare Insertions Across Samples

     Merge Variants Across Samples: Use SURVIVOR to merge and compare the detected insertions in all samples against the WT:
    
     SURVIVOR merge input_vcfs.txt 1000 1 1 1 0 30 merged.vcf
    
         Input: List of VCF files from Sniffles2.
         Output: A consolidated VCF file with shared and unique variants.
    
     Filter WT Insertions:
    
         Identify transposons present only in samples 1–9 by subtracting WT variants using bcftools:
    
             bcftools isec WT.vcf merged.vcf -p comparison_results
  11. Validate and Visualize

     Visualize with IGV: Use IGV to inspect insertion sites in the alignment and confirm quality.
    
     igv.sh
    
     Validate Findings:
         Perform PCR or additional sequencing for key transposon insertion sites to confirm results.
  12. Alternatives to TEPID for Long-Read Data

     If you’re looking for transposon-specific tools for long reads:
    
         REPET: A robust transposon annotation tool compatible with assembled genomes.
         EDTA (Extensive de novo TE Annotator):
             A pipeline to identify, classify, and annotate transposons.
             Works directly on your assembled genomes.
    
             perl EDTA.pl --genome WT.fasta --type all
  13. The WT.vcf file in the pipeline is generated by detecting structural variants (SVs) in the wild-type (WT) genome aligned against itself or using it as a baseline reference. Here’s how you can generate the WT.vcf:

     Steps to Generate WT.vcf
     1. Align WT Reads to the WT Reference Genome
    
     The goal here is to create an alignment of the WT sequencing data to the WT reference genome to detect any self-contained structural variations, such as native insertions, deletions, or duplications.
    
     Command using Minimap2:
    
     minimap2 -ax map-ont WT.fasta WT_reads.fastq | samtools sort -o WT.sorted.bam
    
     Index the BAM file:
    
     samtools index WT.sorted.bam
    
     2. Detect Structural Variants with Sniffles2
    
     Run Sniffles2 on the WT alignment to call structural variants:
    
     sniffles --input WT.sorted.bam --vcf WT.vcf
    
     This step identifies:
    
         Native transposons and insertions present in the WT genome.
         Other structural variants that are part of the reference genome or sequencing artifacts.
    
     Key parameters to consider:
    
         --min_support: Adjust based on your WT sequencing coverage.
         --max_distance: Define proximity for merging variants.
         --min_length: Set a minimum SV size (e.g., >50 bp for transposons).
  14. Clean and Filter the WT.vcf, Variant Filtering: Remove low-confidence variants based on read depth, quality scores, or allele frequency.

     To ensure the WT.vcf only includes relevant transposons or SVs:
    
         Use bcftools or similar tools to filter out low-confidence variants:
    
         bcftools filter -e "QUAL < 20 || INFO/SVTYPE != 'INS'" WT.vcf > WT_filtered.vcf
         bcftools filter -e "QUAL < 1 || INFO/SVTYPE != 'INS'" 1_.vcf > 1_filtered_.vcf
  15. NOTE that in this pipeline, the WT.fasta (reference genome) is typically a high-quality genome sequence from a database or a well-annotated version of your species’ genome. It is not assembled from the WT.fastq sequencing reads in this context. Here’s why:

     Why Use a Reference Genome (WT.fasta) from a Database?
    
         Higher Quality and Completeness:
             Database references (e.g., NCBI, Ensembl) are typically well-assembled, highly polished, and annotated. They serve as a reliable baseline for variant detection.
    
         Consistency:
             Using a standard reference ensures consistent comparisons across your WT and samples (1–9). Variants detected will be relative to this reference, not influenced by possible assembly errors.
    
         Saves Time:
             Assembling a reference genome from WT reads requires significant computational effort. Using an existing reference streamlines the analysis.
    
     Alternative: Assembling WT from FASTQ
    
     If you don’t have a high-quality reference genome (WT.fasta) and must rely on your WT FASTQ reads:
    
         Assemble the genome from your WT.fastq:
             Use long-read assemblers like Flye, Canu, or Shasta to create a draft genome.
    
         flye --nano-raw WT.fastq --out-dir WT_assembly --genome-size 
    Polish the assembly using tools like Racon (with the same reads) or Medaka for higher accuracy. Use the assembled and polished genome as your WT.fasta reference for further steps. Key Takeaways: If you have access to a reliable, high-quality reference genome, use it as the WT.fasta. Only assemble WT.fasta from raw reads (WT.fastq) if no database reference is available for your organism.
  16. Annotate Transposable Elements: Tools like ANNOVAR or SnpEff provide functional insights into the detected variants.

     # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #
     #Using snpEff to annotate the insertion!
     conda activate /home/jhuang/miniconda3/envs/spandx
# --> BUG:
LOCUS       HD46_Ctrl 2707468 bp    DNA     circular BCT
    02-OCT-2025
DEFINITION  Staphylococcus epidermidis strain HD46-ctrl chromosome, whole
            genome shotgun sequence.
ACCESSION
VERSION

# --> DEBUG: adapt the genbank-file header as follows:
LOCUS       HD46_Ctrl 2707468 bp    DNA     circular BCT 02-OCT-2025
DEFINITION  Staphylococcus epidermidis strain HD46-ctrl chromosome, whole
            genome shotgun sequence.
ACCESSION   HD46_Ctrl
VERSION     HD46_Ctrl.1
DBLINK      BioProject: PRJNA1337321
            BioSample: SAMN52215988
KEYWORDS    .
SOURCE      Staphylococcus epidermidis
  ORGANISM  Staphylococcus epidermidis
            Bacteria; Firmicutes; Bacilli; Bacillales; Staphylococcaceae;
            Staphylococcus.
COMMENT     Annotated genome for HD46_Ctrl.
...
    # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #
    mkdir ~/miniconda3/envs/spandx/share/snpeff-5.1-2/data/HD46_Ctrl
    cp HD46_Ctrl_chr.gb ~/miniconda3/envs/spandx/share/snpeff-5.1-2/data/HD46_Ctrl/genes.gbk

    vim ~/miniconda3/envs/spandx/share/snpeff-5.1-2/snpEff.config  #HD46_Ctrl.genome : HD46_Ctrl
    /home/jhuang/miniconda3/envs/spandx/bin/snpEff build -genbank HD46_Ctrl      -d

    sed -i 's/^cluster_001_consensus/HD46_Ctrl.1/' HD46-8_ngmlr+sniffles_filtered.vcf
    sed -i 's/^cluster_001_consensus/HD46_Ctrl.1/' HD46-13_ngmlr+sniffles_filtered.vcf
    #snpEff eff -nodownload -no-downstream -no-intergenic -ud 100 -v HD46_Ctrl HD46-8_ngmlr+sniffles_filtered.vcf > HD46-8_ngmlr+sniffles_filtered.annotated.vcf
    #snpEff eff -nodownload -no-downstream -no-intergenic -ud 100 -v HD46_Ctrl HD46-13_ngmlr+sniffles_filtered.vcf > HD46-13_ngmlr+sniffles_filtered.annotated.vcf

    # HD46-8
    snpEff ann -Xmx8g -v -hgvs -canon -ud 200 \
    -stats HD46-8_snpeff_stats.html \
    HD46_Ctrl \
    HD46-8_ngmlr+sniffles_filtered.vcf \
    > HD46-8_ngmlr+sniffles_filtered.annotated.vcf

    # HD46-13
    snpEff ann -Xmx8g -v -hgvs -canon -ud 200 \
    -stats HD46-13_snpeff_stats.html \
    HD46_Ctrl \
    HD46-13_ngmlr+sniffles_filtered.vcf \
    > HD46-13_ngmlr+sniffles_filtered.annotated.vcf
  1. Summarize the results as a Excel-file

    # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #
    conda activate plot-numpy1
    #python generate_common_vcf.py
    #mv common_variants.xlsx putative_transposons.xlsx
    
    # * Reads each of your VCFs.
    # * Filters variants → only keep those with FILTER == PASS.
    # * Compares the two aligner methods (minimap2+sniffles2 vs ngmlr+sniffles2) per sample.
    # * Keeps only variants that appear in both methods for the same sample.
    # * Outputs: An Excel file with the common variants and a log text file listing which variants were filtered out, and why (not_PASS or not_COMMON_in_two_VCF).
    
    #python generate_fuzzy_common_vcf_v1.py
    #Sample PASS_minimap2   PASS_ngmlr  COMMON
    #  HD46-Ctrl_Ctrl   39  29  28
    #  HD46-1   39  32  29
    #  HD46-2   40  32  28
    #  HD46-3   38  30  27
    #  HD46-4   46  35  32
    #  HD46-5   40  35  31
    #  HD46-6   43  35  30
    #  HD46-7   40  33  28
    #  HD46-8   37  20  11
    #  HD46-13  39  38  27
    
    #Sample PASS_minimap2   PASS_ngmlr  COMMON_FINAL
    #HD46-Ctrl_Ctrl 39  29  6
    #HD46-1 39  32  8
    #HD46-2 40  32  8
    #HD46-3 38  30  6
    #HD46-4 46  35  8
    #HD46-5 40  35  9
    #HD46-6 43  35  10
    #HD46-7 40  33  8
    #HD46-8 37  20  4
    #HD46-13    39  38  5
    
    # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #
    #!!!! Summarize the results of ngmlr+sniffles !!!!
    python merge_ngmlr+sniffles_filtered_results_and_summarize.py
    
    #!!!! Post-Processing !!!!
    #DELETE "2186168    N   

    . PASS” in Sheet HD46-13 and Summary #DELETE “2427785 N CGTCAGAATCGCTGTCTGCGTCCGAGTCACTGTCTGAGTCTGAATCACTATCTGCGTCTGAGTCACTGTCTG . PASS” due to “0/1:169:117” in HD46-13 and Summary #DELETE “2441640 N GCTCATTAAGAATCATTAAATTAC . PASS” due to 0/1:170:152 in HD46-13 and Summary

  2. Source code of merge_ngmlr+sniffles_filtered_results_and_summarize.py

    python add_ann_to_excel.py         --excel merged_ngmlr+sniffles_variants.xlsx         --sheet8 "HD46-8"         --sheet13 "HD46-13"         --vcf8 HD46-8_ngmlr+sniffles_filtered.annotated.vcf         --vcf13 HD46-13_ngmlr+sniffles_filtered.annotated.vcf         --out merged_ngmlr+sniffles_variants_with_ANN.xlsx
    
    #!/usr/bin/env python3
    # -*- coding: utf-8 -*-
    """
    Add SnpEff ANN columns (for SVTYPE=INS) from annotated VCFs into an Excel workbook,
    with detailed debug about why CHROM+POS may not match.
    
    Key improvements:
    - Stronger CHROM/POS normalization (strip 'chr', unify MT naming, coerce numbers).
    - Explicit detection and logging of sheet key columns used.
    - Debug block prints:
    * Unique key counts in sheet vs VCF (before/after normalization)
    * Example non-matching keys from the sheet and from the VCF (top N)
    * Chromosome naming diagnostics (e.g., 'chr' presence, 'MT'/'M' harmonization)
    * Off-by-N diagnostics via --pos_tolerance (counts for would-match @ ±N)
    * Optional preview of the sheet's SV type column, if present
    - Safer ANN parsing and aggregation.
    - Command-line options: --debug_examples, --pos_tolerance
    
    ANN filling is done only for **exact** (CHROM, POS) equality (as before).
    Tolerance is used only for *diagnostics*, not for filling, to avoid incorrect merges.
    """
    
    import argparse
    import gzip
    import io
    import re
    from pathlib import Path
    from typing import List, Tuple, Dict, Iterable, Set
    
    import pandas as pd
    
    FALLBACK_ANN_FIELDS: List[str] = [
        'Allele','Annotation','Annotation_Impact','Gene_Name','Gene_ID',
        'Feature_Type','Feature_ID','Transcript_BioType','Rank','HGVS.c',
        'HGVS.p','cDNA.pos/cDNA.length','CDS.pos/CDS.length','AA.pos/AA.length',
        'Distance','Errors_Warnings_Info'
    ]
    
    def open_text_maybe_gzip(path: Path):
        if str(path).endswith('.gz'):
            return io.TextIOWrapper(gzip.open(path, 'rb'), encoding='utf-8', errors='ignore')
        return open(path, 'r', encoding='utf-8', errors='ignore')
    
    def normalize_chrom(col: pd.Series) -> pd.Series:
        s = col.astype(str).str.strip()
        s = s.str.replace(r'^(chr|CHR)', '', regex=True)
        # Standardize mitochondrial names to "MT"
        s = s.str.replace(r'^(M|MtDNA|MTDNA|Mito|Mitochondrion)$', 'MT', regex=True, case=False)
        return s.str.upper()
    
    def normalize_pos(col: pd.Series) -> pd.Series:
        # Excel can make ints look like floats; coerce then Int64
        # (We keep Int64 nullable for robustness; we never compare NaNs.)
        s = pd.to_numeric(col, errors='coerce')
        # If people had 0-based starts in the sheet (rare for INS), this won't fix it,
        # but the tolerance debug will reveal a +1 shift if present.
        return s.astype('Int64')
    
    def parse_vcf_ann(vcf_path: Path) -> Tuple[pd.DataFrame, List[str]]:
        ann_fields = None
        header_cols = None
        records = []
    
        with open_text_maybe_gzip(vcf_path) as f:
            for line in f:
                if line.startswith('##INFO=<ID=ANN'):
                    m = re.search(r'Format:\s*([^">]+)', line)
                    if m:
                        ann_fields = [s.strip() for s in m.group(1).split('|')]
                if line.startswith('#CHROM'):
                    header_cols = line.strip().lstrip('#').split('\t')
                    break
    
            if not header_cols:
                raise RuntimeError(f"Could not find VCF header line (#CHROM ...) in {vcf_path}")
    
            if not ann_fields:
                ann_fields = FALLBACK_ANN_FIELDS
    
            ann_cols = [f'ANN_{x}' for x in ann_fields]
    
            for line in f:
                if not line or line[0] == '#':
                    continue
                parts = line.rstrip('\n').split('\t')
                if len(parts) < len(header_cols):
                    continue
                row = dict(zip(header_cols, parts))
                info = row.get('INFO', '')
    
                # Only INS
                if not re.search(r'(?:^|;)SVTYPE=INS(?:;|$)', info):
                    continue
    
                chrom = row.get('#CHROM') or row.get('CHROM')
                pos_str = row.get('POS')
                try:
                    pos = int(pos_str)
                except Exception:
                    continue
    
                # Extract ANN entries
                ann_match = re.search(r'(?:^|;)ANN=([^;]+)', info)
                ann_entries = ann_match.group(1).split(',') if ann_match else []
    
                field_values: Dict[str, List[str]] = {k: [] for k in ann_fields}
                for ann in ann_entries:
                    items = ann.split('|')
                    if len(items) < len(ann_fields):
                        items += [''] * (len(ann_fields) - len(items))
                    elif len(items) > len(ann_fields):
                        items = items[:len(ann_fields)]
                    for k, v in zip(ann_fields, items):
                        field_values[k].append(v)
    
                joined = {f'ANN_{k}': (';'.join(v) if v else '') for k, v in field_values.items()}
                records.append({'CHROM': chrom, 'POS': pos, **joined})
    
        df = pd.DataFrame.from_records(records)
        if not df.empty:
            df['POS'] = pd.to_numeric(df['POS'], errors='coerce').astype('Int64')
            df['CHROM'] = normalize_chrom(df['CHROM'])
        return df, ann_cols
    
    def detect_key_columns(df: pd.DataFrame) -> Dict[str, str]:
        chrom_candidates = ['CHROM', '#CHROM', 'Chrom', 'Chromosome', 'chrom', 'chr', 'Chr']
        pos_candidates   = ['POS', 'Position', 'position', 'pos', 'Start', 'start']
        mapping = {}
        for c in chrom_candidates:
            if c in df.columns:
                mapping['CHROM'] = c
                break
        for p in pos_candidates:
            if p in df.columns:
                mapping['POS'] = p
                break
        return mapping
    
    def normalize_chrom_pos_df(df: pd.DataFrame, keys: Dict[str, str]) -> pd.DataFrame:
        out = df.copy()
        out[keys['CHROM']] = normalize_chrom(out[keys['CHROM']])
        out[keys['POS']]   = normalize_pos(out[keys['POS']])
        return out
    
    def summarize_chr_formats(series: pd.Series, label: str):
        raw = series.astype(str)
        has_chr_prefix = raw.str.startswith(('chr','CHR')).sum()
        mt_like = raw.str.fullmatch(r'(M|MtDNA|MTDNA|Mito|Mitochondrion)', case=False).sum()
        print(f"[{label}] CHROM diagnostics:")
        print(f"  total rows: {len(raw)}")
        print(f"  with 'chr'/'CHR' prefix: {has_chr_prefix}")
        print(f"  mitochondrial names like M/MtDNA/etc: {mt_like}")
    
    def keys_set(df: pd.DataFrame, chrom_col: str, pos_col: str) -> Set[Tuple[str, int]]:
        # Drop NA POS, NA CHROM
        sub = df[[chrom_col, pos_col]].dropna()
        # Ensure ints (drop NA after coercion)
        sub = sub[(sub[pos_col].astype('Int64').notna())]
        return set(zip(sub[chrom_col].astype(str), sub[pos_col].astype('int64')))
    
    def tolerance_match_count(sheet_keys: Iterable[Tuple[str,int]],
                            vcf_keys: Set[Tuple[str,int]],
                            tol: int) -> int:
        if tol <= 0:
            return sum(1 for k in sheet_keys if k in vcf_keys)
        cnt = 0
        for chrom, pos in sheet_keys:
            if (chrom, pos) in vcf_keys:
                cnt += 1
            else:
                matched = False
                # check +/- 1..tol
                for d in range(1, tol+1):
                    if (chrom, pos - d) in vcf_keys or (chrom, pos + d) in vcf_keys:
                        matched = True
                        break
                if matched:
                    cnt += 1
        return cnt
    
    def debug_match_report(df_sheet: pd.DataFrame,
                        vcf_df: pd.DataFrame,
                        keys: Dict[str, str],
                        debug_examples: int = 15,
                        pos_tolerance: int = 1):
        print("\n=== DEBUG: Matching overview ===")
        # Raw diagnostics
        summarize_chr_formats(df_sheet[keys['CHROM']], label="SHEET (raw)")
        summarize_chr_formats(vcf_df['CHROM'], label="VCF (normalized)")
    
        # Normalize sheet
        df_norm = normalize_chrom_pos_df(df_sheet, keys)
        print(f"Detected key columns -> CHROM: '{keys['CHROM']}'  POS: '{keys['POS']}'")
        # Basic stats
        n_sheet_all = len(df_sheet)
        n_sheet_key_nonnull = df_norm[keys['CHROM']].notna().sum() - df_norm[keys['CHROM']].isna().sum()
        n_sheet_pos_nonnull = df_norm[keys['POS']].notna().sum()
        print(f"SHEET rows total: {n_sheet_all}")
        print(f"SHEET rows with non-null CHROM: {n_sheet_key_nonnull}, non-null POS: {n_sheet_pos_nonnull}")
    
        # Unique key counts
        sheet_norm_keys_df = df_norm.rename(columns={keys['CHROM']: 'CHROM', keys['POS']: 'POS'})
        sheet_norm_keys_df = sheet_norm_keys_df.dropna(subset=['CHROM','POS'])
        sheet_norm_keys_df['POS'] = sheet_norm_keys_df['POS'].astype('Int64')
        sheet_keys_unique = keys_set(sheet_norm_keys_df, 'CHROM', 'POS')
        vcf_keys_unique   = keys_set(vcf_df, 'CHROM', 'POS')
    
        print(f"Unique (CHROM,POS) keys -> SHEET: {len(sheet_keys_unique)}  VCF(INS): {len(vcf_keys_unique)}")
    
        # Exact match count
        exact_matches = len(sheet_keys_unique & vcf_keys_unique)
        print(f"Exact key matches (SHEET∩VCF): {exact_matches}")
    
        # Tolerance diagnostics (diagnose off-by-one etc.)
        if pos_tolerance > 0:
            approx_matches = tolerance_match_count(sheet_keys_unique, vcf_keys_unique, pos_tolerance)
            print(f"Keys that would match within ±{pos_tolerance}: {approx_matches}")
    
        # Show some examples of non-matching keys from SHEET
        if debug_examples > 0:
            not_in_vcf = sorted(k for k in sheet_keys_unique if k not in vcf_keys_unique)
            not_in_sheet = sorted(k for k in vcf_keys_unique if k not in sheet_keys_unique)
            print(f"\nExamples of SHEET keys not found in VCF (showing up to {debug_examples}):")
            for k in not_in_vcf[:debug_examples]:
                print("  SHEET-only:", k)
            print(f"\nExamples of VCF keys not found in SHEET (showing up to {debug_examples}):")
            for k in not_in_sheet[:debug_examples]:
                print("  VCF-only:", k)
    
        # Try to detect a type column and report counts
        type_cols = [c for c in df_sheet.columns if c.lower() in ('svtype','type','variant_type','sv_type')]
        if type_cols:
            tcol = type_cols[0]
            is_ins = df_sheet[tcol].astype(str).str.upper() == 'INS'
            print(f"\nType column detected: '{tcol}'. SHEET rows with INS: {int(is_ins.sum())} / {len(df_sheet)}")
            # Of the INS rows, how many have keys that match?
            ins_keys = keys_set(df_norm[is_ins], keys['CHROM'], keys['POS'])
            exact_ins_matches = len(ins_keys & vcf_keys_unique)
            print(f"  INS-only exact key matches: {exact_ins_matches} / {len(ins_keys)}")
            if pos_tolerance > 0:
                approx_ins_matches = tolerance_match_count(ins_keys, vcf_keys_unique, pos_tolerance)
                print(f"  INS-only matches within ±{pos_tolerance}: {approx_ins_matches} / {len(ins_keys)}")
        else:
            print("\nNo explicit type column found in SHEET.")
    
    def merge_ann_into_sheet(df_sheet: pd.DataFrame, vcf_df: pd.DataFrame, ann_cols: List[str],
                            pos_tolerance: int = 1, debug_examples: int = 15) -> pd.DataFrame:
        df = df_sheet.copy()
    
        keys = detect_key_columns(df)
        if 'CHROM' not in keys or 'POS' not in keys:
            print("WARNING: Could not detect CHROM/POS columns in sheet; ANN columns will be empty.")
            for c in ann_cols:
                if c not in df.columns:
                    df[c] = ''
            return df
    
        # DEBUG: run a comprehensive match report
        debug_match_report(df, vcf_df, keys, debug_examples=debug_examples, pos_tolerance=pos_tolerance)
    
        # Normalize sheet keys for merge
        df_norm = normalize_chrom_pos_df(df, keys)
    
        # Prepare VCF map (unique by CHROM,POS), aggregate ANN fields
        vcf_use = vcf_df.copy()
        if vcf_use.empty:
            print("NOTE: No INS records found in VCF; ANN columns will be created but empty.")
        else:
            agg = {c: lambda s: ';'.join([x for x in s.astype(str).tolist() if x]) for c in ann_cols}
            vcf_use = vcf_use.groupby(['CHROM', 'POS'], as_index=False).agg(agg)
    
        # Identify potential type column in sheet
        type_cols = [c for c in df.columns if c.lower() in ('svtype','type','variant_type','sv_type')]
        has_type = bool(type_cols)
        if has_type:
            tcol = type_cols[0]
            is_ins = df[tcol].astype(str).str.upper() == 'INS'
            print(f"\nMERGE: using type column '{tcol}' -> rows marked INS: {int(is_ins.sum())} / {len(df)}")
        else:
            is_ins = pd.Series([False]*len(df), index=df.index)
            print("\nMERGE: no type column -> will fill ANN wherever exact (CHROM,POS) matches VCF INS.")
    
        # Left merge on exact keys only (do not use tolerance for filling, just for diagnostics)
        left = df_norm.rename(columns={keys['CHROM']: 'CHROM', keys['POS']: 'POS'})
        merged = left.merge(vcf_use[['CHROM','POS'] + ann_cols], on=['CHROM','POS'], how='left', suffixes=('',''))
    
        # Initialize ANN columns on original df
        for c in ann_cols:
            if c not in df.columns:
                df[c] = ''
    
        # Fill values:
        for c in ann_cols:
            values = merged[c]
            if has_type:
                df.loc[is_ins, c] = values[is_ins].fillna('').astype(str).values
            else:
                df[c] = values.fillna('').astype(str).values
    
        # Report matching stats on the actual merge
        matched = merged[ann_cols].notna().any(axis=1).sum()
        print(f"\nMERGE RESULT: rows with any ANN filled (exact VCF match): {int(matched)} / {len(df)}")
    
        # Additional hint if tolerance suggests many near-misses
        if pos_tolerance > 0:
            sheet_keys = keys_set(left, 'CHROM', 'POS')
            vcf_keys_unique = keys_set(vcf_use, 'CHROM', 'POS')
            approx = tolerance_match_count(sheet_keys, vcf_keys_unique, pos_tolerance)
            if approx > matched:
                print(f"NOTE: There appear to be {approx - matched} additional rows that would match within ±{pos_tolerance}.")
                print("      This often indicates a 0-based vs 1-based position shift or use of END instead of POS in the sheet.")
    
        return df
    
    def main():
        ap = argparse.ArgumentParser()
        ap.add_argument('--excel', default='merged_ngmlr+sniffles_variants.xlsx', help='Input Excel workbook')
        ap.add_argument('--sheet8', default='HD46-8', help='Sheet name for HD46-8 sample')
        ap.add_argument('--sheet13', default='HD46-13', help='Sheet name for HD46-13 sample')
        ap.add_argument('--vcf8', default='HD46-8_ngmlr+sniffles_filtered.annotated.vcf', help='Annotated VCF for HD46-8')
        ap.add_argument('--vcf13', default='HD46-13_ngmlr+sniffles_filtered.annotated.vcf', help='Annotated VCF for HD46-13')
        ap.add_argument('--out', default='merged_ngmlr+sniffles_variants_with_ANN.xlsx', help='Output Excel path')
        ap.add_argument('--debug_examples', type=int, default=15, help='How many non-match examples to print from each side')
        ap.add_argument('--pos_tolerance', type=int, default=1, help='Diagnostic tolerance (±N bp) for off-by-N checks (used for debug only)')
        args = ap.parse_args()
    
        excel_path = Path(args.excel)
        vcf8_path = Path(args.vcf8)
        vcf13_path = Path(args.vcf13)
        out_path = Path(args.out)
    
        # Load sheets (resolve case-insensitive names)
        xls = pd.ExcelFile(excel_path)
        def resolve_sheet(name: str) -> str:
            if name in xls.sheet_names:
                return name
            lower_map = {s.lower(): s for s in xls.sheet_names}
            return lower_map.get(name.lower(), name)
    
        sheet8 = resolve_sheet(args.sheet8)
        sheet13 = resolve_sheet(args.sheet13)
    
        df8 = pd.read_excel(excel_path, sheet_name=sheet8)
        df13 = pd.read_excel(excel_path, sheet_name=sheet13)
    
        # Parse VCFs (INS only)
        vcf8_df, ann_cols = parse_vcf_ann(vcf8_path)
        vcf13_df, _ = parse_vcf_ann(vcf13_path)
    
        print(f"VCF8 INS variants: {len(vcf8_df)}; VCF13 INS variants: {len(vcf13_df)}")
        print(f"ANN subfields ({len(ann_cols)}): {', '.join(ann_cols)}")
    
        # Merge with diagnostics
        df8_out = merge_ann_into_sheet(df8, vcf8_df, ann_cols,
                                    pos_tolerance=args.pos_tolerance,
                                    debug_examples=args.debug_examples)
        df13_out = merge_ann_into_sheet(df13, vcf13_df, ann_cols,
                                        pos_tolerance=args.pos_tolerance,
                                        debug_examples=args.debug_examples)
    
        # Save
        with pd.ExcelWriter(out_path, engine='xlsxwriter') as writer:
            df8_out.to_excel(writer, sheet_name=sheet8, index=False)
            df13_out.to_excel(writer, sheet_name=sheet13, index=False)
    
        print(f"\nDone. Wrote: {out_path.resolve()}")
    
    if __name__ == '__main__':
        main()
  3. Manually merge all contents of ANN=? to a seperate column ‘ANN’ in the isolate-specific sheets in the Excel-file.

    #Add CHROM and HD46_Ctrl.1 to first column of the input Excel-file
    (plot-numpy1) jhuang@WS-2290C:~/DATA/Data_Patricia_Transposon_2025$ python add_ann_to_excel.py         --excel merged_ngmlr+sniffles_variants.xlsx         --sheet8 "HD46-8"         --sheet13 "HD46-13"         --vcf8 HD46-8_ngmlr+sniffles_filtered.annotated.vcf         --vcf13 HD46-13_ngmlr+sniffles_filtered.annotated.vcf         --out merged_ngmlr+sniffles_variants_with_ANN.xlsx
    #DEL some columns (INFO, NN_Allele, ANN_Rank, ANN_Errors_Warnings_Info, from the table, and COPY the summary-sheet to the final table.
  4. Run nextflow bacass

    # -- samplesheet_bacass.tsv --
    #ID R1  R2  LongFastQ   Fast5   GenomeSize
    #HD46_Ctrl          HD46_Ctrl.fastq.gz  NA  NA
    #HD46_1         HD46_1.fastq.gz NA  NA
    #An6    /mnt/md1/DATA/Data_Tam_DNAseq_2026_An6_BG5/X101SC26036392-Z01-J002/clean_data/An6/An6_L1_1.clean.rd.fq.gz   /mnt/md1/DATA/Data_Tam_DNAseq_2026_An6_BG5/X101SC26036392-Z01-J002/clean_data/An6/An6_L1_2.clean.rd.fq.gz   /mnt/md1/DATA/Data_Tam_DNAseq_2026_An6_BG5/X101SC26036392-Z01-J003/Release-X101SC26036392-Z01-J003-20260513_01/Data-X101SC26036392-Z01-J003/An6/2157_4C_PBK79106_7ec05c46/merged_An6_longreads.fastq.gz NA  2.7m
    #BG5    /mnt/md1/DATA/Data_Tam_DNAseq_2026_An6_BG5/X101SC26036392-Z01-J002/clean_data/BG5/BG5_L1_1.clean.rd.fq.gz   /mnt/md1/DATA/Data_Tam_DNAseq_2026_An6_BG5/X101SC26036392-Z01-J002/clean_data/BG5/BG5_L1_2.clean.rd.fq.gz   /mnt/md1/DATA/Data_Tam_DNAseq_2026_An6_BG5/X101SC26036392-Z01-J003/Release-X101SC26036392-Z01-J003-20260513_01/Data-X101SC26036392-Z01-J003/BG5/2157_4C_PBK79106_7ec05c46/merged_BG5_longreads.fastq.gz NA  6.5m
    
    conda deactivate
    # DEBUG: --kmerfinderdb /mnt/nvme1n1p1/REFs/kmerfinder/bacteria/ not working, maybe due to the version, since 2.5.0 was working (see below)!
    #nextflow run nf-core/bacass -r 2.5.0 -profile docker \
    #--input samplesheet.tsv \
    #--outdir bacass_out \
    #--assembly_type long \
    #--kraken2db /mnt/nvme1n1p1/REFs/k2_standard_08_GB_20251015.tar.gz \
    #--kmerfinderdb /mnt/nvme1n1p1/REFs/kmerfinder/bacteria/ \
    #-resume
    
    #For hybrid-assembly: --assembly_type hybrid --assembler unicycler,dragonflye [unicycler,autocycler,canu,dragonflye,flye,miniasm,raven,megahit] \
    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 long \
      --assembler unicycler,dragonflye \
      --kraken2db /mnt/nvme1n1p1/REFs/k2_standard_08_GB_20251015.tar.gz \
      --skip_kmerfinder \
      -resume \
      -c unicycler.config \
      -work-dir bacass_out/work
    
    #SAVE bacass_out/Kmerfinder/kmerfinder_summary.csv to bacass_out/Kmerfinder/An6?/An6?_kmerfinder_results.xlsx
    
    #busco example results:
    Input_file      Dataset Complete        Single  Duplicated      Fragmented      Missing n_markers       Scaffold N50    Contigs N50     Percent gaps    Number of scaffolds
    wt_cef.scaffolds.fa     bacteria_odb10  98.4    98.4    0.0     1.6     0.0     124     285852  285852  0.000%  45
    wt_cipro.scaffolds.fa   bacteria_odb10  90.3    89.5    0.8     8.1     1.6     124     7434    7434    0.000%  1699
  5. Detecting the next closest genome

    mamba activate gtdbtk
    
    # 验证环境变量是否加载成功
    echo $GTDBTK_DATA_PATH
    # 应输出:/mnt/nvme4n1p1/gtdb_data/release232
    
    # 3. 运行分类(你提供的命令 + 实用参数)
    gtdbtk classify_wf \
      --genome_dir ./bacass_out/Medaka \
      --out_dir gtdb_out \
      --cpus 64 \
      --extension .fa \
      --prefix mygenome
    
    # 4. 查看结果
    cat gtdb_out/classify/mygenome.bac120.summary.tsv   # 细菌结果
  6. Structural variant calling

    conda activate sv_assembly
    
    # MLST calling
    for sample in HD46_Ctrl HD46_1 HD46_2 HD46_3 HD46_4 HD46_5 HD46_6 HD46_7 HD46_8 HD46_13 _WT _1 _2 _3 _4 _5 _7 _8 _9 _10; do
        mlst bacass_out/Medaka/${sample}-unicycler-medaka_polished_genome.fa >> mlst_res
    done
    
    # After running MLST and genome taxonomy checks, I found that CP020463 is only a suitable reference for the first dataset (_WT, _1–_10), because all of these samples share ST 86 — the same sequence type as CP020463. For the second dataset (HD46 series), CP020463 is not an appropriate reference. MLST and GTDB-Tk classification results (see mygenome.bac120.summary2.xlsx and mlst_res.xlsx) show that these genomes are genetically distinct.
    for sample in _WT _1 _2 _3 _4 _5 _7 _8 _9 _10; do
            nucmer --maxmatch -l 100 -c 500 CP020463.fasta bacass_out/Medaka/${sample}-dragonflye-medaka_polished_genome.fa -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
    samtools faidx bacass_out/Medaka/HD46_Ctrl-dragonflye-medaka_polished_genome.fa contig00001 > bacass_out/Medaka/HD46_Ctrl_chrom.fa
    for sample in HD46_1 HD46_2 HD46_5 HD46_6 HD46_7; do
            nucmer --maxmatch -l 100 -c 500 bacass_out/Medaka/HD46_Ctrl_chrom.fa bacass_out/Medaka/${sample}-dragonflye-medaka_polished_genome.fa -p ${sample};
            delta-filter -1 -q ${sample}.delta > ${sample}.filtered.delta;
            #Usage: Assemblytics delta output_prefix    unique_length_required    min_size    max_size; Note that we use a large threshold 500,000 nt.
            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 into two parts: one is HD46-series and one is _*series.
    
    unicycler -l HD46_Ctrl.fastq.gz --mode normal -t 40 -o HD46_Ctrl_unicycler_normal
    unicycler -l HD46_3.fastq.gz --mode normal -t 40 -o HD46_3_unicycler_normal
    unicycler -l HD46_4.fastq.gz --mode normal -t 40 -o HD46_4_unicycler_normal
    unicycler -l HD46_13.fastq.gz --mode normal -t 40 -o HD46_13_unicycler_normal

Processing Data_Ute_RNAseq_FINAL/Data_RNA-Seq_MKL-1+WaGa

Here is the optimized version. Key improvements: (1) it verifies that the two Geneid/gene_name columns from the pasted files are identical and keeps only one copy; (2) it never blindly overwrites column names — R’s automatic header sanitization (check.names=TRUE) produces exactly the names you were assigning by hand (e.g. 042_MKL-1_wt_EVX042_MKL.1_wt_EV, MKL-1_EV-RNA_87MKL.1_EV.RNA_87), so we use those; (3) the final friendly rename is done via an explicit name→name map instead of position-based colnames(x) <- c(...), so a sample can never be silently mislabeled — any mismatch throws an error.

# ---------------------------------------------------------------
# 0) Read the merged count table
#    - Do NOT set row.names=1 yet: we first need to compare the two
#      Geneid / gene_name columns that come from the two pasted files.
#    - check.names=TRUE (default) sanitizes the header deterministically:
#        042_MKL-1_wt_EV        -> X042_MKL.1_wt_EV
#        MKL-1_EV-RNA_87        -> MKL.1_EV.RNA_87
#        2nd Geneid / gene_name -> Geneid.1 / gene_name.1
#    (i.e. exactly the names you previously assigned by hand)
# ---------------------------------------------------------------
d.full <- read.delim2("merged_gene_counts_40samples.txt", sep="\t", header=TRUE)

stopifnot(ncol(d.full) == 44)   # 2 x (Geneid + gene_name) + 40 samples
colnames(d.full)                # eyeball-check the auto-generated names

# ---------------------------------------------------------------
# 1) Check that the two Geneid / gene_name columns are identical
#    If this fails, the two files were pasted in different gene order
#    -> stop and fix (diagnosis: which(d.full$Geneid != d.full$Geneid.1))
# ---------------------------------------------------------------
stopifnot(identical(d.full$Geneid,    d.full$Geneid.1))
stopifnot(identical(d.full$gene_name, d.full$gene_name.1))

# identical -> keep only one copy of each
d.full$Geneid.1    <- NULL
d.full$gene_name.1 <- NULL

# Ensembl IDs as row names; gene_name was also dropped before DESeq2 previously
stopifnot(!anyDuplicated(d.full$Geneid))
rownames(d.full) <- d.full$Geneid
d.full$Geneid    <- NULL
d.full$gene_name <- NULL

stopifnot(ncol(d.full) == 40)   # only sample count columns remain

# ---------------------------------------------------------------
# 2) Reorder columns by biological group, using the AUTO-GENERATED names
#    (order matches the condition/donor/batch vectors of the old script)
# ---------------------------------------------------------------
col_order <- c(
  # MKL-1 RNA
  "MKL.1_RNA","MKL.1_RNA_118","MKL.1_RNA_147",
  # MKL-1 wt EV
  "MKL.1_EV.RNA","MKL.1_EV.RNA_2","MKL.1_EV.RNA_118","MKL.1_EV.RNA_87","MKL.1_EV.RNA_27",
  "X042_MKL.1_wt_EV",
  # MKL-1 EV DMSO / Dox
  "X042_MKL.1_sT_DMSO","X0505_MKL.1_sT_DMSO_EV",
  "X042_MKL.1_scr_DMSO_EV","X0505_MKL.1_scr_DMSO_EV",
  "X042_MKL.1_sT_Dox","X0505_MKL.1_sT_Dox_EV",
  "X042_MKL.1_scr_Dox_EV","X0505_MKL.1_scr_Dox_EV",
  # WaGa RNA
  "WaGa_RNA","WaGa_RNA_118","WaGa_RNA_147",
  # WaGa wt EV
  "WaGa_EV.RNA","WaGa_EV.RNA_2","WaGa_EV.RNA_118","WaGa_EV.RNA_147","WaGa_EV.RNA_226",
  "X1107_WaGa_wt_EV","X1605_WaGa_wt_EV","X2706_WaGa_wt_EV",
  # WaGa EV DMSO / Dox
  "X1107_WaGa_sT_DMSO_EV","X1605_WaGa_sT_DMSO_EV","X2706_WaGa_sT_DMSO_EV",
  "X1107_WaGa_scr_DMSO_EV","X1605_WaGa_scr_DMSO_EV","X2706_WaGa_scr_DMSO_EV",
  "X1107_WaGa_sT_Dox_EV","X1605_WaGa_sT_Dox_EV","X2706_WaGa_sT_Dox_EV",
  "X1107_WaGa_scr_Dox_EV","X1605_WaGa_scr_Dox_EV","X2706_WaGa_scr_Dox_EV")

stopifnot(length(col_order) == 40,
          !anyDuplicated(col_order),
          all(col_order %in% colnames(d.full)))   # errors instead of mislabeling

reordered.raw <- d.full[, col_order]

# ---------------------------------------------------------------
# 3) Friendly sample names via an EXPLICIT map (position-independent).
#    Kept because all downstream code refers to e.g. "MKL-1 EV sT DMSO 042".
# ---------------------------------------------------------------
name_map <- c(
  "MKL.1_RNA"               = "MKL-1 RNA",
  "MKL.1_RNA_118"           = "MKL-1 RNA 118",
  "MKL.1_RNA_147"           = "MKL-1 RNA 147",
  "MKL.1_EV.RNA"            = "MKL-1 EV",
  "MKL.1_EV.RNA_2"          = "MKL-1 EV 2",
  "MKL.1_EV.RNA_118"        = "MKL-1 EV 118",
  "MKL.1_EV.RNA_87"         = "MKL-1 EV 87",
  "MKL.1_EV.RNA_27"         = "MKL-1 EV 27",
  "X042_MKL.1_wt_EV"        = "MKL-1 EV 042",
  "X042_MKL.1_sT_DMSO"      = "MKL-1 EV sT DMSO 042",
  "X0505_MKL.1_sT_DMSO_EV"  = "MKL-1 EV sT DMSO 0505",
  "X042_MKL.1_scr_DMSO_EV"  = "MKL-1 EV scr DMSO 042",
  "X0505_MKL.1_scr_DMSO_EV" = "MKL-1 EV scr DMSO 0505",
  "X042_MKL.1_sT_Dox"       = "MKL-1 EV sT Dox 042",
  "X0505_MKL.1_sT_Dox_EV"   = "MKL-1 EV sT Dox 0505",
  "X042_MKL.1_scr_Dox_EV"   = "MKL-1 EV scr Dox 042",
  "X0505_MKL.1_scr_Dox_EV"  = "MKL-1 EV scr Dox 0505",
  "WaGa_RNA"                = "WaGa RNA",
  "WaGa_RNA_118"            = "WaGa RNA 118",
  "WaGa_RNA_147"            = "WaGa RNA 147",
  "WaGa_EV.RNA"             = "WaGa EV",
  "WaGa_EV.RNA_2"           = "WaGa EV 2",
  "WaGa_EV.RNA_118"         = "WaGa EV 118",
  "WaGa_EV.RNA_147"         = "WaGa EV 147",
  "WaGa_EV.RNA_226"         = "WaGa EV 226",
  "X1107_WaGa_wt_EV"        = "WaGa EV 1107",
  "X1605_WaGa_wt_EV"        = "WaGa EV 1605",
  "X2706_WaGa_wt_EV"        = "WaGa EV 2706",
  "X1107_WaGa_sT_DMSO_EV"   = "WaGa EV sT DMSO 1107",
  "X1605_WaGa_sT_DMSO_EV"   = "WaGa EV sT DMSO 1605",
  "X2706_WaGa_sT_DMSO_EV"   = "WaGa EV sT DMSO 2706",
  "X1107_WaGa_scr_DMSO_EV"  = "WaGa EV scr DMSO 1107",
  "X1605_WaGa_scr_DMSO_EV"  = "WaGa EV scr DMSO 1605",
  "X2706_WaGa_scr_DMSO_EV"  = "WaGa EV scr DMSO 2706",
  "X1107_WaGa_sT_Dox_EV"    = "WaGa EV sT Dox 1107",
  "X1605_WaGa_sT_Dox_EV"    = "WaGa EV sT Dox 1605",
  "X2706_WaGa_sT_Dox_EV"    = "WaGa EV sT Dox 2706",
  "X1107_WaGa_scr_Dox_EV"   = "WaGa EV scr Dox 1107",
  "X1605_WaGa_scr_Dox_EV"   = "WaGa EV scr Dox 1605",
  "X2706_WaGa_scr_Dox_EV"   = "WaGa EV scr Dox 2706")

stopifnot(all(colnames(reordered.raw) %in% names(name_map)),
          !anyDuplicated(unname(name_map[colnames(reordered.raw)])))
colnames(reordered.raw) <- unname(name_map[colnames(reordered.raw)])

# ---------------------------------------------------------------
# 4) Write out and filter (same as before)
# ---------------------------------------------------------------
write.csv(reordered.raw, file="counts.txt")

# IMPORTANT: filter low-count genes at this step!
d <- reordered.raw[rowSums(reordered.raw > 3) > 2, ]

What changed and why it is safer

Old code New code Why
row.names=1 at import Import first, compare, then set row names Allows checking the two Geneid/gene_name columns before one is discarded
No check of duplicated id columns identical() + stopifnot Catches a mis-pasted file (different gene order) immediately
Manual colnames(d.raw) <- c(...43 names...) R’s automatic sanitization (which yields the same names) No forced overwrite; names always come from the file header
Position-based colnames(reordered.raw) <- c(...) Explicit name_map keyed by actual column name A sample can never be mislabeled; any unknown/duplicate name throws an error
Filename "merged_gene_counts_40samples.txt " (trailing space) Corrected The stray space would make read.delim2 fail

Note: if stopifnot(identical(...)) fails, the two pasted files do not have the same gene order — in that case do not proceed; re-merge them properly (e.g. with a join on Geneid rather than paste/cbind).



Yes — I checked all six metadata vectors against the reordered column order (the 17 MKL-1 + 23 WaGa samples from your col_order/rename step). All vectors have exactly 40 entries, and every entry maps to the correct sample. The code is correct.

Verification summary

Vector MKL-1 block (17) WaGa block (23) Match?
condition_for_pca RNA×3, EV×6, sT.DMSO×2, scr.DMSO×2, sT.Dox×2, scr.Dox×2 RNA×3, EV×8, sT.DMSO×3, scr.DMSO×3, sT.Dox×3, scr.Dox×3
condition MKL1.RNA×3, MKL1.EV×6, then 2+2+2+2 WaGa.RNA×3, WaGa.EV×8, then 3+3+3+3
donor 1,118,147 / 1,2,118,87,27,042 / 042,0505 pairs 1,118,147 / 1,2,118,147,226,1107,1605,2706 / 1107,1605,2706 ×4
batch 2021.08, 2021.09×2 / 2021.08×2, 2021.09×3 / 2022.08×9 2021.08, 2021.09×2 / 2021.08×2, 2021.09×3 / 2022.11×15
cell.line MKL-1 ×17 WaGa ×23
ids unique per sample, order matches unique per sample, order matches

Spot checks that confirm the alignment:

  • Position 9: MKL-1 EV 042MKL1.EV, donor 042, batch 2022.08
  • Positions 26–28: WaGa EV 1107/1605/2706WaGa.EV, donors 1107/1605/2706, batch 2022.11
  • The skipped sample 0505_MKL-1_wt_EV is correctly absent (only one MKL-1 wt EV: 042) ✓

Run these assertions to be safe

stopifnot(ncol(d) == 40,
          length(condition) == 40, length(condition_for_pca) == 40,
          length(donor) == 40, length(batch) == 40,
          length(cell.line) == 40, length(ids) == 40)

# metadata consistent with sample names
stopifnot(all(grepl("^MKL-1", colnames(d)[cell.line == "MKL-1"])),
          all(grepl("^WaGa",  colnames(d)[cell.line == "WaGa"])))
stopifnot(all(paste0(ifelse(cell.line == "MKL-1", "MKL1", "WaGa"), ".") ==
              substr(as.character(condition), 1, 5)))

# visual check — eyeball once before DESeq2
data.frame(sample = colnames(d), condition, donor, batch, cell.line)

Three notes (not errors, but worth knowing)

  1. Design ~batch+condition is full rank here — good. The DMSO/Dox conditions exist only in batch 2022.08 (MKL-1) / 2022.11 (WaGa), but the wt-EV bridge samples (MKL-1 EV 042, WaGa EV 1107/1605/2706) connect those batches to MKL1.EV/WaGa.EV in the earlier batches, so all contrasts remain estimable. DESeq2 will stop with an error if this weren’t the case.
  2. donor is recorded but not used in the design. The sT/scr × DMSO/Dox samples are paired by donor (042/0505 and 1107/1605/2706), so a paired design (~batch+donor+condition) could add power — but given your note that you’re unsure whether “rep1” is the same donor across the 2021 and 2022 experiments, keeping ~batch+condition (as in the original analysis) is the safer, consistent choice.
  3. condition_for_pca is currently unused (it’s only referenced in a commented-out design line). Harmless, but you can drop it if you want to tidy up.

Everything downstream (vst(dds), estimateSizeFactors, the results(dds, name=...) contrast names like MKL1.sT.DMSO_vs_MKL1.scr.DMSO) will work exactly as in your established workflow.

中国版「生物DeepSeek」诞生:GeneLLM多组学大模型技术解析

https://www.vava8.com/index.php?app=index&act=view&id=122462

核心突破

由4位牛津大学归国博士创办的津渡生科,自主研发了生命科学垂类大模型 GeneLLM。该模型是全球首个直接基于组学原始数据进行预训练的多组学大模型,相关成果已发表于《Nature Communications》与《Advanced Science》,被视为继AlphaFold(结构预测)和EVO2(基因组理解)之后,向生命”系统级”理解迈进的关键一步。

技术架构与训练范式

Token化策略: GeneLLM将约150bp长度的RNA测序片段,通过7碱基滑动窗口(7-mer)切分为”生命Token”,以RNA四种碱基(A/U/G/C)作为基本语义单元,构建生命语言的离散表示。

无监督预训练: 模型采用Transformer架构,在无基因注释、无人工标签的条件下,直接对原始测序数据执行下一碱基预测(Next-Base Prediction)。训练过程处理了约数十万亿条RNA reads,在百卡NVIDIA A100集群上完成,使模型摆脱对已知基因注释的依赖,从原始信号中自主发现潜在生物模式。

两阶段训练流程:

  1. 无监督预训练与原型挖掘——从海量原始组学数据中学习通用生命表征;
  2. 患者级疾病微调(Disease Tuning)——面向具体疾病任务进行下游适配。

模型规模: 基础版已完成15亿参数、3.5万亿碱基序列预训练;XLarge版本扩展至300亿参数,持续扩大技术壁垒。

多组学覆盖: 训练数据涵盖RNA组、蛋白质组、代谢组等多组学原始数据,具备跨模态生命信息的联合理解能力。

推理效率: 传统方法依赖6Gb深度测序,GeneLLM在1Gb极浅深度测序下仍保持AUC > 0.8,测序成本降低约83%,为大规模临床筛查与普惠精准医疗提供工程可行性。

BioFord物理AI科研平台:从模型到闭环

以GeneLLM为认知底座,津渡生科构建了连接AI与物理实验的完整AI for Science闭环:

认知层——五大智能体协同网络: 文献检索智能体(自动文献综述与假设生成)、实验设计智能体(将实验设计周期从数月压缩至一周)、科学智能体(推理与假设验证)、实验调度智能体(基于通用仪器抽象层统一纳管PCR仪、酶标仪、流式细胞仪、自动化移液工作站等异构设备,内置动态调度算法实现自动排程与冲突规避)、数据分析智能体(结果解读与模型反馈)。

执行层——BioFord Harness: 将实验室转化为可编译、可调度、可观测、可追溯的工程系统,完成三项关键任务:将科学意图或实验DSL编译为跨设备可执行指令;在多设备间完成调度、资源约束管理与异常处理;将实验结果、设备日志和环境参数回流至模型,驱动下一轮迭代。

数据闭环——DBTL循环: 通过Design-Build-Test-Learn闭环,每一次实验(包括失败实验)的参数逻辑、环境记录与错误路径均沉淀为训练数据,形成持续进化的科研数据飞轮。

全球主要路线对比

公司/平台 路线 核心能力 当前阶段
FutureHouse AI科学大脑 文献理解、科学推理 数字科研Agent
DeepMind 基础科学模型 生物结构预测 科学基础模型
Insilico Medicine AI药物研发 靶点发现、分子设计 AI制药平台
XtalPi晶泰科技 AI+机器人实验 自动化药物研发 实验闭环
Lila Sciences AI Science Factory 自动化科学工厂 重资产实验室
Recursion 生物数据工业化 大规模细胞实验 数据驱动研发
Isomorphic Labs AI药物设计 AlphaFold路线延伸 分子发现
津渡生科 生命大模型+科研Agent AI理解生命+执行实验 全栈AI科研OS

在这一格局中,津渡生科选择轻量化物理AI路线:不做纯数字AI科学家,不做重资产科学工厂,不做端到端制药管线,而是专注打通模型与实体实验室之间的”最后一公里”基础设施,以实验轨迹、设备接口协议、失败经验库和跨实验室执行网络构建核心护城河。

资本

公司一年内完成4轮融资,投资方包括红杉中国种子基金、创东方投资、南山战新投及高特佳投资(A轮近亿元领投)。

“AI for Science真正的分水岭,不是模型回答得多像科学家,而是实验室能否开始像一个持续学习的系统。”

Comprehensive Guide for GEO Submission (Data_Ute_smallRNA_via_exceRpt_workspace_FINAL)

For MKL-1, the two files are complementary, not interchangeable.

  • exceRpt_biotypeCounts.txt = biotype abundance / composition table
  • exceRpt_mapping_heatmaps_MKL-1.xlsx = mapping/QC summary table

Comparison table for MKL-1 files

Feature exceRpt_biotypeCounts.txt exceRpt_mapping_heatmaps_MKL-1.xlsx
Main purpose Shows how many reads/abundance estimates were assigned to different RNA biotypes Shows how reads progressed through QC, trimming, alignment, and mapping categories
Data type Biotype-level abundance table Mapping/QC fraction table, likely normalized to input reads
MKL-1 samples included Same MKL-1 sample set: 2404_MKL1_wt_EVs, 2608_MKL1_scr_DMSO, 2608_MKL1_scr_Dox, 2608_MKL1_sT_DMSO, 2608_MKL1_sT_Dox, 2608_MKL1_wt_EVs, 2701_MKL1_scr_DMSO, 2701_MKL1_scr_Dox, 2701_MKL1_sT_DMSO, 2701_MKL1_sT_Dox, 2802_MKL1_scr_DMSO, 2802_MKL1_scr_Dox, 2802_MKL1_sT_DMSO, 2802_MKL1_sT_Dox, plus nf780, nf796, nf797 Same sample list as the biotypeCounts file
Format Plain text / tab-delimited table Excel file
Values Numeric abundance values for RNA biotypes. Some values are fractional, suggesting normalized or fractional assignment rather than simple integer raw counts Values appear to be fractions/proportions, with input = 1
Main biological categories miRNA, tRNA, piRNA, snRNA, snoRNA, rRNA, protein_coding, lincRNA, retained_intron, processed_transcript, antisense, misc_RNA, exogenous_genomes, exogenous_miRNA, exogenous_rRNA, circularRNA, etc. input, successfully_clipped, failed_quality_filter, failed_homopolymer_filter, UniVec_contaminants, rRNA, reads_used_for_alignment, genome, miRNA_sense/antisense, tRNA_sense/antisense, piRNA_sense/antisense, gencode_sense/antisense, circularRNA_sense/antisense, not_mapped_to_genome_or_libs, repetitiveElements, exogenous_genomes, etc.
Best used for Small RNA composition plots, e.g. miRNA / tRNA / piRNA / long RNA biotype percentages Mapping efficiency, QC filtering, alignment statistics, and reproducibility of how reads were distributed
Most relevant manuscript panel Figure 4A and Supplementary Figure S5A: small RNA biotype composition Mapping/QC text, e.g. percentage of reads mapped, and supplementary QC information
Strength Directly supports the biological small RNA composition results Directly supports mapping quality and reproducibility
Limitation Does not show QC/mapping steps such as adapter clipping, quality filtering, unmapped reads, etc. Does not provide the full biotype abundance table needed to reproduce Figure 4A/S5A
Repository suitability High, if converted to clean CSV/TSV and accompanied by sample metadata Moderate to high, but should be converted from Excel to CSV/TSV and clearly labeled as mapping/QC summary
Enough for miRNA-level figures? No. It gives total miRNA abundance, but not individual miRNA counts/RPM No. It gives mapping fractions, not individual miRNA abundance

Which file should you submit?

Best recommendation

Submit both, but with different roles:

File Submit? Role
exceRpt_biotypeCounts.txt Yes — main processed small RNA biotype file Supports small RNA biotype composition, e.g. Figure 4A / Suppl. Fig. S5A
exceRpt_mapping_heatmaps_MKL-1.xlsx Yes — secondary mapping/QC file Supports mapping efficiency and QC reproducibility

If you can upload only one file, submit:

exceRpt_biotypeCounts.txt

because it is closer to the actual biological small RNA composition results shown in the manuscript.

However, the best processed-data package would be:

smallRNA_exceRpt_biotypeCounts_MKL-1.csv
smallRNA_exceRpt_mapping_summary_MKL-1.csv
smallRNA_sample_metadata_MKL-1.txt

Important practical suggestions before submission

  1. Convert the Excel mapping file to CSV/TSV

    Many repositories prefer plain text files. For example:

    exceRpt_mapping_heatmaps_MKL-1.xlsx

    could become:

    smallRNA_exceRpt_mapping_summary_MKL-1.csv
  2. Clarify the value type in exceRpt_biotypeCounts.txt

    The values are not all integers. Before submission, check whether this file contains:

    • raw read counts,
    • normalized counts,
    • RPM/CPM,
    • or exceRpt fractional assignment values.

    Add a short README or column description, for example:

    Values are exceRpt-derived biotype abundance estimates.

    or, if confirmed:

    Values are read counts assigned to RNA biotypes by exceRpt.
  3. Define the mapping categories

    In the mapping heatmap file, categories such as:

    • reads_used_for_alignment
    • genome
    • not_mapped_to_genome_or_libs
    • gencode_sense
    • miRNA_sense

    should be explained in a README. Otherwise reviewers may not know which row corresponds to “reads mapped to the human genome”.

  4. Decide what to do with nf780, nf796, and nf797

    These samples are included in both files but are not MKL-1 EV samples. If they are not part of the manuscript figures, you should either:

    • remove them from the MKL-1 processed table, or
    • keep them but clearly annotate them in the metadata as non-MKL-1/control samples.
  5. You still need an individual miRNA table

    Neither of these two files is sufficient for the individual miRNA Manhattan plots or miRNA-level analyses, e.g. Figure 4B or Supplementary Figure S5B. For those, you should also provide:

    smallRNA_miRNA_counts_all_samples.txt
    smallRNA_miRNA_RPM_normalized_all_samples.txt

    if those were used for the figures.


Final suggestion

For MKL-1:

Primary processed file to submit:
exceRpt_biotypeCounts.txt
→ rename to smallRNA_exceRpt_biotypeCounts_MKL-1.csv

Secondary processed file to submit:
exceRpt_mapping_heatmaps_MKL-1.xlsx
→ convert to smallRNA_exceRpt_mapping_summary_MKL-1.csv

If only one file can be submitted, choose exceRpt_biotypeCounts.txt.
If you want full reproducibility, submit both, plus a sample metadata file.



Yes — but you do not necessarily need to submit every intermediate pipeline file.
You should submit the final processed files that are required to reproduce the quantitative results/figures in the manuscript.

By “raw sequencing data” I assume you mean FASTQ/BAM. Files such as *_raw_counts*.txt are already processed data relative to FASTQ, and they are usually expected for GEO/journal submission if they underlie the figures.

Below is a careful breakdown based on the manuscript text.


1. Processed data used in the manuscript that you should consider submitting

# Processed data type Where it is used in the manuscript Original manuscript sentence / relevant text Suggested file(s) to submit
1 RNA-seq gene-level count matrix Methods 4.10; Results on EV RNA cargo; Fig. 3; Suppl. Fig. S4, S9, S11 “A total of 40 RNA-seq libraries were processed using the nf-core/rnaseq pipeline…”
“Gene-level read counts were generated using featureCounts…”
“Raw count data were analyzed using DESeq2.”
RNAseq_raw_counts_all_samples.txt
2 RNA-seq normalized / VST-transformed counts Methods 4.10; heatmaps/PCA/clustering; Fig. 6D; Suppl. Fig. S9 “For visualization, count data were normalized and variance-stabilized using the variance stabilizing transformation (VST).” RNAseq_VST_normalized_counts_all_samples.txt or similar
3 RNA-seq differential abundance tables Fig. 3B; Fig. 6D; Suppl. Fig. S9; Results on EV vs parental cells and sT knockdown “Approximately 26,000 transcripts showed a higher relative abundance in EVs, whereas approximately 2,800 transcripts exhibited a lower relative abundance compared with the parental cells.”
“Fifteen transcripts showed significantly higher relative abundance following sT knockdown…”
RNAseq_DESeq2_EV_vs_parental_results.txt
RNAseq_DESeq2_sT_knockdown_results.txt
4 Viral MCPyV transcript counts / normalized abundance Methods 4.10; Suppl. Fig. S11 “Sequencing reads were aligned using STAR against a combined reference comprising the human genome (GRCh38) and the corresponding Merkel cell polyomavirus (MCPyV) genome…”
Suppl. Fig. S11: “Red dots indicate MCPyV-derived viral transcripts…”
RNAseq_MCPyV_viral_transcript_counts.txt or include viral rows in the main RNA-seq count table
5 small RNA-seq miRNA count matrix Methods 4.11; Fig. 4; Fig. 5; Suppl. Fig. S5, S7 “Known human miRNAs were annotated according to miRBase, and read counts for individual miRNAs were generated using the COMPSRA pipeline.”
“Raw miRNA count data were analyzed using DESeq2.”
smallRNA_miRNA_counts_all_samples.txt
6 small RNA-seq normalized miRNA abundance Fig. 4B; Suppl. Fig. S5B; possibly S7 Fig. 4B legend: “Manhattan plots showing normalized miRNA abundance(log10 reads per million; RPM)…” smallRNA_miRNA_RPM_normalized_all_samples.txt or VST/DESeq2-normalized miRNA table
7 small RNA biotype composition table Fig. 4A; Suppl. Fig. S5A; Results on small RNA composition “In parental cells, 59% of mapped reads were annotated as miRNAs…”
“In WaGa-derived EVs, 28% of mapped reads were annotated as miRNAs… tRNA-derived reads increased to approximately 29% and piRNA-derived reads accounted for 0.65%.”
smallRNA_biotype_counts_summary.txt or the source count table used to generate Fig. 4A/S5A
8 small RNA mapping summary Results on small RNA mapping; Fig. 4A context “Approximately 98% of reads obtained from WaGa cells mapped to the human genome, whereas approximately 73% of reads from EVs could be mapped to the human genome.” smallRNA_mapping_summary.txt
9 RBP motif enrichment results Results on RNA-binding protein motifs; Fig. 3 / Fig. S4 “Analysis of Motif Enrichment(AME) was performed using the ATtRACT database.”
“Several significantly enriched sequence motifs and their corresponding RBPs were identified.”
RNAseq_RBP_motif_enrichment_results.txt
10 miRNA target network table Results 2.5; Fig. 5; Suppl. Fig. S6 “…experimentally validated target genes of the 15 most abundant miRNAs identified in WaGa-derived EVs were retrieved from miRTarBase and used to construct a miRNA-target interaction network.”
“The resulting network comprised 196 target genes.”
miRNA_target_network_nodes.txt
miRNA_target_network_edges.txt
11 Proteomics protein identification and quantification table Results proteomics; Fig. 2; Fig. 6; Suppl. Fig. S3, S10; Data Availability “To characterize the protein cargo of WaGa-derived EVs, mass spectrometry was performed, identifying 608 proteins consistently detected across all biological replicates.”
“The proteomics data have been deposited with the ProteomeXchange Consortium via the PRIDE partner repository…”
Proteomics_protein_identifications_quantification.txt
optionally peptide-level table
12 Proteomics differential abundance table Fig. 6C; Suppl. Fig. S10; Results on sT knockdown proteome “Differential abundance analysis identified 26 proteins whose abundance was significantly altered following sT knockdown, comprising 24 proteins with higher abundance and 2 proteins with lower abundance in EVs.” Proteomics_DE Proteins_sT_knockdown_WaGa.txt
Proteomics_DE Proteins_sT_knockdown_MKL-1.txt

2. Do you really need to submit the six files you listed?

File My recommendation Reason
RNAseq_raw_counts_all_samples.txt Yes — submit This is the count matrix used for DESeq2 analysis and underlies Fig. 3, Fig. 6D, Suppl. Fig. S9, and likely Suppl. Fig. S11. The manuscript explicitly says: “Raw count data were analyzed using DESeq2.”
smallRNA_miRNA_counts_all_samples.txt Yes — submit This is the main small RNA count matrix used for miRNA abundance, differential miRNA analysis, Fig. 4, Fig. 5, Suppl. Fig. S5 and S7. The manuscript explicitly says miRNA read counts were generated and analyzed with DESeq2.
smallRNA_gencode_counts_all_samples.txt Conditional Submit only if this file was used to generate the small RNA biotype composition shown in Fig. 4A / Suppl. Fig. S5A, e.g. miRNA, tRNA, piRNA, long RNA, etc. The manuscript does not explicitly mention “GENCODE counts”, so it is not automatically required. If this file is only an intermediate file and the biotype percentages were generated from another summary table, you can submit the summary table instead.
smallRNA_tRNA_counts_all_samples.txt Conditional / optional The manuscript only reports tRNAs as a class: “tRNAs represented 29% of mapped reads…” It does not appear to analyze individual tRNA-derived fragments in detail. If this file is needed to calculate the tRNA percentage in Fig. 4A/S5A, submit it or include it in a combined biotype table. If not, a summarized biotype count table is enough.
smallRNA_piRNA_counts_all_samples.txt Conditional / optional Same logic as tRNA. The manuscript only reports piRNAs as a percentage: “piRNAs accounted for 0.65%.” If this file is required to reproduce that number, submit it or include it in a biotype summary. Otherwise, individual piRNA counts are not strictly needed unless they are analyzed elsewhere.
smallRNA_mapping_summary.txt Recommended, but not always mandatory This supports the mapping rates and possibly the biotype composition reported in the manuscript: “Approximately 98% of reads… mapped… whereas approximately 73%…” It is very useful for reproducibility. I would include it as a supplementary processed file or as part of a small RNA QC/summary table.

3. My practical recommendation for the minimal processed-data package

I would submit the following as the clean, final processed files:

RNA-seq

File Needed? Why
RNAseq_raw_counts_all_samples.txt Yes Underlies DESeq2 analysis
RNAseq_normalized_counts_all_samples.txt Strongly recommended The manuscript says VST-normalized data were used for visualization
RNAseq_DESeq2_results_EV_vs_parental.txt Yes Underlies volcano plots / differential RNA cargo
RNAseq_DESeq2_results_sT_knockdown.txt Yes Underlies Fig. 6D / Suppl. Fig. S9
RNAseq_MCPyV_transcript_counts_or_RPM.txt Recommended Needed if Suppl. Fig. S11 is shown
RNAseq_sample_metadata.txt Yes Needed to understand conditions, replicates, EV vs cells, Dox/DMSO, WaGa/MKL-1

small RNA-seq

File Needed? Why
smallRNA_miRNA_counts_all_samples.txt Yes Main miRNA count matrix
smallRNA_miRNA_RPM_normalized_all_samples.txt Strongly recommended Fig. 4B and Suppl. Fig. S5B show normalized miRNA abundance in RPM
smallRNA_biotype_counts_summary.txt Yes / strongly recommended Underlies Fig. 4A and Suppl. Fig. S5A
smallRNA_mapping_summary.txt Recommended Supports mapping rates reported in the text
smallRNA_DESeq2_miRNA_results.txt Yes if differential miRNA analysis is shown Suppl. Fig. S7 / related text
smallRNA_sample_metadata.txt Yes Needed to interpret samples

If you create one clean file called, for example:

smallRNA_biotype_counts_summary.txt

with columns like:

sample
cell_line
sample_type
treatment
total_reads
mapped_reads
miRNA_counts
miRNA_percent
tRNA_counts
tRNA_percent
piRNA_counts
piRNA_percent
longRNA_counts
longRNA_percent
other_counts
other_percent

then you may not need to submit separate smallRNA_tRNA_counts_all_samples.txt and smallRNA_piRNA_counts_all_samples.txt, unless those files are the direct source of the figure.


4. Proteomics: do not forget the processed MS files

The manuscript says:

“The proteomics data have been deposited with the ProteomeXchange Consortium via the PRIDE partner repository under accession number PXDXXXXX.”

And your manuscript notes also say:

“Bente has to do this: Upload processed files (e.g., peptide/protein identifications, quantification tables). Include metadata describing the experimental design.”

So for proteomics, you should submit at least:

Proteomics processed file Needed?
Protein identification table Yes
Protein quantification table Yes
Peptide identification table Recommended / often expected by PRIDE
Differential protein abundance table for sT knockdown Yes
Sample metadata for MS runs Yes

5. Important consistency issues before submission

You should also check these points, because they affect which processed files are correct to submit.

A. Small RNA pipeline inconsistency

In the Results section, the manuscript says:

“Sequencing data were analyzed using the exceRpt pipeline, which is optimized for extracellular RNA analysis.”

But in Methods 4.11, the manuscript says:

“Raw FASTQ files generated by small RNA sequencing were processed using Cutadapt…”
“High-quality reads were aligned… using COMPSRA with the STAR aligner.”

You need to decide which pipeline actually produced the final count files:

  • If the final files came from exceRpt, the Methods should mention exceRpt.
  • If the final files came from COMPSRA/Cutadapt, remove or correct the exceRpt sentence.

This matters because the submitted processed files must match the described pipeline.

B. Small RNA reference genome inconsistency

Methods 4.11 says:

“…combined reference comprising the human genome (GRCh38) and the Merkel cell polyomavirus genome(JN707599)…”

But for RNA-seq, Methods 4.10 says:

“…combined reference comprising the human genome(GRCh38) and the corresponding Merkel cell polyomavirus(MCPyV) genome(KJ128379.1 for WaGa samples and FJ173815.1 for MKL-1 samples).”

Check whether small RNA-seq really used only JN707599 for all samples, or whether it should match the WaGa/MKL-1 references.

C. Figure 3 panel order

In the Results text, Figure 3C and 3D may not match the Figure legend perfectly. The Results text describes GO enrichment and then RBP motif enrichment, while the Figure legend appears to assign them differently. Before submission, make sure the text, figure panels, and supplementary files match.


Bottom line

For your listed files:

RNAseq_raw_counts_all_samples.txt              → Yes, submit
smallRNA_miRNA_counts_all_samples.txt          → Yes, submit
smallRNA_gencode_counts_all_samples.txt        → Only if used to generate the reported small RNA biotype composition
smallRNA_tRNA_counts_all_samples.txt           → Optional unless needed to reproduce Fig. 4A/S5A
smallRNA_piRNA_counts_all_samples.txt          → Optional unless needed to reproduce Fig. 4A/S5A
smallRNA_mapping_summary.txt                   → Recommended, useful for reproducibility

The safest approach is to submit a clean final processed-data package containing:

  1. RNA-seq raw counts
  2. RNA-seq normalized counts / VST
  3. RNA-seq differential expression tables
  4. small RNA miRNA raw counts
  5. small RNA miRNA normalized RPM
  6. small RNA biotype summary
  7. small RNA mapping summary
  8. proteomics protein/peptide identification and quantification tables
  9. sample metadata for all omics datasets


1. Do You Need to Submit to GEO and Processed Data?

Yes, absolutely. Based on the manuscript’s Data Availability Statement and the email from Nicole Fischer:

Data Type Repository Accession Format
RNA-seq (total RNA) NCBI GEO (linked to SRA) GSEXXXXX / SRXXXXX
Small RNA-seq NCBI GEO (linked to SRA) GSEXXXXX / SRXXXXX
Proteomics (LC-MS/MS) PRIDE / ProteomeXchange PXDXXXXX

What GEO requires:

  • Raw data: All .fastq.gz files → uploaded to SRA (Sequence Read Archive) via GEO
  • Processed data: YES, required. GEO mandates at least one processed data file per sample or a summary matrix. For RNA-seq this means:
    • Gene/transcript count matrix (raw counts from featureCounts)
    • Optionally: normalized counts (VST/TPM/FPKM)
    • For small RNA-seq: miRNA count matrix
  • Metadata: Sample attributes, experimental design, platform info, protocols

⚠️ Critical: The email says data must be uploaded but NOT released publicly yet. In GEO you set a future release date (e.g., 1–2 years from now) so you get the accession number immediately but data stays private until manuscript publication.


数据可用性声明

本研究生成的RNA测序(RNA-seq)和小RNA测序(small RNA-seq)数据集已存入基因表达综合数据库(GEO),登录号为GSEXXXXX。蛋白质组学数据已存入ProteomeXchange联盟的PRIDE合作存储库,登录号为PXDXXXXX。


3. Step-by-Step GEO Submission Guide (Einleitung)

Phase 0: Preparation Checklist

Before you start, gather:

Item Details
NCBI Account Register at ncbi.nlm.nih.gov/account
GEO Account After NCBI login, request GEO submission access at geo@ncbi.nlm.nih.gov
Raw FASTQ files All .fastq.gz files (RNA-seq + small RNA-seq)
Processed data Count matrices (featureCounts output for RNA-seq; COMPSRA miRNA counts for small RNA-seq)
Sample metadata Cell line, condition, replicate, sample type (EV vs parental cell)
Protocol info Library prep kits, sequencing platform, read length

Phase 1: Create GEO Submitter Account

  1. Go to https://www.ncbi.nlm.nih.gov/ → Sign in / Register
  2. Once logged in, email geo@ncbi.nlm.nih.gov with:
    • Subject: “GEO Submitter Account Request”
    • Your name, institution, email
    • State you want to submit RNA-seq and small RNA-seq data
  3. Wait for confirmation (usually 1–2 business days). You will receive a GEO submitter login.

Phase 2: Organize Your Samples & Metadata

Based on the manuscript, organize samples into a clear table. Here is the suggested sample organization:

RNA-seq Samples (WaGa)

Sample Name Cell Line Sample Type Condition Replicate File
WaGa_wt_EV_rep1 WaGa EV untreated (wt) 1 1107_WaGa_wt_EV.fastq.gz
WaGa_wt_EV_rep2 WaGa EV untreated (wt) 2 1605_WaGa_wt_EV.fastq.gz
WaGa_wt_EV_rep3 WaGa EV untreated (wt) 3 2706_WaGa_wt_EV.fastq.gz
WaGa_scr_DMSO_EV_rep1 WaGa EV scr + DMSO 1 1107_WaGa_scr_DMSO_EV.fastq.gz
WaGa_scr_Dox_EV_rep1 WaGa EV scr + Dox 1 1107_WaGa_scr_Dox_EV.fastq.gz
WaGa_sT_DMSO_EV_rep1 WaGa EV sT + DMSO 1 1107_WaGa_sT_DMSO_EV.fastq.gz
WaGa_sT_Dox_EV_rep1 WaGa EV sT + Dox (sT KD) 1 1107_WaGa_sT_Dox_EV.fastq.gz
… (rep2, rep3)
WaGa_cell_RNA_rep1 WaGa Parental cell untreated 1 WaGa_RNA.fastq.gz
WaGa_cell_RNA_rep2 WaGa Parental cell untreated 2 WaGa_RNA_118.fastq.gz
WaGa_cell_RNA_rep3 WaGa Parental cell untreated 3 WaGa_RNA_147.fastq.gz

RNA-seq Samples (MKL-1)

Same structure as WaGa.

Small RNA-seq Samples (WaGa)

Map nf774, nf930nf939, nf961, nf962, nf971nf974 to their corresponding conditions. You need to check your lab records to confirm which nf-number corresponds to which sample.

Small RNA-seq Samples (MKL-1)

Map 2404_MKL1_wt_EVs, 2608_*, 2701_*, 2802_*, nf780, nf796, nf797 similarly.

⚠️ Important: The file naming is inconsistent (some have EV suffix, some don’t; some use dates like 042/0505, others use nf-numbers). You must verify the sample-to-file mapping with Ute or your lab notebook before submission.

Phase 3: Start GEO Submission via GEO Submitter Portal

  1. Log in to the GEO Submitter Portal: https://www.ncbi.nlm.nih.gov/geo/submitter/
  2. Click “Submit” → Select “High-Throughput Sequencing” (for RNA-seq and small RNA-seq)
  3. You will create a GEO Series (GSE) that contains:
    • BioProject (auto-created)
    • BioSample entries (one per biological sample)
    • SRA entries (linked FASTQ files)

Phase 4: Fill in the Submission Metadata

4a. Series (GSE) Level

Field What to Enter
Title Multi-omics characterization of extracellular vesicles derived from virus-positive Merkel cell carcinoma cells – RNA-seq and small RNA-seq
Summary Copy/adapt from manuscript abstract
Overall design Two MCPyV-positive MCC cell lines (WaGa and MKL-1) with doxycycline-inducible sT knockdown. EVs isolated by differential ultracentrifugation. Total RNA-seq and small RNA-seq of EVs and parental cells under sT knockdown and control conditions. Three biological replicates per condition.
Experiment type RNA-Seq; small RNA-seq
Release date Set to a future date (e.g., 2027-08-01 or later) to keep data private

4b. BioSample Attributes (per sample)

For each sample, provide:

  • organism: Homo sapiens
  • cell_line: WaGa or MKL-1
  • sample_type: extracellular vesicle / parental cell
  • treatment: doxycycline / DMSO / untreated
  • genotype: shRNA-sT / shRNA-scramble / wild-type
  • molecule: total RNA / small RNA
  • source_name: e.g., “WaGa EV sT knockdown replicate 1”

4c. SRA / Platform Info

Field Value
Instrument Illumina NextSeq 500
Read length 75 bp single-end
Library strategy RNA-Seq / miRNA-Seq
Library source transcriptomic
Library selection cDNA / size fractionation (small RNA)
Library kit CORALL Total RNA-Seq V2 Kit / LEXOGEN Small RNA-Seq Library Prep Kit

Phase 5: Upload Raw FASTQ Files

  1. GEO submission will direct you to upload files to SRA via FTP or Aspera
  2. You will receive an FTP upload folder (e.g., ftp://ftp-private.ncbi.nlm.nih.gov/uploads/geo/...)
  3. Upload all .fastq.gz files:
# Example using FTP (you can also use Aspera or web browser)
# Connect to the FTP server provided by GEO
ftp ftp-private.ncbi.nlm.nih.gov
# Login with provided credentials
cd uploads/geo/your_folder/

# Upload WaGa RNA-seq
put Data_RNA-Seq_WaGa/1107_WaGa_wt_EV.fastq.gz
put Data_RNA-Seq_WaGa/1107_WaGa_scr_DMSO_EV.fastq.gz
# ... repeat for all files

# Upload small RNA-seq
put Data_smallRNA_WaGa/nf774.fastq.gz
# ... repeat for all files

Or use Aspera (faster for large files):

ascp -k 1 -T -l 300m \
  Data_RNA-Seq_WaGa/*.fastq.gz \
  subasp@upload.ncbi.nlm.nih.gov:uploads/geo/your_folder/RNASeq_WaGa/

ascp -k 1 -T -l 300m \
  Data_smallRNA_WaGa/*.fastq.gz \
  subasp@upload.ncbi.nlm.nih.gov:uploads/geo/your_folder/smallRNA_WaGa/

Phase 6: Upload Processed Data

GEO requires processed data. Prepare:

  1. RNA-seq count matrix (from featureCounts):

    • A tab-delimited file: rows = genes, columns = samples
    • Raw counts (not normalized)
  2. Small RNA-seq miRNA count matrix (from COMPSRA):

    • A tab-delimited file: rows = miRNAs, columns = samples
  3. Upload these as supplementary files in the GEO submission portal:

    • Format: .txt or .csv (tab-delimited)
    • Label clearly: e.g., WaGa_MKL1_RNAseq_raw_counts.txt, WaGa_MKL1_smallRNA_miRNA_counts.txt

Phase 7: Review and Submit

  1. Review all metadata, sample descriptions, and file assignments
  2. Set the release date to a future date (this keeps data private)
  3. Click Submit
  4. You will receive:
    • GSE accession number (e.g., GSE123456) → put this in the manuscript
    • SRA accession numbers for individual samples
    • A token for reviewer access (optional, for peer review)

Phase 8: Update the Manuscript

Replace GSEXXXXX in the Data Availability Statement with the actual GSE number:

“The RNA sequencing (RNA-seq) and small RNA sequencing (small RNA-seq) datasets generated in this study have been deposited in the Gene Expression Omnibus (GEO) under accession number GSE[ACTUAL NUMBER].”


Quick Checklist Summary

Step Action Status
☐ 1 Register NCBI + GEO account
☐ 2 Verify sample-to-file mapping (especially nf-numbers!)
☐ 3 Prepare metadata spreadsheet
☐ 4 Generate processed count matrices
☐ 5 Start GEO submission (High-Throughput Sequencing)
☐ 6 Fill in Series, BioSample, SRA metadata
☐ 7 Upload FASTQ files via FTP/Aspera
☐ 8 Upload processed data files
☐ 9 Set future release date (NOT public)
☐ 10 Submit → receive GSE accession number
☐ 11 Send GSE number to Nicole/Ute for manuscript

💡 Tip: The proteomics data (PRIDE/ProteomeXchange) is a separate submission handled by the mass spectrometry team (Bente Siebels). You only need to handle the RNA-seq and small RNA-seq GEO submission.