长读测序错误纠正策略¶

# === 方案1：Medaka 纠错 ONT 组装结果（最常用） ===

# 1. 将ONT原始读段比对到组装结果
minimap2 -ax map-ont \                 # ONT比对模式
  assembly.fasta \                     # 组装结果
  ont_reads.fq.gz | \                  # ONT原始读段
  samtools sort -o aligned.bam         # 排序输出BAM
samtools index aligned.bam             # 建立索引

# 2. 用 Medaka 纠错（ONT官方推荐）
medaka_polish \                        # Medaka纠错命令
  -i ont_reads.fq.gz \                # 输入读段
  -d assembly.fasta \                  # 输入组装
  -o medaka_out/ \                     # 输出目录
  -m r1041_e82_400bps_sup_v5.0.0 \    # 模型（匹配测序化学版本）
  -t 16                               # 线程数

# === 方案2：短读 polish（Pilon） ===

# 3. 将短读比对到组装结果
bwa index assembly.fasta               # 建立BWA索引
bwa mem -t 16 \                        # BWA比对
  assembly.fasta \                     # 参考序列
  R1.fq.gz R2.fq.gz | \              # 配对短读
  samtools sort -o short_aligned.bam   # 排序

samtools index short_aligned.bam       # 建立索引

# 4. 运行 Pilon 纠错
java -Xmx64G -jar pilon.jar \         # 分配64GB内存
  --genome assembly.fasta \            # 输入组装
  --frags short_aligned.bam \          # 短读比对结果
  --output pilon_polished \            # 输出前缀
  --changes \                          # 输出修改记录
  --threads 16                         # 线程数

# === 方案3：多轮迭代纠错（推荐流程） ===

# 推荐顺序：长读 polish 2轮 → 短读 polish 2轮
# 第1轮 Medaka
medaka_polish -i reads.fq -d assembly.fa -o round1/ -t 16
# 第2轮 Medaka
medaka_polish -i reads.fq -d round1/consensus.fasta -o round2/ -t 16
# 第3轮 Pilon（用短读）
# ... 同上方案2的流程，输入改为 round2 的结果

# 评估纠错前后的组装质量
import subprocess                      # 系统命令调用

def compare_assemblies(before, after, ref):
    """比较纠错前后的组装质量"""
    # 使用 QUAST 评估
    cmd = (
        f"quast.py "                   # QUAST评估工具
        f"{before} {after} "           # 纠错前和后的组装
        f"-r {ref} "                   # 参考基因组
        f"-o quast_comparison/ "       # 输出目录
        f"--threads 8 "               # 线程数
        f"--labels before,after"       # 标签
    )
    subprocess.run(cmd, shell=True)    # 执行命令
    print("QUAST 报告已生成: quast_comparison/report.html")

# 统计 Pilon 修改数量
def count_pilon_changes(changes_file):
    """统计 Pilon 做了多少修改"""
    snps = 0                           # SNP修正数
    insertions = 0                     # 插入修正数
    deletions = 0                      # 删除修正数
    with open(changes_file) as f:      # 打开变更文件
        for line in f:                 # 逐行读取
            if "snp" in line.lower():  # SNP类型
                snps += 1
            elif "ins" in line.lower():  # 插入类型
                insertions += 1
            elif "del" in line.lower():  # 删除类型
                deletions += 1
    print(f"SNP修正: {snps}")          # 打印结果
    print(f"插入修正: {insertions}")
    print(f"删除修正: {deletions}")
    print(f"总修正: {snps + insertions + deletions}")