ONT 自适应采样技术¶

# ===== ONT 自适应采样配置与分析 =====

# 方法1：使用 MinKNOW 内置的 Adaptive Sampling
# 在 MinKNOW 软件中配置（GUI 操作）：
# 1. 开始新测序实验
# 2. 在 "Adaptive Sampling" 选项中上传 BED 文件
# 3. 选择 "Enrich"（富集目标）或 "Deplete"（去除非目标）
# 4. 提供参考基因组 FASTA 文件

# 准备目标区域 BED 文件（示例：BRCA1 + BRCA2 区域）
cat > target_regions.bed << 'EOF'
chr17   43044295    43170245    BRCA1
chr13   32315474    32400266    BRCA2
EOF

# 方法2：使用 Readfish 进行自适应采样（更灵活）
# Readfish 是 ONT 官方支持的自适应采样工具

# 安装 readfish
pip install readfish

# 准备 TOML 配置文件
cat > adaptive_config.toml << 'TOML'
[caller_settings]
config_name = "dna_r10.4.1_e8.2_400bps_sup.cfg"
host = "localhost"
port = 8000

[conditions.target_regions]
name = "target"
control = false
min_chunks = 1
max_chunks = 4
targets = "target_regions.bed"
single_on = "stop_receiving"    # 目标序列：继续读
single_off = "unblock"          # 非目标：弹出
multi_on = "stop_receiving"
multi_off = "unblock"
no_seq = "proceed"
no_map = "proceed"
TOML

# 运行自适应采样（需要与 MinKNOW 同时运行）
readfish targets \
    --toml adaptive_config.toml \
    --device MN12345 \    # 设备序列号
    --experiment "adaptive_brca"

# ===== 数据分析 =====

# 1. Basecalling（如果未实时完成）
dorado basecaller sup \
    pod5_input/ \
    --reference hg38.fa \
    --modified-bases 5mCG_5hmCG > calls.bam  # 同时检测甲基化

# 2. 比对到参考基因组
minimap2 -a -x map-ont -t 16 \
    hg38.fa calls.fastq.gz | \
    samtools sort -@ 8 -o aligned.bam
samtools index aligned.bam

# 3. 统计富集效果
# 计算目标区域的覆盖度
samtools depth -b target_regions.bed aligned.bam | \
    awk '{sum+=$3; n++} END {print "目标区域平均覆盖度:", sum/n, "x"}'

# 计算全基因组平均覆盖度
samtools depth aligned.bam | \
    awk '{sum+=$3; n++} END {print "全基因组平均覆盖度:", sum/n, "x"}'

# 富集倍数 = 目标区域覆盖度 / 全基因组覆盖度

# 4. 统计 accept/reject reads
# 检查 BAM 文件中的标签
python3 << 'PYEOF'
import pysam

bam = pysam.AlignmentFile("aligned.bam", "rb")
on_target = 0   # 目标区域的 reads
off_target = 0  # 非目标区域的 reads
rejected = 0    # 被弹出的短 reads（<1kb 可能是 rejected）

for read in bam:
    if read.query_length < 1000:
        rejected += 1
    elif read.reference_name in ["chr17", "chr13"]:
        # 简化判断：实际应检查精确坐标
        on_target += 1
    else:
        off_target += 1

total = on_target + off_target + rejected
print(f"目标区域 reads: {on_target} ({on_target/total*100:.1f}%)")
print(f"非目标 reads: {off_target} ({off_target/total*100:.1f}%)")
print(f"Rejected (短) reads: {rejected} ({rejected/total*100:.1f}%)")
PYEOF