利用samtools基于染色体位置信息批量生成pileup文件的脚本需求

嘿，我完全懂你现在的困扰——刚接触Linux又没编程基础，要把这个批量处理的逻辑写成脚本确实有点懵。别慌，我给你写一个清晰的脚本，每一步都加注释，保证你能看懂也能用上。

先理清楚需求逻辑 #

我们要做的是：

• 遍历两种组织（colon/muscle）和22条染色体（1-22），找到对应的44个文件
• 对每个文件里的每一行，提取位置x，计算区间[x-5, x+5]
• 针对这个区间运行samtools mpileup命令

完整脚本示例 #

先把脚本写出来，后面再逐行解释：

#!/bin/bash

# 1. 定义需要用到的变量
# 替换成你的参考基因组文件的绝对路径，比如/hg38/hg38.fa
REFERENCE="/path/to/your/reference/genome.fa"
# 组织列表
TISSUES=("colon" "muscle")
# 染色体范围1-22
CHROMOSOMES=$(seq 1 22)

# 2. 外层循环：遍历组织和染色体
for TISSUE in "${TISSUES[@]}"; do
  for c in $CHROMOSOMES; do
    # 拼接当前要处理的文件名
    INPUT_FILE="chr${c}.${TISSUE}_run1_en2hic_PE1.bam.sorted.bam.breaks_COL1_and_COL2_ONLY"
    
    # 先检查文件是否存在，避免报错
    if [ ! -f "$INPUT_FILE" ]; then
      echo "警告：文件 $INPUT_FILE 不存在，跳过..."
      continue
    fi
    
    echo "正在处理文件：$INPUT_FILE"
    
    # 3. 内层循环：读取文件的每一行
    # while read 用来逐行读取文件内容，chr_col是第一列（染色体），pos是第二列（位置）
    while read -r chr_col pos; do
      # 计算区间的起始和结束位置
      a=$((pos - 5))
      b=$((pos + 5))
      
      # 拼接samtools需要的区域参数
      REGION="chr${c}:${a}-${b}"
      
      # 定义输出文件名，避免不同区间的结果互相覆盖，比如加个pos后缀
      OUTPUT_FILE="mpileup_chr${c}_${TISSUE}_pos${pos}.txt"
      
      # 执行samtools mpileup命令，把结果输出到指定文件
      echo "正在运行：samtools mpileup -f $REFERENCE File1 File2 -r $REGION > $OUTPUT_FILE"
      samtools mpileup -f "$REFERENCE" File1 File2 -r "$REGION" > "$OUTPUT_FILE"
      
    done < "$INPUT_FILE"  # 这里指定读取的输入文件
    
  done
done

echo "所有任务完成！"

脚本关键部分解释 #

• 变量定义：先把固定的内容（参考基因组路径、组织列表、染色体范围）定义成变量，方便后续修改
• 外层循环：先遍历组织，再遍历染色体，拼接出每个要处理的文件名，并且先检查文件是否存在，避免找不到文件的错误
• 内层循环：用while read -r chr_col pos逐行读取文件，chr_col对应第一列（染色体），pos对应第二列（位置）
• 区间计算：用$(( ))做算术运算，计算x-5和x+5
• 输出文件命名：给每个区间的结果单独命名，比如mpileup_chr1_colon_pos12345.txt，这样不会覆盖之前的结果
• 命令执行：最后运行samtools mpileup，把结果重定向到输出文件里

使用前要做的修改 #

1. 替换参考基因组路径：把REFERENCE="/path/to/your/reference/genome.fa"改成你实际的参考基因组文件路径，比如/data/hg38/hg38.fa
2. 确认File1和File2的路径：如果File1和File2不在当前目录，要写绝对路径，比如/data/File1.bam和/data/File2.bam
3. 检查文件名是否匹配：确保你的文件名和脚本里的INPUT_FILE拼接规则完全一致，如果有差异，要调整拼接的字符串

怎么运行脚本 #

1. 把上面的脚本复制到一个文件里，比如run_mpileup.sh
2. 给脚本添加执行权限：chmod +x run_mpileup.sh
3. 在终端里运行脚本：./run_mpileup.sh

小提示 #

• 如果你的文件很大，这个循环可能会运行比较久，可以考虑在后台运行：nohup ./run_mpileup.sh > run_log.txt 2>&1 &，这样即使终端关闭，脚本也会继续跑，日志会存在run_log.txt里
• 如果运行时提示samtools找不到，说明samtools不在你的PATH里，可以用绝对路径，比如/usr/bin/samtools或者你安装samtools的路径

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