长读长测序数据分析（Nanopore / PacBio）¶

原始电信号（POD5/FAST5）
   ↓ Basecalling（Dorado）
FASTQ/BAM
   ↓ 质控（NanoPlot / NanoQC）
   ↓ 比对（minimap2）
sorted BAM
   ↓ 变异检测（Clair3 / PEPPER-DeepVariant）
   ↓ 结构变异（Sniffles2 / SVIM）
   ↓ 基因组组装（Flye → polish with Medaka）
   ↓ 甲基化检测（Dorado 内置 / modkit）

# ============================================================
# Dorado basecalling：把电信号（POD5格式）转换为碱基序列
# ============================================================

# 安装 Dorado（官方提供预编译二进制，需要 GPU）
# 下载地址: https://github.com/nanoporetech/dorado/releases

# 查看可用模型（不同flowcell/化学版本用不同模型）
dorado download --list              # 列出所有可下载的模型

# 下载模型（以 R10.4.1 flowcell + 最高精度模型为例）
dorado download --model dna_r10.4.1_e8.2_400bps_sup@v5.0.0
# sup = super accuracy 最高精度模型（最慢但最准）
# hac = high accuracy 高精度模型（速度和准确率平衡）
# fast = 快速模型（最快但准确率稍低）

# 执行 basecalling
dorado basecaller \
    dna_r10.4.1_e8.2_400bps_sup@v5.0.0 \   # 指定模型
    ./pod5_pass/ \                           # 输入：POD5 文件目录
    --emit-fastq \                           # 输出 FASTQ 格式（默认输出 BAM）
    --device cuda:0 \                        # 使用第一块 GPU
    --min-qscore 9 \                         # 最低质量分数过滤（Q9 ≈ 87%准确率）
    > basecalled.fastq                       # 输出文件

# 如果同时需要甲基化信息（5mC 和 6mA）：
dorado basecaller \
    dna_r10.4.1_e8.2_400bps_sup@v5.0.0 \
    ./pod5_pass/ \
    --modified-bases 5mCG_5hmCG \            # 检测 CpG 位点的 5mC 和 5hmC 修饰
    --device cuda:all \                      # 使用所有 GPU
    --reference ref.fasta \                  # 直接比对到参考基因组
    > aligned_with_mods.bam                  # 输出含甲基化标签的 BAM

# ============================================================
# NanoPlot：长读长专用质控工具，生成读长分布、质量分布等图表
# ============================================================

# 安装
pip install NanoPlot  # 或 conda install -c bioconda nanoplot

# 从 FASTQ 文件生成质控报告
NanoPlot \
    --fastq basecalled.fastq \    # 输入 FASTQ 文件
    --outdir nanoplot_qc/ \       # 输出目录
    --threads 8 \                 # 使用 8 个线程
    --plots hex dot \             # 生成的图表类型
    --title "Sample_01_QC" \      # 报告标题
    --loglength                   # 用 log 刻度展示读长（长读长数据分布跨度大）

# 关键质控指标解读：
# - Mean read length: 平均读长（期望 >10kb）
# - N50: 把所有读长从长到短排列，累积到总碱基数 50% 时的读长
#         （N50 越大说明长 reads 越多，组装越有利）
# - Mean read quality: 平均质量（Q10=90%准确率, Q20=99%准确率）
# - Number of reads: 总读数

# 从 BAM 文件生成质控（如果 Dorado 直接输出 BAM）
NanoPlot --bam aligned.bam --outdir nanoplot_bam_qc/ --threads 8

# ============================================================
# minimap2：长读长比对的标准工具（速度快，支持各种长读长数据）
# ============================================================

# 安装
conda install -c bioconda minimap2 samtools

# Nanopore 数据比对到参考基因组
minimap2 \
    -ax map-ont \              # -a: 输出 SAM 格式; -x map-ont: Nanopore 预设参数
    -t 16 \                    # 使用 16 个线程
    --secondary=no \           # 不输出次优比对（减少冗余）
    --MD \                     # 输出 MD tag（方便后续变异检测）
    ref.fasta \                # 参考基因组
    basecalled.fastq \         # 输入 reads
    | samtools sort \          # 管道：排序
    -@ 8 \                     # 排序用 8 个线程
    -o aligned_sorted.bam      # 输出排序后的 BAM

# 建立索引（后续查看/变异检测需要）
samtools index aligned_sorted.bam

# minimap2 的 -x 预设参数说明：
# map-ont   → Nanopore 基因组比对
# map-hifi  → PacBio HiFi 基因组比对
# map-pb    → PacBio CLR（旧版连续长读长）比对
# splice    → 长读长转录组比对（RNA-seq，考虑剪接）
# asm5      → 基因组-基因组比对（序列差异 <5%）

# ============================================================
# Clair3：长读长 SNP/Indel 检测（深度学习模型）
# ============================================================

# 推荐用 Docker/Singularity 运行（依赖复杂）
# Docker 方式：
docker run -it \
    -v "$(pwd):/input" \                    # 挂载当前目录
    -v "$(pwd)/output:/output" \            # 输出目录
    hkubal/clair3:latest \
    /opt/bin/run_clair3.sh \
    --bam_fn=/input/aligned_sorted.bam \    # 输入 BAM
    --ref_fn=/input/ref.fasta \             # 参考基因组
    --output=/output/ \                     # 输出目录
    --threads=16 \                          # 线程数
    --platform="ont" \                      # 平台：ont / hifi / ilmn
    --model_path="/opt/models/ont" \        # 模型路径（ont专用模型）
    --sample_name="sample01"                # 样本名

# 输出文件说明：
# merge_output.vcf.gz      → 合并的变异结果（SNP + Indel）
# pileup.vcf.gz            → pileup 阶段初筛结果
# full_alignment.vcf.gz    → 全比对阶段精细结果

# Conda 方式安装（替代 Docker）：
# conda create -n clair3 -c bioconda clair3
# run_clair3.sh --bam_fn ... --ref_fn ... --platform ont --output ...

# ============================================================
# Flye：长读长 de novo 基因组组装
# ============================================================

# 安装
conda install -c bioconda flye

# Nanopore 数据组装
flye \
    --nano-hq basecalled.fastq \    # --nano-hq: 高质量 Nanopore reads（Q20+）
                                    # --nano-raw: 旧版低质量 reads
                                    # --pacbio-hifi: PacBio HiFi reads
    --out-dir flye_assembly/ \      # 输出目录
    --threads 32 \                  # 线程数
    --genome-size 5m \              # 预估基因组大小（5m=5Mb，细菌基因组）
                                    # 人类基因组用 3g
    --iterations 2                  # polishing 迭代次数

# Flye 输出文件说明：
# assembly.fasta           → 组装结果
# assembly_graph.gfa       → 组装图（可用 Bandage 可视化）
# assembly_info.txt        → 各 contig 的统计信息

# ============================================================
# Medaka：Nanopore 专用 polishing 工具（ONT官方，基于神经网络）
# ============================================================

# 安装
# pip install medaka  # 或 conda install -c bioconda medaka

# 用原始 reads 对组装结果进行 polish（纠错）
medaka_polish \
    -i basecalled.fastq \           # 原始 reads
    -d flye_assembly/assembly.fasta \  # 待 polish 的组装结果
    -o medaka_polished/ \           # 输出目录
    -t 16 \                         # 线程数
    -m r1041_e82_400bps_sup_v5.0.0  # 模型（需匹配 basecalling 模型）

# Polish 后的结果：medaka_polished/consensus.fasta
# 白话解释：组装 = 先拼个大概框架，polish = 逐字校对纠错

原始 Subreads（BAM）
   ↓ CCS 生成（pbccs / on-instrument）
HiFi Reads（BAM/FASTQ, Q20+）
   ↓ 比对（minimap2 -x map-hifi / pbmm2）
sorted BAM
   ↓ 变异检测（DeepVariant）
   ↓ 结构变异（pbsv / Sniffles2）
   ↓ 基因组组装（hifiasm）

# ============================================================
# CCS（Circular Consensus Sequencing）：从 subreads 生成 HiFi reads
# 现在大多数情况 Revio 仪器已经在机上自动完成 CCS
# 如果需要手动生成：
# ============================================================

# 安装 pbccs（PacBio SMRT Link 套件的一部分）
# conda install -c bioconda pbccs

# 从 subreads 生成 HiFi reads
ccs \
    subreads.bam \              # 输入：原始 subreads
    hifi_reads.bam \            # 输出：HiFi reads
    --min-rq 0.99 \             # 最低读取质量（0.99 = Q20）
    --min-passes 3 \            # 最少环形读取次数（越多越准，但读长越短）
    --num-threads 32            # 线程数

# 白话解释 CCS 原理：
# DNA 分子被连成一个环，测序仪绕着这个环反复读
# 读 10 圈产生 10 条 subreads，对同一位置取多数表决
# 就像让 10 个人各抄一遍课文，对比取共识 → 极大降低错误率

# 转换为 FASTQ（如果后续工具需要）
samtools fastq hifi_reads.bam > hifi_reads.fastq

# ============================================================
# PacBio HiFi 比对（两种工具可选）
# ============================================================

# 方法一：minimap2（通用，推荐）
minimap2 \
    -ax map-hifi \             # HiFi 专用预设
    -t 16 \
    --secondary=no \
    ref.fasta \
    hifi_reads.fastq \
    | samtools sort -@ 8 -o hifi_aligned.bam

samtools index hifi_aligned.bam

# 方法二：pbmm2（PacBio 官方封装的 minimap2）
# pbmm2 align ref.fasta hifi_reads.bam hifi_aligned.bam \
#     --sort --preset HiFi --num-threads 16

# ============================================================
# DeepVariant：PacBio HiFi 变异检测（Google 深度学习模型）
# ============================================================

# 推荐 Docker 方式运行
docker run \
    -v "$(pwd):/input" \
    -v "$(pwd)/output:/output" \
    google/deepvariant:1.10.0 \
    /opt/deepvariant/bin/run_deepvariant \
    --model_type=PACBIO \                   # 模型类型：PACBIO / ONT / WGS(Illumina)
    --ref=/input/ref.fasta \                # 参考基因组
    --reads=/input/hifi_aligned.bam \       # 输入 BAM
    --output_vcf=/output/deepvariant.vcf.gz \  # 输出 VCF
    --output_gvcf=/output/deepvariant.g.vcf.gz \ # 输出 gVCF（含非变异位点）
    --num_shards=16                         # 并行数（≈线程数）

# DeepVariant 在 PacBio HiFi 数据上的 SNP F1-score 通常 >99.7%
# 这比传统工具（如 GATK HaplotypeCaller）在长读长数据上表现好很多

# ============================================================
# hifiasm：PacBio HiFi 基因组组装（也支持 ONT 超长reads辅助）
# ============================================================

# 安装
conda install -c bioconda hifiasm

# 基本 HiFi 组装
hifiasm \
    -o sample_asm \            # 输出前缀
    -t 32 \                    # 线程数
    hifi_reads.fastq.gz        # 输入 HiFi reads

# 输出文件（GFA 格式，需转换为 FASTA）：
# sample_asm.bp.p_ctg.gfa  → primary assembly（主要组装）
# sample_asm.bp.a_ctg.gfa  → alternate assembly（等位基因版本）

# GFA 转 FASTA
awk '/^S/{print ">"$2;print $3}' sample_asm.bp.p_ctg.gfa > primary_assembly.fasta

# HiFi + ONT 超长 reads 混合组装（获得更好的连续性）
hifiasm \
    -o hybrid_asm \
    -t 32 \
    --ul ont_ultralong.fastq.gz \   # ONT 超长 reads（>50kb）辅助
    hifi_reads.fastq.gz

# HiFi + Hi-C 数据组装（获得染色体级别分型）
hifiasm \
    -o hic_asm \
    -t 32 \
    --h1 hic_R1.fastq.gz \    # Hi-C 正向 reads
    --h2 hic_R2.fastq.gz \    # Hi-C 反向 reads
    hifi_reads.fastq.gz

# hifiasm 特点：
# - 不需要先 polish（HiFi 数据本身就很准）
# - 可以输出 phased（分型的）组装结果
# - 人类基因组组装约 1-2 小时（32线程）

# ============================================================
# Sniffles2：长读长结构变异检测（最主流工具）
# ============================================================

# 安装
conda install -c bioconda sniffles

# 单样本 SV 检测
sniffles \
    --input aligned_sorted.bam \    # 输入排序后的 BAM
    --vcf sv_calls.vcf.gz \         # 输出 VCF
    --reference ref.fasta \         # 参考基因组
    --threads 8 \                   # 线程数
    --minsvlen 50 \                 # 最小 SV 长度（默认 50bp）
    --minsupport auto               # 最小支持 reads 数（auto=自动根据覆盖度决定）

# 多样本联合分析（population-level SV calling）
# 第一步：每个样本先生成 snf 文件
sniffles --input sample1.bam --snf sample1.snf --reference ref.fasta
sniffles --input sample2.bam --snf sample2.snf --reference ref.fasta

# 第二步：合并多个样本
sniffles --input sample1.snf sample2.snf \
    --vcf population_sv.vcf.gz \
    --reference ref.fasta

# ============================================================
# SVIM：另一个长读长 SV 检测工具（特点：灵敏度高）
# ============================================================

# 安装
conda install -c bioconda svim

# SV 检测
svim alignment \
    svim_output/ \                  # 输出目录
    aligned_sorted.bam \            # 输入 BAM
    ref.fasta \                     # 参考基因组
    --min_sv_size 50 \              # 最小 SV 长度
    --minimum_depth 5               # 最小覆盖深度

# 输出：svim_output/variants.vcf（需按质量分过滤）
# 过滤建议：QUAL >= 10 用于 Nanopore，QUAL >= 15 用于 HiFi
bcftools view -i 'QUAL>=10' svim_output/variants.vcf > svim_filtered.vcf

# ============================================================
# Nanopore 直接 RNA 测序分析流程
# ============================================================

# Basecalling（RNA 模式）
dorado basecaller \
    rna004_130bps_sup@v5.1.0 \      # RNA 测序专用模型
    ./pod5_rna/ \                    # RNA 测序产生的 POD5 文件
    --modified-bases m6A inosine \   # 检测 m6A 和肌苷修饰
    --reference transcriptome.fasta \
    > rna_aligned.bam

# RNA 比对（如果不在 basecalling 时比对）
minimap2 \
    -ax splice \                    # splice 模式：考虑 RNA 剪接
    -uf \                           # -u f: 强制正链比对（RNA 有方向性）
    --junc-bed known_junctions.bed \ # 已知剪接位点（可选，提高比对准确率）
    ref.fasta \
    rna_reads.fastq \
    | samtools sort -o rna_sorted.bam

# ============================================================
# modkit：ONT 官方修饰碱基分析工具
# ============================================================

# 从 BAM 文件提取甲基化/修饰信息
modkit pileup \
    rna_aligned.bam \               # 含修饰标签的 BAM
    m6A_pileup.bed \                # 输出 BED 格式
    --ref ref.fasta \               # 参考基因组
    --mod-thresholds m:0.75 \       # 修饰判定阈值（概率 >0.75 认为有修饰）
    --threads 8

# 输出 BED 文件包含每个位置的：
# - 覆盖深度
# - 修饰比例（methylation frequency）
# - 统计置信度

# ============================================================
# 长读长宏基因组组装（以 Nanopore 数据为例）
# ============================================================

# 方法一：Flye 宏基因组模式
flye \
    --nano-hq metagenome_reads.fastq \
    --meta \                        # 关键参数：启用宏基因组模式
    --out-dir meta_assembly/ \      #（处理不均匀覆盖度、多物种混合）
    --threads 64 \
    --genome-size 500m              # 宏基因组预估总大小

# 方法二：metaFlye 直接调用（和上面等价，--meta 即是 metaFlye）

# ============================================================
# 宏基因组 Binning：从组装结果中分离各物种的基因组（MAG）
# ============================================================

# 长读长 MAG 质量通常更高，因为 contig 更长、更连续
# 工具：MetaBAT2 / SemiBin2

# 先计算覆盖度
minimap2 -ax map-ont -t 16 meta_assembly/assembly.fasta \
    metagenome_reads.fastq | samtools sort -o meta_mapped.bam
samtools index meta_mapped.bam

# MetaBAT2 binning
jgi_summarize_bam_contig_depths --outputDepth depth.txt meta_mapped.bam
metabat2 \
    -i meta_assembly/assembly.fasta \   # 组装结果
    -a depth.txt \                      # 覆盖度文件
    -o bins/bin \                       # 输出前缀
    -m 1500 \                           # 最小 contig 长度
    -t 16

# MAG 质量评估
checkm2 predict \
    --input bins/ \                     # binning 结果目录
    --output-directory checkm2_out/ \
    --threads 16

# 面试关联：如果被问"该项目如果用长读长测序会有什么优势"
# → 可以获得完整的微生物基因组（高质量 MAG），
#   更准确地识别 T2D 相关菌株的功能基因和毒力因子