如何批量对多分辨率.mcool文件迭代运行EagleC的predictSV工具

需求说明 #

• .mcool文件包含多分辨率矩阵数据
• 对./input/*.mcool中的每个文件，当cooler ls "$mcool_file"输出结果的最后一个/后为5000、10000或50000时，运行EagleC的predictSV命令
• 单个文件的正确命令示例：

predictSV --hic-5k SKNAS-MboI-allReps-filtered.mcool::/resolutions/5000 
          --hic-10k SKNAS-MboI-allReps-filtered.mcool::/resolutions/10000 
          --hic-50k SKNAS-MboI-allReps-filtered.mcool::/resolutions/50000 
          -O SK-N-AS -g hg38 --balance-type CNV --output-format full 
          --prob-cutoff-5k 0.8 --prob-cutoff-10k 0.8 --prob-cutoff-50k 0.99999

错误脚本分析 #

尝试脚本1 #

for mcool_file in ./input/*.mcool; do
  while IFS= read -r id; do
    id_suffix=${id##*/}
    case $id_suffix in 5000|10000|50000)
      num=${id_suffix:0:1}
      predictSV \
        --hic-"${num}k" "$id_suffix" \
        -g hg38 \
        -O basename "${id%%.*}" \
        --balance-type CNV \
        --output-format full \
        --prob-cutoff-5k 0.8 \
        --prob-cutoff-10k 0.8 \
        --prob-cutoff-50k 0.99999
      ;;
    esac
  done < <(cooler ls "$mcool_file")
done

报错信息

usage: predictSV [-h] [-v] [--hic-5k HIC_5K] [--hic-10k HIC_10K]
                 [--hic-50k HIC_50K] [-O OUTPUT_PREFIX]
                 [-g {hg38,hg19,chm13,other}] [-C [CHROMS ...]]
                 [--balance-type {ICE,CNV,Raw}]
                 [--output-format {full,NeoLoopFinder}]
                 [--prob-cutoff-5k PROB_CUTOFF_5K]
                 [--prob-cutoff-10k PROB_CUTOFF_10K]
                 [--prob-cutoff-50k PROB_CUTOFF_50K]
predictSV: error: unrecognized arguments:
usage: predictSV [-h] [-v] [--hic-5k HIC_5K] [--hic-10k HIC_10K]
                 [--hic-50k HIC_50K] [-O OUTPUT_PREFIX]
                 [-g {hg38,hg19,chm13,other}] [-C [CHROMS ...]]
                 [--balance-type {ICE,CNV,Raw}]
                 [--output-format {full,NeoLoopFinder}]
                 [--prob-cutoff-5k PROB_CUTOFF_5K]
                 [--prob-cutoff-10k PROB_CUTOFF_10K]
                 [--prob-cutoff-50k PROB_CUTOFF_50K]
predictSV: error: unrecognized arguments: --hic-1k 10000
usage: predictSV [-h] [-v] [--hic-5k HIC_5K] [--hic-10k HIC_10K]
                 [--hic-50k HIC_50K] [-O OUTPUT_PREFIX]
                 [-g {hg38,hg19,chm13,other}] [-C [CHROMS ...]]
                 [--balance-type {ICE,CNV,Raw}]
                 [--output-format {full,NeoLoopFinder}]
                 [--prob-cutoff-5k PROB_CUTOFF_5K]
                 [--prob-cutoff-10k PROB_CUTOFF_10K]
                 [--prob-cutoff-50k PROB_CUTOFF_50K]
predictSV: error: unrecognized arguments:
usage: predictSV [-h] [-v] [--hic-5k HIC_5K] [--hic-10k HIC_10K]
                 [--hic-50k HIC_50K] [-O OUTPUT_PREFIX]
                 [-g {hg38,hg19,chm13,other}] [-C [CHROMS ...]]
                 [--balance-type {ICE,CNV,Raw}]
                 [--output-format {full,NeoLoopFinder}]
                 [--prob-cutoff-5k PROB_CUTOFF_5K]
                 [--prob-cutoff-10k PROB_CUTOFF_10K]
                 [--prob-cutoff-50k PROB_CUTOFF_50K]
predictSV: error: unrecognized arguments:
usage: predictSV [-h] [-v] [--hic-5k HIC_5K] [--hic-10k HIC_10K]
                 [--hic-50k HIC_50K] [-O OUTPUT_PREFIX]
                 [-g {hg38,hg19,chm13,other}] [-C [CHROMS ...]]
                 [--balance-type {ICE,CNV,Raw}]
                 [--output-format {full,NeoLoopFinder}]
                 [--prob-cutoff-5k PROB_CUTOFF_5K]
                 [--prob-cutoff-10k PROB_CUTOFF_10K]
                 [--prob-cutoff-50k PROB_CUTOFF_50K]
..

问题点

1. num=${id_suffix:0:1}对10000取第一个字符得到1，生成无效参数--hic-1k，predictSV仅支持--hic-5k/--hic-10k/--hic-50k
2. --hic-"${num}k"后的参数用了$id_suffix，但实际需要完整的分辨率路径（即mcool_file::/resolutions/xxx），而非仅数字
3. -O basename "${id%%.*}"语法错误，basename需用$()包裹执行，否则会被当成输出前缀的一部分

尝试脚本2 #

for mcool_file in ./input/*.mcool; do
  while IFS= read -r id; do
    id_suffix=${id##*/}
    case $id_suffix in 5000|10000|50000)
      num=${id_suffix:0:1}
      predictSV \
        if [ $id_suffix=5000 ];
          then --hic-5k "$id_suffix";
        elif [ $id_suffix=10000 ];
          then --hic-10k "$id_suffix";
        else:
          --hic-50k "$id_suffix";
        fi \
        -g hg38 \
        -O basename "${id%%.*}" \
        --balance-type CNV \
        --output-format full \
        --prob-cutoff-5k 0.8 \
        --prob-cutoff-10k 0.8 \
        --prob-cutoff-50k 0.99999
      ;;
    esac
  done < <(cooler ls "$mcool_file")
done

问题点

1. 在predictSV命令行中直接写入if语句属于语法错误，shell无法识别，需用参数拼接方式动态生成选项
2. 同样存在$id_suffix代替完整路径、basename未正确执行的问题
3. 每个分辨率单独执行一次predictSV，但实际要求一次命令同时传入三个分辨率参数

正确脚本实现 #

思路：对每个mcool文件，先收集所有符合条件的分辨率路径，再一次性拼接成predictSV的命令参数后执行。

for mcool_file in ./input/*.mcool; do
    # 获取输出前缀：提取mcool文件名（去掉路径和后缀）
    output_prefix=$(basename "${mcool_file%.mcool}")
    # 初始化三个分辨率的变量
    hic_5k=""
    hic_10k=""
    hic_50k=""

    # 遍历当前mcool文件的所有分辨率
    while IFS= read -r id; do
        id_suffix=${id##*/}
        case $id_suffix in
            5000)
                hic_5k="--hic-5k $id"
                ;;
            10000)
                hic_10k="--hic-10k $id"
                ;;
            50000)
                hic_50k="--hic-50k $id"
                ;;
        esac
    done < <(cooler ls "$mcool_file")

    # 检查三个分辨率是否都存在，存在则执行命令
    if [ -n "$hic_5k" ] && [ -n "$hic_10k" ] && [ -n "$hic_50k" ]; then
        predictSV \
            $hic_5k \
            $hic_10k \
            $hic_50k \
            -O "$output_prefix" \
            -g hg38 \
            --balance-type CNV \
            --output-format full \
            --prob-cutoff-5k 0.8 \
            --prob-cutoff-10k 0.8 \
            --prob-cutoff-50k 0.99999
    else
        echo "警告：文件 $mcool_file 缺少5000/10000/50000中的部分分辨率，跳过处理"
    fi
done

脚本说明 #

1. 对每个mcool文件，先提取文件名作为输出前缀，避免basename的语法错误
2. 遍历cooler ls的结果，分别匹配三个分辨率，保存对应的完整参数
3. 确保三个分辨率都存在时才执行predictSV，避免因参数不全报错
4. 使用完整的分辨率路径作为--hic-xxk的参数，符合命令要求

内容的提问来源于stack exchange，提问作者Anon