Nanopore 16S 全长测序分析¶

# ===== Nanopore 16S 全长分析流程 =====

# 1. Basecalling（如果使用 live basecalling 则跳过）
dorado basecaller sup \
    pod5_directory/ \
    --kit-name SQK-16S114.24 > basecalled.bam

# 转换为 FASTQ
samtools fastq basecalled.bam > raw_reads.fastq

# 2. 质控：过滤低质量和异常长度的读
# Nanopore 16S 应该在 1300-1600bp
chopper --quality 10 \
    --minlength 1300 \
    --maxlength 1600 \
    --input raw_reads.fastq > filtered_reads.fastq

# 统计过滤后的数据
echo "过滤后的 reads 数:"
grep -c "^@" filtered_reads.fastq
echo "reads 长度分布:"
awk 'NR%4==2 {print length($0)}' filtered_reads.fastq | \
    sort -n | \
    awk 'BEGIN{print "Min\tMedian\tMax"} {a[NR]=$0} END{print a[1]"\t"a[int(NR/2)]"\t"a[NR]}'

# 3. 方法1：EMU — 种水平物种丰度估计（推荐）
# EMU 使用期望最大化算法处理 Nanopore 的错误率
emu abundance \
    --db /ref/emu_database/ \
    --type map-ont \
    --threads 16 \
    --output-dir emu_output/ \
    filtered_reads.fastq

# EMU 输出：每个样本的种水平丰度表
# species_name  abundance  tax_id
cat emu_output/filtered_reads_rel-abundance.tsv

# 4. 方法2：minimap2 比对到 SILVA 数据库
minimap2 -a -x map-ont \
    --secondary=no \
    -t 16 \
    /ref/SILVA_138.2_SSU_NR99.fasta \
    filtered_reads.fastq | \
    samtools sort -@ 8 -o silva_aligned.bam
samtools index silva_aligned.bam

# 统计比对到各物种的 reads 数
samtools idxstats silva_aligned.bam | \
    sort -k3 -rn | head -20  # 前20个最丰富的物种

# 5. 方法3：使用 QIIME2 + q2-long-read 插件
# 导入数据到 QIIME2
qiime tools import \
    --type 'SampleData[SequencesWithQuality]' \
    --input-path manifest.csv \
    --output-path reads.qza \
    --input-format SingleEndFastqManifestPhred33V2

# 使用 VSEARCH 去嵌合体和聚类
qiime vsearch cluster-features-de-novo \
    --i-sequences reads.qza \
    --i-table feature-table.qza \
    --p-perc-identity 0.99 \
    --o-clustered-table clustered-table.qza \
    --o-clustered-sequences clustered-seqs.qza

# 分类注释（使用 SILVA 分类器）
qiime feature-classifier classify-consensus-vsearch \
    --i-query clustered-seqs.qza \
    --i-reference-reads silva-138-ssu-nr99-seqs.qza \
    --i-reference-taxonomy silva-138-ssu-nr99-tax.qza \
    --o-classification taxonomy.qza \
    --p-threads 16

# ===== R: 全长 16S 结果可视化与比较 =====
library(ggplot2)
library(dplyr)

# 读取 EMU 输出的丰度表
emu_results <- read.delim("emu_output/filtered_reads_rel-abundance.tsv")

# 按丰度排序，取 Top 15 物种
top15 <- emu_results %>%
    arrange(desc(abundance)) %>%
    head(15)

# 物种丰度条形图
ggplot(top15, aes(x = reorder(species, abundance), y = abundance)) +
    geom_bar(stat = "identity", fill = "steelblue") +
    coord_flip() +   # 横向展示，方便看物种名
    labs(x = "物种", y = "相对丰度 (%)",
         title = "Nanopore 16S 全长：Top 15 物种") +
    theme_minimal()

# 比较全长 vs V3-V4 的分辨率
# 全长结果
full_length <- data.frame(
    level = c("门", "纲", "目", "科", "属", "种"),
    resolution = c(99, 98, 97, 95, 92, 85))  # 示例分辨率

# V3-V4 结果
v3v4 <- data.frame(
    level = c("门", "纲", "目", "科", "属", "种"),
    resolution = c(99, 97, 95, 90, 80, 50))  # 种水平分辨率低

comparison <- rbind(
    cbind(full_length, method = "全长 16S"),
    cbind(v3v4, method = "V3-V4"))

ggplot(comparison, aes(x = level, y = resolution,
                       color = method, group = method)) +
    geom_line(linewidth = 1.2) +
    geom_point(size = 3) +
    scale_x_discrete(limits = c("门","纲","目","科","属","种")) +
    labs(x = "分类水平", y = "分辨率 (%)",
         title = "全长 16S vs V3-V4：分类分辨率比较") +
    theme_minimal()