WSL环境下Python调用bcftools+awk无输出问题求助

1. awk计算无输出的常见原因 #

• 字段匹配错误：先手动执行bcftools query -f '%QUAL\t%INFO/DP\n' input.vcf，确认输出的列序是否和awk脚本中使用的$1（QUAL）、$2（DP）对应，避免列位置搞错导致计算失效。
• 缺失值未处理：VCF中QUAL或DP字段常出现.（缺失值），awk会将其视为空字符串，直接累加会导致计算异常。需在awk中加入判断逻辑，仅处理有效值：

  bcftools query -f '%QUAL\t%INFO/DP\n' input.vcf | awk '{if ($1 != "." && $2 != ".") {sum_qual += $1; sum_dp += $2; count++}} END {if (count > 0) print sum_qual/count, sum_dp/count}'
  

• awk脚本逻辑疏漏：忘记初始化累加变量（sum_qual、count等），或者未在END块中打印结果，都会导致无输出。确保脚本包含完整的累加+最终计算逻辑。
• bcftools输出为空：如果VCF经过过滤后无符合条件的变异，bcftools没有输出内容，awk自然也不会返回结果，可手动执行bcftools命令验证输出。

2. Python、WSL与bash的引号/转义问题 #

这是此类场景的高频问题，主要因subprocess调用复杂bash命令时，引号嵌套解析错误导致。

错误示例（易触发问题）： #

import subprocess
cmd = "bcftools query -f '%QUAL\t%INFO/DP\n' input.vcf | awk '{sum += $1; count++} END {print sum/count}'"
result = subprocess.run(cmd, shell=True, capture_output=True, text=True)

Python字符串中的单引号会与awk脚本的单引号冲突，导致bash解析命令失败。

解决方案： #

• 使用raw字符串+双引号包裹awk脚本：

  cmd = r"""bcftools query -f '%QUAL\t%INFO/DP\n' input.vcf | awk "{sum_qual += $1; sum_dp += $2; count++} END {if (count > 0) print sum_qual/count, sum_dp/count}" """
  result = subprocess.run(cmd, shell=True, capture_output=True, text=True, executable="/bin/bash")
  

  注意指定executable="/bin/bash"，避免WSL默认用sh导致语法兼容问题。
• 更安全的无shell模式：拆分命令为列表形式调用，避免引号转义：

  import subprocess
  
  # 提取QUAL和DP字段
  bcf_cmd = ["bcftools", "query", "-f", "%QUAL\t%INFO/DP\n", "input.vcf"]
  bcf_output = subprocess.run(bcf_cmd, capture_output=True, text=True).stdout
  
  # awk计算平均值
  awk_cmd = ["awk", "{if ($1 != \".\" && $2 != \".\") {sum_qual += $1; sum_dp += $2; count++}} END {if (count > 0) print sum_qual/count, sum_dp/count}"]
  awk_result = subprocess.run(awk_cmd, input=bcf_output, capture_output=True, text=True).stdout.strip()
  

  这种方式无需处理复杂引号，还能避免shell注入风险。

3. %INFO/DP是否为正确的深度字段？ #

• 如果目标是整个变异位点的总测序深度（所有样本在该位点的reads数总和），%INFO/DP是正确的，它是VCF INFO区域的标准字段，记录位点的整体深度。
• 如果目标是单个样本的测序深度，则需使用%FORMAT/DP，同时指定样本（如-s sample_name），命令示例：

  bcftools query -s sample1 -f '%QUAL\t%FORMAT/DP\n' input.vcf
  

  注意区分INFO区域的位点总深度和FORMAT区域的单样本深度，避免混淆。

内容的提问来源于stack exchange，提问作者Rahul Kumar Singh