以下是您提供内容的完整中文翻译，已针对生物信息学与表观遗传学语境进行专业术语校准，并保留原有排版结构：

---

✅ 关键说明：您的样本是细菌！
这极大地改变了结果的解释方向——实际上让您的发现变得更加有趣，而非毫无意义。

🔬 更新后的生物学背景

维度	先前假设（植物）	正确背景（细菌）
REBASE 相关性	低（原核数据库 vs. 真核数据）	✅ 高（原核数据库 vs. 原核数据）
预期匹配率	5-15%（偶然匹配）	若甲基化由 R-M 系统驱动，匹配率应更高
甲基化类型	植物中以 5mC 为主	✅ 细菌中 5mC、4mC 和 6mA 均十分常见
生物学驱动因素	转录因子结合影响甲基化	✅ 限制-修饰系统 或 孤儿甲基化酶

---

📊 重新解读您的结果（细菌背景）

5mC_CG:  18 motifs → 2 REBASE matches (11.1%) → 1 methylase match ⚠️
4mC_CG:  23 motifs → 3 REBASE matches (13.0%) → 0 methylase matches
5mC_nonCG: 15 motifs → 1 REBASE match (6.7%) → 0 methylase matches

为什么匹配率仍然适中（约 10%）？

1. HOMER 识别的是类似转录因子结合的基序：即使在细菌中，DNA结合蛋白（不仅是R-M系统）也会造成序列富集。
2. 孤儿甲基化酶（Orphan methylases）：许多细菌甲基转移酶缺乏对应的限制性内切酶，且可能尚未被完全收录至数据库。
3. 菌株特异性新发现：Acinetobacter harbinensis 和 Pedobacter sp. BG5 可能拥有尚未收录到 REBASE 中的新型 R-M 系统。
4. 序列上下文很重要：您的基序在甲基化位点处富集，但驱动甲基化的序列特异性可能比催化核心序列更宽泛。

---

🎯 5mC_CG 数据中唯一的甲基化酶匹配具有高度重要性

由于您的研究对象是细菌，REBASE 中的甲基化酶匹配提示：

该基序处的甲基化信号可能直接由已知的细菌甲基转移酶催化。

提取并验证该匹配：

# 找出匹配甲基化酶的基序
grep -v "^motif" 5mC_CG_results.tsv | awk -F'\t' '$6 > 0 {print}'

随后在 REBASE 中核实：

# 在 withrefm.txt 中搜索该酶
grep -B2 -A10 "
<1>M\.酶名称" withrefm.txt | grep -E "(<1>|<3>|<4>)"

**在 `

3>` 或 “ 字段中查找甲基化标注：** “` GATATC2(6) → M.EcoRV 在第2位腺嘌呤上添加 N6-甲基 (6mA) CCWGG3(5) → 在第3位胞嘧啶上添加 5-甲基 (5mC) GANTC2(4) → 在第2位胞嘧啶上添加 N4-甲基 (4mC) “` ⚠️ **重要提示**：如果匹配到的是 `6mA` 甲基化酶（标注为 `(6)`），但您研究的是 `5mC/4mC`，则可能表明： – 检测交叉反应或实验假象 – 同一基因组位点上多种甲基化类型共存 – 不同甲基化系统之间的生物学串扰（Crosstalk） — ### 🔍 推荐的后续步骤 **1. 提取所有匹配的基序及其详细信息：** “`python # quick_extract_matches.py import csv with open(‘5mC_CG_results.tsv’, ‘r’) as f: reader = csv.DictReader(f, delimiter=’\t’) for row in reader: if int(row[‘n_rebase_matches’]) > 0: print(f”\n基序: {row[‘motif’]}”) print(f” HOMER最佳匹配TF: {row[‘bestguess_tf’]}”) print(f” REBASE匹配酶: {row[‘matched_enzymes’]}”) if row[‘methylase_matches’]: print(f” ⚠️ 甲基化酶: {row[‘methylase_matches’]}”) print(f” 备注: {row[‘notes’]}”) “` **2. 针对每个甲基化酶匹配，核查其甲基化类型：** “`bash # 批量查看匹配酶的修饰位点标注 for enzyme in $(grep -v “^motif” 5mC_CG_results.tsv | awk -F’\t’ ‘$6 > 0 {print $8}’ | tr ‘;’ ‘\n’ | sort -u); do echo “=== $enzyme ===” grep ” $enzyme” withrefm.txt | grep -oP ‘(?<=)[^<]+|(?<=)[^<]+' | head -2 done “` **3. 考虑扩展至 6mA 分析：** 由于细菌基因组中普遍存在 N6-甲基腺嘌呤，若您有 6mA 测序数据，建议同步运行： “`bash python annotate_motifs_rebase.py \ –homer 6mA_CG/nonRedundant.motifs \ –rebase withrefm.txt \ –output 6mA_CG_results.tsv “` — ### 🧠 生物学解读框架（细菌） | 观察结果 | 可能的生物学解释 | |———-|—————-| | **基序仅匹配限制性内切酶** | 偶然的序列相似性；甲基化可能来自孤儿甲基化酶或其他系统 | | **基序匹配标注为 `3(5)` 的甲基化酶** | ✅ 第3位胞嘧啶发生 **5-甲基化** —— 与您的 5mC 研究直接相关 | | **基序匹配标注为 `5(4)` 的甲基化酶** | ✅ 第5位胞嘧啶发生 **N4-甲基化** —— 与您的 4mC 研究直接相关 | | **基序匹配标注为 `(6)` 的甲基化酶** | **N6-甲基腺嘌呤** —— 可能提示 6mA 串扰或检测交叉反应 | | **无 REBASE 匹配** | 可能是细菌转录调控蛋白的结合位点；甲基化可能是次要效应或由新型未知酶催化 | — ### ✅ 快速验证清单 – [ ] 运行 `quick_extract_matches.py` 识别 `5mC_CG` 中匹配甲基化酶的基序 – [ ] 在 `withrefm.txt` 中核实甲基化类型标注（`(4)`、`(5)` 或 `(6)`） – [ ] 交叉比对匹配酶的来源菌株是否与 *Acinetobacter* 或 *Pedobacter* 存在亲缘关系 – [ ] 评估短基序（如 `TCGA`，4bp）的匹配是否具有生物学意义，还是随机假阳性 – [ ] 对于未匹配的基序：可在 JASPAR/RegulonDB 中检索细菌转录因子结合位点 — ### 🎯 核心结论 您的结果在细菌背景下具有生物学合理性且可被清晰解读： – 约 **10% 的 REBASE 匹配率** 表明，大多数富集基序反映的是细菌 DNA 结合蛋白的特异性，而不仅仅是 R-M 系统。 – `5mC_CG` 中**唯一的甲基化酶匹配值得重点验证**——它可能直接定位负责您甲基化信号的关键酶。 – `4mC/nonCG` 数据集中匹配率较低，可能真实反映了甲基化机制的生物学差异。 需要我帮您从 `5mC_CG_results.tsv` 中提取并详细解读那个特定的甲基化酶匹配基序吗？🎯

🚀 DNA甲基化分析一体化脚本

这是一个完整的”过滤 + 分离 + 分析”自动化脚本：

#!/bin/bash
#===============================================================================
#  脚本名称: dna_methylation_analysis.sh
#  功能描述: DNA甲基化数据过滤、分离、motif分析和可视化一体化流程
#  作者: Bioinformatics Pipeline
#  版本: 1.0
#  日期: 2026
#  用法: bash dna_methylation_analysis.sh -i input.bed -r reference.fa -o output_dir
#===============================================================================

set -euo pipefail  # 严格模式：遇到错误立即退出

#-------------------------------------------------------------------------------
#  配置参数
#-------------------------------------------------------------------------------

# 默认参数
INPUT_BED=""
REF_GENOME=""
OUTPUT_DIR="methylation_analysis"
MIN_COVERAGE=10
MIN_MOD_5MC=70
MIN_MOD_4MC=50
FLANK_SIZE=50
THREADS=8
RUN_HOMER=true
RUN_MEME=false
RUN_VISUALIZATION=true

# 颜色定义
RED='\033[0;31m'
GREEN='\033[0;32m'
YELLOW='\033[1;33m'
BLUE='\033[0;34m'
NC='\033[0m' # No Color

# 日志函数
log_info() {
    echo -e "${GREEN}[INFO]${NC} $(date '+%Y-%m-%d %H:%M:%S') - $1"
}

log_warn() {
    echo -e "${YELLOW}[WARN]${NC} $(date '+%Y-%m-%d %H:%M:%S') - $1"
}

log_error() {
    echo -e "${RED}[ERROR]${NC} $(date '+%Y-%m-%d %H:%M:%S') - $1"
    exit 1
}

log_step() {
    echo -e "${BLUE}[STEP]${NC} $(date '+%Y-%m-%d %H:%M:%S') - $1"
    echo "━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━"
}

#-------------------------------------------------------------------------------
#  参数解析
#-------------------------------------------------------------------------------

usage() {
    cat << EOF
用法: $0 -i <input.bed> -r <reference.fa> [选项]

必需参数:
  -i, --input FILE        输入的bed文件 (modkit pileup输出)
  -r, --reference FILE    参考基因组fasta文件

可选参数:
  -o, --outdir DIR        输出目录 (默认: methylation_analysis)
  -c, --coverage INT      最小覆盖度 (默认: 10)
  -m, --mod5mc INT        5mC最小修饰率% (默认: 70)
  -n, --mod4mc INT        4mC最小修饰率% (默认: 50)
  -f, --flank INT         提取序列的侧翼大小 (默认: 50)
  -t, --threads INT       并行线程数 (默认: 8)
  --no-homer              不运行HOMER motif分析
  --meme                  运行MEME motif分析 (默认关闭)
  --no-viz                不生成可视化图表
  -h, --help              显示帮助信息

示例:
  $0 -i An6_4mC_5mC_raw.bed -r genome.fa -o results -c 10 -m 70
  $0 -i input.bed -r genome.fa --mod5mc 80 --mod4mc 60 --threads 16

EOF
    exit 0
}

# 解析命令行参数
while [[ $# -gt 0 ]]; do
    case $1 in
        -i|--input)
            INPUT_BED="$2"
            shift 2
            ;;
        -r|--reference)
            REF_GENOME="$2"
            shift 2
            ;;
        -o|--outdir)
            OUTPUT_DIR="$2"
            shift 2
            ;;
        -c|--coverage)
            MIN_COVERAGE="$2"
            shift 2
            ;;
        -m|--mod5mc)
            MIN_MOD_5MC="$2"
            shift 2
            ;;
        -n|--mod4mc)
            MIN_MOD_4MC="$2"
            shift 2
            ;;
        -f|--flank)
            FLANK_SIZE="$2"
            shift 2
            ;;
        -t|--threads)
            THREADS="$2"
            shift 2
            ;;
        --no-homer)
            RUN_HOMER=false
            shift
            ;;
        --meme)
            RUN_MEME=true
            shift
            ;;
        --no-viz)
            RUN_VISUALIZATION=false
            shift
            ;;
        -h|--help)
            usage
            ;;
        *)
            log_error "未知参数: $1"
            usage
            ;;
    esac
done

#-------------------------------------------------------------------------------
#  参数验证
#-------------------------------------------------------------------------------

log_step "验证输入参数"

[[ -z "$INPUT_BED" ]] && log_error "未指定输入bed文件 (-i)"
[[ -z "$REF_GENOME" ]] && log_error "未指定参考基因组 (-r)"
[[ ! -f "$INPUT_BED" ]] && log_error "输入文件不存在: $INPUT_BED"
[[ ! -f "$REF_GENOME" ]] && log_error "参考基因组不存在: $REF_GENOME"

# 检查必需工具
check_tool() {
    if ! command -v "$1" &> /dev/null; then
        log_warn "未找到工具: $1 (某些功能可能不可用)"
        return 1
    fi
    return 0
}

log_info "检查必需工具..."
check_tool samtools || SAMTOOLS_FOUND=false
check_tool bedtools || BEDTOOLS_FOUND=false
check_tool awk || AWK_FOUND=false

if [[ "$RUN_HOMER" == "true" ]]; then
    check_tool findMotifsGenome.pl || HOMER_FOUND=false
fi

#-------------------------------------------------------------------------------
#  创建输出目录结构
#-------------------------------------------------------------------------------

log_step "创建输出目录"

mkdir -p "$OUTPUT_DIR"/{filtered,sequences,motif_analysis/homer,motif_analysis/meme,figures,reports,logs}

# 记录配置
cat > "$OUTPUT_DIR/reports/analysis_config.txt" << EOF
=====================================
DNA甲基化分析配置
=====================================
分析时间: $(date)
输入文件: $INPUT_BED
参考基因组: $REF_GENOME
输出目录: $OUTPUT_DIR

过滤参数:
  最小覆盖度: $MIN_COVERAGE
  5mC最小修饰率: ${MIN_MOD_5MC}%
  4mC最小修饰率: ${MIN_MOD_4MC}%
  侧翼序列大小: ${FLANK_SIZE}bp

分析选项:
  运行HOMER: $RUN_HOMER
  运行MEME: $RUN_MEME
  生成可视化: $RUN_VISUALIZATION
  线程数: $THREADS
=====================================
EOF

log_info "配置已保存到: $OUTPUT_DIR/reports/analysis_config.txt"

#-------------------------------------------------------------------------------
#  步骤1: 数据质量评估
#-------------------------------------------------------------------------------

log_step "步骤1/7: 数据质量评估"

RAW_COUNT=$(wc -l < "$INPUT_BED")
log_info "原始数据总行数: $RAW_COUNT"

# 统计各类型位点数量
log_info "统计各修饰类型..."

COUNT_5MC_CG=$(awk -F'\t' '$4 ~ /^m,CG,/' "$INPUT_BED" | wc -l)
COUNT_5MC_NONCG=$(awk -F'\t' '$4 ~ /^m,C,/ && $4 !~ /^m,CG,/' "$INPUT_BED" | wc -l)
COUNT_4MC_CG=$(awk -F'\t' '$4 ~ /^21839,CG,/' "$INPUT_BED" | wc -l)
COUNT_4MC_NONCG=$(awk -F'\t' '$4 ~ /^21839,C,/ && $4 !~ /^21839,CG,/' "$INPUT_BED" | wc -l)

cat > "$OUTPUT_DIR/reports/raw_statistics.txt" << EOF
原始数据统计 (无过滤)
=====================================
5mC @ CG:      $COUNT_5MC_CG
5mC @ non-CG:  $COUNT_5MC_NONCG
4mC @ CG:      $COUNT_4MC_CG
4mC @ non-CG:  $COUNT_4MC_NONCG
-------------------------------------
总计:          $RAW_COUNT
=====================================
EOF

# 修饰率分布
log_info "计算修饰率分布..."
awk -F'\t' '{print $11}' "$INPUT_BED" | \
    awk '{
        if($1>=90) h["90-100"]++;
        else if($1>=70) h["70-90"]++;
        else if($1>=50) h["50-70"]++;
        else if($1>=30) h["30-50"]++;
        else h["0-30"]++
    } END {
        print "修饰率分布:";
        for(i in h) print i"%: "h[i]
    }' > "$OUTPUT_DIR/reports/mod_rate_distribution.txt"

log_info "质量评估完成！"

#-------------------------------------------------------------------------------
#  步骤2: 数据过滤和分离
#-------------------------------------------------------------------------------

log_step "步骤2/7: 数据过滤和分离"

# 定义过滤函数
filter_and_separate() {
    local type_name=$1
    local pattern=$2
    local min_mod=$3
    local output_file="$OUTPUT_DIR/filtered/${type_name}_filtered.bed"

    log_info "过滤 $type_name (修饰率≥${min_mod}%, 覆盖度≥$MIN_COVERAGE)..."

    awk -F'\t' -v pat="$pattern" -v min_mod="$min_mod" -v min_cov="$MIN_COVERAGE" '
        $4 ~ pat && $11 >= min_mod && $5 >= min_cov {
            print $0
        }' "$INPUT_BED" > "$output_file"

    local count=$(wc -l < "$output_file")
    log_info "  → 保留 $count 个位点"
    echo "$count"
}

# 过滤和分离4种类型
COUNT_FILT_5MC_CG=$(filter_and_separate "5mC_CG" "^m,CG," "$MIN_MOD_5MC")
COUNT_FILT_5MC_NONCG=$(filter_and_separate "5mC_nonCG" "^m,C," "$MIN_MOD_5MC" | \
    awk -v pat="^m,CG," '$4 !~ pat' > "$OUTPUT_DIR/filtered/5mC_nonCG_filtered.bed" && \
    wc -l < "$OUTPUT_DIR/filtered/5mC_nonCG_filtered.bed")

# 重新正确过滤5mC non-CG（排除CG）
awk -F'\t' -v min_mod="$MIN_MOD_5MC" -v min_cov="$MIN_COVERAGE" '
    $4 ~ /^m,C,/ && $4 !~ /^m,CG,/ && $11 >= min_mod && $5 >= min_cov {
        print $0
    }' "$INPUT_BED" > "$OUTPUT_DIR/filtered/5mC_nonCG_filtered.bed"
COUNT_FILT_5MC_NONCG=$(wc -l < "$OUTPUT_DIR/filtered/5mC_nonCG_filtered.bed")

COUNT_FILT_4MC_CG=$(filter_and_separate "4mC_CG" "^21839,CG," "$MIN_MOD_4MC")
COUNT_FILT_4MC_NONCG=$(filter_and_separate "4mC_nonCG" "^21839,C," "$MIN_MOD_4MC" | \
    awk -v pat="^21839,CG," '$4 !~ pat' > "$OUTPUT_DIR/filtered/4mC_nonCG_filtered.bed" && \
    wc -l < "$OUTPUT_DIR/filtered/4mC_nonCG_filtered.bed")

# 重新正确过滤4mC non-CG
awk -F'\t' -v min_mod="$MIN_MOD_4MC" -v min_cov="$MIN_COVERAGE" '
    $4 ~ /^21839,C,/ && $4 !~ /^21839,CG,/ && $11 >= min_mod && $5 >= min_cov {
        print $0
    }' "$INPUT_BED" > "$OUTPUT_DIR/filtered/4mC_nonCG_filtered.bed"
COUNT_FILT_4MC_NONCG=$(wc -l < "$OUTPUT_DIR/filtered/4mC_nonCG_filtered.bed")

# 生成过滤统计报告
cat > "$OUTPUT_DIR/reports/filtering_summary.txt" << EOF
过滤结果统计
=====================================
过滤条件:
  最小覆盖度: $MIN_COVERAGE
  5mC修饰率阈值: ${MIN_MOD_5MC}%
  4mC修饰率阈值: ${MIN_MOD_4MC}%

过滤后位点数量:
  5mC @ CG:      $COUNT_FILT_5MC_CG (原始: $COUNT_5MC_CG, 保留: $(awk "BEGIN {printf \"%.1f\", $COUNT_FILT_5MC_CG/$COUNT_5MC_CG*100}")%)
  5mC @ non-CG:  $COUNT_FILT_5MC_NONCG (原始: $COUNT_5MC_NONCG, 保留: $(awk "BEGIN {printf \"%.1f\", $COUNT_FILT_5MC_NONCG/$COUNT_5MC_NONCG*100}")%)
  4mC @ CG:      $COUNT_FILT_4MC_CG (原始: $COUNT_4MC_CG, 保留: $(awk "BEGIN {printf \"%.1f\", $COUNT_FILT_4MC_CG/$COUNT_4MC_CG*100}")%)
  4mC @ non-CG:  $COUNT_FILT_4MC_NONCG (原始: $COUNT_4MC_NONCG, 保留: $(awk "BEGIN {printf \"%.1f\", $COUNT_FILT_4MC_NONCG/$COUNT_4MC_NONCG*100}")%)
=====================================
EOF

log_info "过滤统计已保存到: $OUTPUT_DIR/reports/filtering_summary.txt"
cat "$OUTPUT_DIR/reports/filtering_summary.txt"

#-------------------------------------------------------------------------------
#  步骤3: 提取序列上下文
#-------------------------------------------------------------------------------

log_step "步骤3/7: 提取序列上下文"

# 检查bedtools
if [[ "${BEDTOOLS_FOUND:-true}" != "true" ]]; then
    log_warn "bedtools未安装，跳过序列提取"
    SEQUENCES_EXTRACTED=false
else
    SEQUENCES_EXTRACTED=true

    # 索引参考基因组
    if [[ ! -f "${REF_GENOME}.fai" ]]; then
        log_info "索引参考基因组..."
        samtools faidx "$REF_GENOME"
    fi

    # 创建染色体大小文件
    cut -f1,2 "${REF_GENOME}.fai" > "$OUTPUT_DIR/reports/chrom.sizes"

    # 提取序列函数
    extract_sequences() {
        local bed_file=$1
        local type_name=$2
        local output_fa="$OUTPUT_DIR/sequences/${type_name}_sequences.fa"

        if [[ ! -s "$bed_file" ]]; then
            log_warn "$type_name 位点数为0，跳过序列提取"
            return
        fi

        log_info "提取 $type_name 序列 (±${FLANK_SIZE}bp)..."

        # 生成BED文件（带侧翼序列）
        awk -F'\t' -v flank="$FLANK_SIZE" '{
            start = ($2 - flank > 0) ? $2 - flank : 0;
            end = $3 + flank;
            print $1"\t"start"\t"end"\t"$4"_"NR
        }' "$bed_file" > "$OUTPUT_DIR/sequences/${type_name}_regions.bed"

        # 提取fasta
        bedtools getfasta -fi "$REF_GENOME" \
                          -bed "$OUTPUT_DIR/sequences/${type_name}_regions.bed" \
                          -fo "$output_fa" \
                          -nameOnly

        local seq_count=$(grep -c "^>" "$output_fa" || echo 0)
        log_info "  → 提取 $seq_count 条序列"
    }

    # 为每个类型提取序列
    extract_sequences "$OUTPUT_DIR/filtered/5mC_CG_filtered.bed" "5mC_CG"
    extract_sequences "$OUTPUT_DIR/filtered/5mC_nonCG_filtered.bed" "5mC_nonCG"
    extract_sequences "$OUTPUT_DIR/filtered/4mC_CG_filtered.bed" "4mC_CG"
    extract_sequences "$OUTPUT_DIR/filtered/4mC_nonCG_filtered.bed" "4mC_nonCG"
fi

#-------------------------------------------------------------------------------
#  步骤4: 生成背景序列
#-------------------------------------------------------------------------------

log_step "步骤4/7: 生成背景序列"

if [[ "$SEQUENCES_EXTRACTED" == "true" ]]; then
    # 合并所有target区域用于生成背景
    cat "$OUTPUT_DIR/sequences/"*_regions.bed 2>/dev/null | \
        bedtools shuffle -g "$OUTPUT_DIR/reports/chrom.sizes" | \
        bedtools getfasta -fi "$REF_GENOME" -bed stdin \
        -fo "$OUTPUT_DIR/sequences/background_sequences.fa" \
        -nameOnly

    BG_COUNT=$(grep -c "^>" "$OUTPUT_DIR/sequences/background_sequences.fa" || echo 0)
    log_info "生成 $BG_COUNT 条背景序列"
else
    log_warn "跳过背景序列生成"
fi

#-------------------------------------------------------------------------------
#  步骤5: Motif富集分析 (HOMER)
#-------------------------------------------------------------------------------

log_step "步骤5/7: Motif富集分析 (HOMER)"

if [[ "$RUN_HOMER" == "true" ]]; then
    if [[ "${HOMER_FOUND:-true}" != "true" ]]; then
        log_warn "HOMER未安装，跳过motif分析"
        RUN_HOMER=false
    else
        run_homer_motif() {
            local type_name=$1
            local bed_file="$OUTPUT_DIR/sequences/${type_name}_regions.bed"
            local output_dir="$OUTPUT_DIR/motif_analysis/homer/${type_name}"

            if [[ ! -s "$bed_file" ]]; then
                log_warn "$type_name 无位点，跳过HOMER分析"
                return
            fi

            local region_count=$(wc -l < "$bed_file")
            if [[ $region_count -lt 10 ]]; then
                log_warn "$type_name 位点数过少 ($region_count)，跳过HOMER分析"
                return
            fi

            log_info "运行HOMER: $type_name..."
            mkdir -p "$output_dir"

            findMotifsGenome.pl "$bed_file" \
                "$REF_GENOME" \
                "$output_dir" \
                -bg "$OUTPUT_DIR/sequences/background_sequences.fa" \
                -len 6,8,10,12 \
                -p "$THREADS" \
                > "$output_dir/homer.log" 2>&1

            log_info "  → 结果保存至: $output_dir"
        }

        run_homer_motif "5mC_CG"
        run_homer_motif "5mC_nonCG"
        run_homer_motif "4mC_CG"
        run_homer_motif "4mC_nonCG"

        # 汇总HOMER结果
        log_info "汇总HOMER结果..."
        echo "HOMER Motif富集结果摘要" > "$OUTPUT_DIR/reports/homer_summary.txt"
        echo "=====================================" >> "$OUTPUT_DIR/reports/homer_summary.txt"

        for type in 5mC_CG 5mC_nonCG 4mC_CG 4mC_nonCG; do
            results_file="$OUTPUT_DIR/motif_analysis/homer/${type}/knownResults.txt"
            if [[ -f "$results_file" ]]; then
                echo -e "\n$type 前5个富集motif:" >> "$OUTPUT_DIR/reports/homer_summary.txt"
                head -6 "$results_file" >> "$OUTPUT_DIR/reports/homer_summary.txt"
            fi
        done
    fi
else
    log_info "HOMER分析已禁用"
fi

#-------------------------------------------------------------------------------
#  步骤6: Motif富集分析 (MEME - 可选)
#-------------------------------------------------------------------------------

log_step "步骤6/7: Motif富集分析 (MEME)"

if [[ "$RUN_MEME" == "true" ]]; then
    if ! command -v dreme &> /dev/null; then
        log_warn "MEME Suite未安装，跳过MEME分析"
        RUN_MEME=false
    else
        run_meme_dreme() {
            local type_name=$1
            local seq_file="$OUTPUT_DIR/sequences/${type_name}_sequences.fa"
            local output_dir="$OUTPUT_DIR/motif_analysis/meme/${type_name}"

            if [[ ! -s "$seq_file" ]]; then
                return
            fi

            local seq_count=$(grep -c "^>" "$seq_file" || echo 0)
            if [[ $seq_count -lt 10 ]]; then
                return
            fi

            log_info "运行DREME: $type_name..."
            mkdir -p "$output_dir"

            dreme -p "$seq_file" \
                  -n "$OUTPUT_DIR/sequences/background_sequences.fa" \
                  -oc "$output_dir" \
                  -e 0.05 \
                  -p "$THREADS" \
                  > "$output_dir/dreme.log" 2>&1

            log_info "  → 结果保存至: $output_dir"
        }

        run_meme_dreme "5mC_CG"
        run_meme_dreme "5mC_nonCG"
        run_meme_dreme "4mC_CG"
        run_meme_dreme "4mC_nonCG"
    fi
else
    log_info "MEME分析已禁用"
fi

#-------------------------------------------------------------------------------
#  步骤7: 生成可视化图表
#-------------------------------------------------------------------------------

log_step "步骤7/7: 生成可视化图表"

if [[ "$RUN_VISUALIZATION" == "true" ]]; then
    # 创建Python可视化脚本
    cat > "$OUTPUT_DIR/generate_plots.py" << 'PYTHON_SCRIPT'
#!/usr/bin/env python3
"""DNA甲基化数据可视化脚本"""

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import os
import sys

# 设置风格
sns.set_style("whitegrid")
sns.set_context("paper", font_scale=1.2)
plt.rcParams['font.sans-serif'] = ['Arial', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

def load_data(bed_file):
    """加载bed文件"""
    columns = ['chrom', 'start', 'end', 'name', 'score', 'strand', 
               'thickStart', 'thickEnd', 'color', 'Nvalid_cov', 
               'percent_modified', 'Nmod', 'Ncanonical', 'Nother_mod',
               'Ndelete', 'Nfail', 'Ndiff', 'Nnocall']
    return pd.read_csv(bed_file, sep='\t', header=None, names=columns)

def classify_modifications(df):
    """分类修饰类型"""
    df = df.copy()
    df['mod_type'] = df['name'].apply(lambda x: '5mC' if x.startswith('m,') else '4mC')
    df['motif_type'] = df['name'].apply(lambda x: 'CG' if ',CG,' in x else 'non-CG')
    df['category'] = df['mod_type'] + '_' + df['motif_type']
    return df

def plot_counts_comparison(df_raw, df_filt, output_file):
    """图1: 位点数量对比"""
    fig, ax = plt.subplots(figsize=(10, 6))

    categories = ['5mC_CG', '5mC_nonCG', '4mC_CG', '4mC_nonCG']
    counts_raw = [len(df_raw[df_raw['category']==c]) for c in categories]
    counts_filt = [len(df_filt[df_filt['category']==c]) for c in categories]

    x = np.arange(len(categories))
    width = 0.35

    bars1 = ax.bar(x - width/2, counts_raw, width, label='原始数据', 
                   color='lightgray', edgecolor='black')
    bars2 = ax.bar(x + width/2, counts_filt, width, label=f'过滤后', 
                   color='skyblue', edgecolor='black')

    ax.set_ylabel('位点数量', fontsize=12)
    ax.set_title('A: 甲基化位点数量统计', fontsize=14, fontweight='bold')
    ax.set_xticks(x)
    ax.set_xticklabels(['5mC@CG', '5mC@non-CG', '4mC@CG', '4mC@non-CG'])
    ax.legend()
    ax.grid(axis='y', alpha=0.3)

    # 添加数值标签
    for bars in [bars1, bars2]:
        for bar in bars:
            height = bar.get_height()
            if height > 0:
                ax.annotate(f'{int(height)}',
                           xy=(bar.get_x() + bar.get_width() / 2, height),
                           xytext=(0, 3),
                           textcoords="offset points",
                           ha='center', va='bottom', fontsize=9)

    plt.tight_layout()
    plt.savefig(output_file, dpi=300, bbox_inches='tight')
    plt.close()
    print(f"  ✓ 保存: {output_file}")

def plot_modification_distribution(df, output_file):
    """图2: 修饰率分布"""
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    axes = axes.flatten()

    categories = ['5mC_CG', '5mC_nonCG', '4mC_CG', '4mC_nonCG']
    colors = {'5mC_CG': '#1f77b4', '5mC_nonCG': '#ff7f0e', 
              '4mC_CG': '#2ca02c', '4mC_nonCG': '#d62728'}
    titles = {'5mC_CG': '5mC @ CG', '5mC_nonCG': '5mC @ non-CG',
              '4mC_CG': '4mC @ CG', '4mC_nonCG': '4mC @ non-CG'}

    for idx, cat in enumerate(categories):
        data = df[df['category']==cat]['percent_modified']
        if len(data) > 0:
            sns.histplot(data=data, ax=axes[idx], color=colors[cat], 
                        bins=50, kde=True, edgecolor='black')
            axes[idx].axvline(x=50, color='red', linestyle='--', alpha=0.7, label='阈值 50%')
            axes[idx].axvline(x=70, color='darkred', linestyle='--', alpha=0.7, label='阈值 70%')
            axes[idx].set_xlabel('修饰率 (%)')
            axes[idx].set_ylabel('位点数量')
            axes[idx].set_title(f'{titles[cat]} (n={len(data)})')
            axes[idx].legend(fontsize=8)
            axes[idx].grid(axis='y', alpha=0.3)
        else:
            axes[idx].text(0.5, 0.5, '无数据', ha='center', va='center', 
                          transform=axes[idx].transAxes, fontsize=14)

    plt.suptitle('B: 修饰率分布', fontsize=14, fontweight='bold', y=0.995)
    plt.tight_layout()
    plt.savefig(output_file, dpi=300, bbox_inches='tight')
    plt.close()
    print(f"  ✓ 保存: {output_file}")

def plot_coverage_distribution(df, output_file):
    """图3: 覆盖度分布"""
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

    # 5mC
    data_5mC_CG = df[df['category']=='5mC_CG']['Nvalid_cov']
    data_5mC_nonCG = df[df['category']=='5mC_nonCG']['Nvalid_cov']

    if len(data_5mC_CG) > 0:
        sns.kdeplot(data=data_5mC_CG, ax=axes[0], label='5mC@CG', 
                   color='#1f77b4', linewidth=2)
    if len(data_5mC_nonCG) > 0:
        sns.kdeplot(data=data_5mC_nonCG, ax=axes[0], label='5mC@non-CG', 
                   color='#ff7f0e', linewidth=2)

    axes[0].set_xlabel('覆盖度 (Reads)')
    axes[0].set_ylabel('密度')
    axes[0].set_title('C: 5mC 覆盖度分布')
    axes[0].legend()
    axes[0].grid(axis='y', alpha=0.3)

    # 4mC
    data_4mC_CG = df[df['category']=='4mC_CG']['Nvalid_cov']
    data_4mC_nonCG = df[df['category']=='4mC_nonCG']['Nvalid_cov']

    if len(data_4mC_CG) > 0:
        sns.kdeplot(data=data_4mC_CG, ax=axes[1], label='4mC@CG', 
                   color='#2ca02c', linewidth=2)
    if len(data_4mC_nonCG) > 0:
        sns.kdeplot(data=data_4mC_nonCG, ax=axes[1], label='4mC@non-CG', 
                   color='#d62728', linewidth=2)

    axes[1].set_xlabel('覆盖度 (Reads)')
    axes[1].set_ylabel('密度')
    axes[1].set_title('D: 4mC 覆盖度分布')
    axes[1].legend()
    axes[1].grid(axis='y', alpha=0.3)

    plt.tight_layout()
    plt.savefig(output_file, dpi=300, bbox_inches='tight')
    plt.close()
    print(f"  ✓ 保存: {output_file}")

def generate_statistics_table(df_raw, df_filt, output_file):
    """图4: 统计汇总表"""
    fig, ax = plt.subplots(figsize=(10, 4))
    ax.axis('off')

    categories = ['5mC_CG', '5mC_nonCG', '4mC_CG', '4mC_nonCG']
    stats_data = []

    for cat in categories:
        cat_raw = df_raw[df_raw['category']==cat]
        cat_filt = df_filt[df_filt['category']==cat]

        if len(cat_raw) > 0:
            retention = len(cat_filt) / len(cat_raw) * 100
        else:
            retention = 0

        stats_data.append({
            '类型': cat.replace('_', ' @ '),
            '原始位点数': len(cat_raw),
            '过滤后位点数': len(cat_filt),
            '保留率(%)': f"{retention:.1f}",
            '中位修饰率(%)': f"{cat_filt['percent_modified'].median():.2f}" if len(cat_filt)>0 else "N/A",
            '平均修饰率(%)': f"{cat_filt['percent_modified'].mean():.2f}" if len(cat_filt)>0 else "N/A",
            '中位覆盖度': f"{cat_filt['Nvalid_cov'].median():.0f}" if len(cat_filt)>0 else "N/A",
        })

    df_stats = pd.DataFrame(stats_data)
    table = ax.table(cellText=df_stats.values, colLabels=df_stats.columns, 
                    cellLoc='center', loc='center', bbox=[0, 0, 1, 1])
    table.auto_set_font_size(False)
    table.set_fontsize(9)
    table.scale(1.2, 1.8)

    # 设置表头样式
    for i in range(len(df_stats.columns)):
        table[(0, i)].set_facecolor('#4472C4')
        table[(0, i)].set_text_props(color='white', fontweight='bold')

    plt.savefig(output_file, dpi=300, bbox_inches='tight')
    plt.close()
    print(f"  ✓ 保存: {output_file}")

def main():
    input_bed = sys.argv[1]
    output_dir = sys.argv[2]

    print("📊 生成可视化图表...")

    # 加载数据
    df_raw = load_data(input_bed)
    df_raw = classify_modifications(df_raw)

    # 加载过滤后数据
    filtered_files = {
        '5mC_CG': f"{output_dir}/filtered/5mC_CG_filtered.bed",
        '5mC_nonCG': f"{output_dir}/filtered/5mC_nonCG_filtered.bed",
        '4mC_CG': f"{output_dir}/filtered/4mC_CG_filtered.bed",
        '4mC_nonCG': f"{output_dir}/filtered/4mC_nonCG_filtered.bed"
    }

    # 合并过滤后的数据
    df_list = []
    for cat, fpath in filtered_files.items():
        if os.path.exists(fpath) and os.path.getsize(fpath) > 0:
            df_tmp = load_data(fpath)
            df_tmp['mod_type'] = '5mC' if '5mC' in cat else '4mC'
            df_tmp['motif_type'] = 'CG' if 'CG' in cat else 'non-CG'
            df_tmp['category'] = cat
            df_list.append(df_tmp)

    if df_list:
        df_filt = pd.concat(df_list, ignore_index=True)
    else:
        df_filt = df_raw.iloc[0:0]  # 空DataFrame

    # 生成图表
    os.makedirs(f"{output_dir}/figures", exist_ok=True)

    plot_counts_comparison(df_raw, df_filt, f"{output_dir}/figures/Fig1_counts.png")
    plot_modification_distribution(df_filt, f"{output_dir}/figures/Fig2_mod_distribution.png")
    plot_coverage_distribution(df_filt, f"{output_dir}/figures/Fig3_coverage.png")
    generate_statistics_table(df_raw, df_filt, f"{output_dir}/figures/Fig4_statistics.png")

    print("✅ 所有图表生成完成！")

if __name__ == "__main__":
    if len(sys.argv) != 3:
        print(f"用法: {sys.argv[0]} <input.bed> <output_dir>")
        sys.exit(1)
    main()
PYTHON_SCRIPT

    chmod +x "$OUTPUT_DIR/generate_plots.py"

    # 运行Python脚本
    if command -v python3 &> /dev/null; then
        log_info "运行Python可视化脚本..."
        python3 "$OUTPUT_DIR/generate_plots.py" "$INPUT_BED" "$OUTPUT_DIR" 2>&1 | tee "$OUTPUT_DIR/logs/visualization.log"
    else
        log_warn "Python3未安装，跳过可视化生成"
    fi
else
    log_info "可视化已禁用"
fi

#-------------------------------------------------------------------------------
#  生成最终报告
#-------------------------------------------------------------------------------

log_step "生成分析总结报告"

cat > "$OUTPUT_DIR/reports/FINAL_REPORT.txt" << EOF
================================================================================
                    DNA甲基化分析总结报告
================================================================================
生成时间: $(date)
分析脚本: $0

📁 输入文件:
   - BED文件: $INPUT_BED
   - 参考基因组: $REF_GENOME

⚙️  分析参数:
   - 最小覆盖度: $MIN_COVERAGE
   - 5mC修饰率阈值: ${MIN_MOD_5MC}%
   - 4mC修饰率阈值: ${MIN_MOD_4MC}%
   - 侧翼序列: ±${FLANK_SIZE}bp
   - 线程数: $THREADS

📊 数据统计:
$(cat "$OUTPUT_DIR/reports/filtering_summary.txt")

📂 输出文件结构:
   $OUTPUT_DIR/
   ├── filtered/              # 过滤后的bed文件
   │   ├── 5mC_CG_filtered.bed
   │   ├── 5mC_nonCG_filtered.bed
   │   ├── 4mC_CG_filtered.bed
   │   └── 4mC_nonCG_filtered.bed
   ├── sequences/             # 提取的fasta序列
   ├── motif_analysis/        # Motif富集分析结果
   │   └── homer/
   ├── figures/               # 可视化图表
   └── reports/               # 统计报告

🔬 主要发现:
   - 5mC主要富集在: $(if [[ -f "$OUTPUT_DIR/motif_analysis/homer/5mC_nonCG/knownResults.txt" ]]; then head -2 "$OUTPUT_DIR/motif_analysis/homer/5mC_nonCG/knownResults.txt" | tail -1 | cut -f1; else "待分析"; fi)
   - 4mC主要富集在: $(if [[ -f "$OUTPUT_DIR/motif_analysis/homer/4mC_nonCG/knownResults.txt" ]]; then head -2 "$OUTPUT_DIR/motif_analysis/homer/4mC_nonCG/knownResults.txt" | tail -1 | cut -f1; else "待分析"; fi)

================================================================================
                              分析完成！
================================================================================
EOF

log_info "最终报告: $OUTPUT_DIR/reports/FINAL_REPORT.txt"
cat "$OUTPUT_DIR/reports/FINAL_REPORT.txt"

#-------------------------------------------------------------------------------
#  清理临时文件
#-------------------------------------------------------------------------------

log_step "清理临时文件"

# 可选：删除大型中间文件
# rm -rf "$OUTPUT_DIR/sequences/"*_regions.bed

log_info "分析流程完成！✅"
log_info "查看最终报告: cat $OUTPUT_DIR/reports/FINAL_REPORT.txt"

exit 0

---

🚀 使用方法

1️⃣ 保存脚本

# 保存为 dna_methylation_analysis.sh
nano dna_methylation_analysis.sh
# 粘贴上面的代码，保存退出

# 赋予执行权限
chmod +x dna_methylation_analysis.sh

2️⃣ 运行分析

# 基本用法
bash dna_methylation_analysis.sh \
    -i An6_4mC_5mC_raw.bed \
    -r bacass_out/Medaka/An6-unicycler-medaka_polished_genome.fa \
    -o An6_methylation_results

# 自定义参数
bash dna_methylation_analysis.sh \
    -i An6_4mC_5mC_raw.bed \
    -r genome.fa \
    -o results \
    -c 10 \          # 最小覆盖度
    -m 70 \          # 5mC阈值
    -n 50 \          # 4mC阈值
    -f 50 \          # 侧翼序列
    -t 16            # 线程数

# 禁用某些功能
bash dna_methylation_analysis.sh \
    -i input.bed -r genome.fa \
    --no-homer       # 不运行HOMER
    --no-viz         # 不生成图表

3️⃣ 查看结果

# 查看最终报告
cat An6_methylation_results/reports/FINAL_REPORT.txt

# 查看图表
ls -lh An6_methylation_results/figures/

# 查看motif结果
ls -lh An6_methylation_results/motif_analysis/homer/

---

📦 依赖安装

# Conda环境（推荐）
conda create -n methylation python=3.9
conda activate methylation

# 必需工具
conda install -c bioconda bedtools samtools

# 可选：HOMER
conda install -c bioconda homer

# 可选：MEME Suite
conda install -c bioconda meme

# Python库
pip install pandas numpy matplotlib seaborn

---

📋 输出文件说明

An6_methylation_results/
├── filtered/                    # 过滤后的bed文件
│   ├── 5mC_CG_filtered.bed     # 5mC @ CG位点
│   ├── 5mC_nonCG_filtered.bed  # 5mC @ non-CG位点
│   ├── 4mC_CG_filtered.bed     # 4mC @ CG位点
│   └── 4mC_nonCG_filtered.bed  # 4mC @ non-CG位点
├── sequences/                   # 提取的序列
│   ├── *_sequences.fa          # 目标序列
│   └── background_sequences.fa # 背景序列
├── motif_analysis/
│   └── homer/
│       ├── 5mC_CG/             # HOMER结果
│       ├── 5mC_nonCG/
│       ├── 4mC_CG/
│       └── 4mC_nonCG/
├── figures/                     # 可视化图表
│   ├── Fig1_counts.png         # 位点数量对比
│   ├── Fig2_mod_distribution.png # 修饰率分布
│   ├── Fig3_coverage.png       # 覆盖度分布
│   └── Fig4_statistics.png     # 统计表
└── reports/                     # 统计报告
    ├── analysis_config.txt     # 分析配置
    ├── filtering_summary.txt   # 过滤统计
    └── FINAL_REPORT.txt        # 最终报告

需要我帮您调整任何参数或添加其他功能吗？🔧✨

🔍 nf-core/methylong 内置 DMR 选项：能替代哪些独立工具？

✅ 简短回答

您之前需要的功能	methylong 内置选项	能否替代独立工具？
位点水平差异分析 (对标 TRANSIT ANOVA)	--population_dmrer dss + --dmr_a/--dmr_b	✅ 完全可以替代独立 DSS
区域水平聚合比较	--population_dmrer dss (自动区域聚合)	✅ 部分替代 methylKit
连续低/高甲基化区域检测 (对标 Tn5Gaps)	❌ 无专门功能	❌ 仍需 metilene 或自定义脚本
基因/操纵子注释关联	❌ 不直接整合 GFF	❌ 仍需 bedtools/ChIPseeker
出版级可视化	❌ 只输出 bed/bigWig	❌ 仍需 pyGenomeTracks/IGV

🎯 核心结论：methylong 的 --population_dmrer dss 可以完全替代您独立运行 DSS 进行组间差异甲基化分析，无需额外工具！

---

📋 methylong DMR 选项详解 (对标 TRANSIT)

🔹 参数功能映射表

methylong 参数	作用	TRANSIT 对应功能	统计方法
--skip_snvs	跳过 SNP calling + phasing	– (细菌无需)	–
--haplotype_dmrer [dss/modkit]	单倍型水平 DMR 工具	– (细菌无单倍型)	贝叶斯层次模型 (DSS)
--population_dmrer [dss/modkit]	⭐ 群体水平 DMR 工具	✅ ANOVA 多组比较	贝塔 – 二项检验 + FDR
--dmr_population_scale	启用多样本组间比较	✅ 条件间差异分析	同上
--dmr_a, --dmr_b	定义两个比较组的样本名	✅ 指定对比条件	–

---

🚀 用 methylong 内置 DSS 替代独立运行的完整命令

场景：细菌 6mA，An6 vs BG5 差异分析

nextflow run nf-core/methylong \
  -r 2.0.0 \
  -profile docker \
  --input samplesheet_6mA.csv \          # 包含 An6 和 BG5 两行
  --outdir methylome_out_6mA \
  --ont_aligner minimap2 \
  --m6a \
  --skip_snvs \                          # 细菌无需 phasing
  --dmr_population_scale \               # ⭐ 启用群体水平比较
  --population_dmrer dss \               # ⭐ 使用 DSS 进行统计 (默认值，可省略)
  --dmr_a "An6" \                        # ⭐ 定义组 A (samplesheet 中的 group 列)
  --dmr_b "BG5" \                        # ⭐ 定义组 B
  -resume \
  -work-dir methylome_out_6mA/work

🔹 samplesheet_6mA.csv 格式要求

group,sample,path,ref,method
An6,An6,/path/An6_6mA_mapped.mod.bam,/path/genome.fa,ont
BG5,BG5,/path/BG5_6mA_mapped.mod.bam,/path/genome.fa,ont

⚠️ 关键：group 列的值 (An6, BG5) 必须与 --dmr_a/--dmr_b 完全一致！

---

📁 输出文件：methylong 内置 DMR 结果位置

methylome_out_6mA/
└── results/
    └── dmr/
        ├── population_scale/
        │   ├── An6_vs_BG5.DSS_DMLs.tsv      # ⭐ 位点水平差异结果 (对标 TRANSIT anova_out.xls)
        │   ├── An6_vs_BG5.DSS_DMRs.tsv      # 区域水平差异结果
        │   └── DSS_summary.html             # 简要统计摘要
        └── qc/
            └── dmr_multiqc_report.html      # DMR 分析质控

🔹 DSS_DMLs.tsv 列解读 (对标 TRANSIT ANOVA 输出)

列名	含义	TRANSIT 对应列
chr, pos	基因组坐标	Orf, TA_pos
coverage_A, coverage_B	两组覆盖深度	TAs
mod_count_A, mod_count_B	两组修饰计数	Mean_[condition]
methylation_A, methylation_B	两组甲基化概率 (0-100)	–
diff	甲基化差异 (A – B)	LFC_[condition]
pval, p.adj	原始/FDR 校正 p 值	Pval, Padj
significant	是否显著 (FDR20%)	status

---

⚠️ methylong 内置 DMR 的局限性 (仍需外部工具的场景)

❌ 无法替代的功能 1：连续区域检测 (对标 Tn5Gaps)

# methylong 输出的是"离散位点"差异，不是"连续非修饰区域"
# 如果您需要检测类似 Tn5Gaps 的"长片段无 6mA 区域"，仍需：

metilene \
  -i An6.bedgraph BG5.bedgraph \
  -g genome.sizes \
  -o continuous_dmr.tsv \
  -p 0.05 -d 0.2

❌ 无法替代的功能 2：基因/操纵子水平注释

# methylong 不直接整合 GFF 注释
# 如需"哪些基因的 6mA 水平显著改变"，仍需：

bedtools intersect -a An6_vs_BG5.DSS_DMLs.tsv -b genome.gff -wa -wb > genes_with_diff_6mA.tsv

❌ 无法替代的功能 3：高级可视化

# methylong 只输出 bed/bigWig，不生成出版级图
# 如需 Circos 式基因组视图，仍需：

pygenometracks --tracks tracks.ini --output final_figure.png

---

🎯 推荐工作流：methylong 内置 + 最小外部工具

graph LR
    A[modkit pileup] --> B[methylong with --population_dmrer dss]
    B --> C[An6_vs_BG5.DSS_DMLs.tsv]

    C --> D{需要连续区域？}
    D -->|是 | E[metilene for continuous DMRs]
    D -->|否 | F[直接进入注释]

    E & F --> G[bedtools + genome.gff for gene annotation]
    G --> H[pyGenomeTracks for visualization]

✅ 最小命令集：

# 1️⃣ 主分析：methylong 内置 DSS (替代独立 DSS 运行)
nextflow run nf-core/methylong ... --dmr_population_scale --dmr_a An6 --dmr_b BG5 ...

# 2️⃣ (可选) 连续区域检测：仍需 metilene
awk '$4 == "a" {printf "%s\t%d\t%d\t%.2f\n", $1, $2, $3, $11}' An6_6mA_raw.bed > An6.bg
metilene -i An6.bg BG5.bg -g genome.sizes -o continuous_dmr.tsv

# 3️⃣ (可选) 基因注释：仍需 bedtools
bedtools intersect -a methylome_out_6mA/results/dmr/population_scale/An6_vs_BG5.DSS_DMLs.tsv \
  -b genome.gff -wa -wb > annotated_diff_6mA.tsv

---

💡 一句话总结

🔑 methylong 内置 DMR 选项可以：

• ✅ 完全替代独立 DSS 进行位点水平差异分析 (--population_dmrer dss)
• ✅ 自动完成：标准化 + 贝塔 – 二项检验 + FDR 校正 + 结果输出

🔑 但仍需外部工具：

• ❌ 连续区域检测 (→ metilene)
• ❌ 基因注释关联 (→ bedtools/ChIPseeker)
• ❌ 出版级可视化 (→ pyGenomeTracks)

🎯 您的最优路径：先用 methylong 内置 DSS 拿到差异位点列表，再按需添加 1-2 个外部工具完成注释和可视化！

如需我帮您编写 samplesheet 或调试 methylong DMR 参数，请随时告知！🚀

---

---

🔍 结果分析：非 CG 位点更多是完全合理的！

您的结果 完全符合逻辑，甚至反映了细菌基因组的典型特征。让我详细解释：

---

📊 您的结果汇总

修饰类型	位点数量	占比	生物学解释
5mC @ CG	437	33%	经典 CpG 甲基化
5mC @ non-CG	870	67%	⭐ 非 CG 位点更多，正常！
4mC @ CG	7	17%	4mC 很少在 CG 上
4mC @ non-CG	34	83%	⭐ 4mC 偏好非 CG 位点

---

✅ 为什么 non-CG 位点更多？三大原因

🔹 原因 1：序列组成决定数量基数

假设基因组碱基随机分布（各 25%）：

所有可能的二核苷酸：16 种 (AA, AC, AG, AT, CA, CC, CG, CT, ...)

其中含 C 的位点：
├─ CG: 1 种 → 占所有 C 位点的 ~25%
├─ CA, CC, CT: 3 种 → 占所有 C 位点的 ~75% ← non-CG 位点天生更多！

因此：即使甲基化概率相同，non-CG 位点数量也会是 CG 的 ~3 倍

🔹 原因 2：细菌甲基化偏好特定非 CG 基序

细菌的限制修饰系统通常靶向 特定回文序列，而非随机 CG：

常见细菌甲基化基序	修饰类型	示例细菌
GATC	6mA / 4mC	E. coli Dam 甲基化酶
CCWGG	5mC	H. influenzae Dcm 甲基化酶
GANTC	6mA	Caulobacter CcrM 甲基化酶
RCGY	4mC	多种细菌

→ 如果您的细菌有类似系统，非 CG 位点会被优先甲基化！

🔹 原因 3：--motif C 0 捕获了所有胞嘧啶

您使用的命令：

modkit pileup ... --motif CG 0 --motif C 0 ...

• --motif C 0 = 检测 每一个胞嘧啶（无论后面是什么碱基）
• 因此 m,C,0 包含：CG + CA + CC + CT 所有 C 位点
• 排除 CG 后，剩余的 non-CG 位点自然更多 ✅

---

🔬 验证：这些 non-CG 信号是真实的吗？

运行以下命令检查数据质量：

# 1️⃣ 比较 CG vs non-CG 的修饰率分布
echo "=== 5mC 修饰率中位数 ==="
echo -n "CG:    "
awk -F'\t' '$4 ~ /^m,CG,/ {print $11}' An6_4mC_5mC_raw.bed | sort -n | awk '{a[NR]=$1} END{print a[int(NR/2)]}'
echo -n "non-CG:"
awk -F'\t' '$4 ~ /^m,C,/ && $4 !~ /^m,CG,/ {print $11}' An6_4mC_5mC_raw.bed | sort -n | awk '{a[NR]=$1} END{print a[int(NR/2)]}'

# 2️⃣ 检查 non-CG 位点的序列背景（提取上游10bp）
# 需要先有基因组 fasta 索引
samtools faidx bacass_out/Medaka/An6-unicycler-medaka_polished_genome.fa

# 提取 non-CG 5mC 位点的序列上下文
awk -F'\t' '$4 ~ /^m,C,/ && $4 !~ /^m,CG,/ && $11 >= 70 && $5 >= 10 {
  print $1":"$2+1"-"$2+11
}' An6_4mC_5mC_raw.bed | head -20 > nonCG_sites.txt

# 3️⃣ 用 MEME/HOMER 做 motif 富集（可选）
# 如果 870 个位点富集到特定序列，说明是真实生物学信号！

---

🎯 分析建议：如何处理这 4 类数据？

✅ 推荐策略：分层分析

# 📁 第一层：核心分析（高置信度 CG 位点）
# 用于主图、主要结论
awk -F'\t' '$4 ~ /^m,CG,/ && $11 >= 70 && $5 >= 10' An6_4mC_5mC_raw.bed > An6_5mC_CG_high.bed
awk -F'\t' '$4 ~ /^21839,CG,/ && $11 >= 50 && $5 >= 10' An6_4mC_5mC_raw.bed > An6_4mC_CG_high.bed

# 📁 第二层：探索性分析（non-CG 位点）
# 用于补充材料、motif 发现
awk -F'\t' '$4 ~ /^m,C,/ && $4 !~ /^m,CG,/ && $11 >= 70 && $5 >= 10' An6_4mC_5mC_raw.bed > An6_5mC_nonCG.bed
awk -F'\t' '$4 ~ /^21839,C,/ && $4 !~ /^21839,CG,/ && $11 >= 50 && $5 >= 10' An6_4mC_5mC_raw.bed > An6_4mC_nonCG.bed

# 📁 第三层：全量数据（用于基因组浏览器）
# 不过滤，供 IGV 可视化
cp An6_4mC_5mC_raw.bed An6_all_modified.bed

📊 可视化建议

数据类型	推荐图表	工具
5mC @ CG	基因组分布 + 修饰率直方图	deepTools + ggplot2
4mC @ non-CG	motif 富集 + 热图	HOMER + seaborn
两者对比	Venn 图 + 散点图	R/ComplexHeatmap

---

⚠️ 注意事项：避免重复计数！

由于 m,C,0 包含 m,CG,0，不要简单相加：

# ❌ 错误：直接相加会重复计算 CG 位点
total_5mC = CG_count + nonCG_count  # 重复！

# ✅ 正确：只用 C  motif 统计全部 5mC
awk -F'\t' '$4 ~ /^m,C,/ && $11 >= 70 && $5 >= 10' An6_4mC_5mC_raw.bed | wc -l
# 或：CG + nonCG（确保互斥）
awk -F'\t' '($4 ~ /^m,CG,/ || ($4 ~ /^m,C,/ && $4 !~ /^m,CG,/)) && $11 >= 70 && $5 >= 10' An6_4mC_5mC_raw.bed | wc -l

---

🧪 快速诊断脚本（一键运行）

#!/bin/bash
FILE="An6_4mC_5mC_raw.bed"

echo "🔍 数据质量诊断"
echo "================"
echo -e "\n📈 各类型原始行数（无过滤）:"
echo "5mC @ CG:    $(grep -c 'm,CG,' $FILE)"
echo "5mC @ C:     $(grep -c 'm,C,' $FILE)"
echo "4mC @ CG:    $(grep -c '21839,CG,' $FILE)"
echo "4mC @ C:     $(grep -c '21839,C,' $FILE)"

echo -e "\n📊 修饰率统计 (列11):"
for label in "5mC_CG" "5mC_nonCG" "4mC_CG" "4mC_nonCG"; do
  case $label in
    "5mC_CG") pattern="^m,CG," ;;
    "5mC_nonCG") pattern="^m,C,"; exclude="m,CG," ;;
    "4mC_CG") pattern="^21839,CG," ;;
    "4mC_nonCG") pattern="^21839,C,"; exclude="21839,CG," ;;
  esac

  if [ -n "$exclude" ]; then
    stats=$(awk -F'\t' -v p="$pattern" -v e="$exclude" '$4 ~ p && $4 !~ e {print $11}' $FILE | \
      awk '{s+=$1; n++} END{if(n>0) printf "n=%d, median=%.2f, mean=%.2f", n, (n%2?$(int(n/2)+1):($(n/2)+$(n/2+1))/2), s/n}')
  else
    stats=$(awk -F'\t' -v p="$pattern" '$4 ~ p {print $11}' $FILE | \
      awk '{s+=$1; n++} END{if(n>0) printf "n=%d, median=%.2f, mean=%.2f", n, (n%2?$(int(n/2)+1):($(n/2)+$(n/2+1))/2), s/n}')
  fi
  echo "$label: $stats"
done

---

🎯 结论

✅ 870 > 437 是完全合理的！
非 CG 胞嘧啶在基因组中本就更多，加上细菌甲基化系统常靶向特定非 CG 基序，导致 non-CG 修饰位点数量超过 CG 位点。

下一步建议：

1. 用 nonCG_sites.txt 做 motif 富集，看看 870 个位点是否聚集在特定序列
2. 如果富集到已知细菌甲基化基序（如 GATC），说明信号真实可靠 🔬
3. 在论文/报告中明确区分 CG 和 non-CG 分析，避免混淆

需要我帮您写 motif 富集分析脚本，或生成对比图表吗？📈✨