Python处理BAM文件argparse参数报错及列长度不匹配问题求助

报错根因 #

• 第一个TypeError: count_bam() missing 3 required positional arguments: 'bed_file', 'output', and 'genome_sizes'报错，是因为你定义了parseargs函数但从未调用，直接调用count_bam时未传入对应参数，和你后续加的全局硬编码变量无关。
• 第二个ValueError: Length mismatch: Expected axis has 1 elements, new values have 3 elements报错，和argparse完全无关，是你读取的bed文件不符合规范：要么存在空行、要么部分行只有1列、要么分隔符不是制表符，导致你取前三列切片后实际列数不足3，无法匹配chr start end三个列名。
• 你怀疑的argparse参数未注册问题不存在：代码中add_argument已经注册了4个参数，只是你没有调用parseargs()触发解析流程，所以参数不会生效。

修复方案 #

方案1：保留硬编码路径本地测试（无需使用argparse） #

删除冗余的argparse相关代码，增加bed文件格式校验逻辑，代码如下：

import pysam
import pandas as pd

def count_bam(bamfile, bed_file, output, genome_sizes, use_fast_count=True, verbose=True):
    reads = pysam.AlignmentFile(bamfile)
    read_chrs = set(reads.references)
    # 增加bed文件读取校验，指定分隔符，跳过空行
    bed_regions = pd.read_table(bed_file, header=None, sep="\t", skip_blank_lines=True)
    # 校验列数是否符合要求
    if bed_regions.shape[1] < 3:
        raise ValueError(f"输入bed文件列数不足3，当前仅{bed_regions.shape[1]}列，请检查文件格式")
    bed_regions = bed_regions[bed_regions.columns[:3]]
    bed_regions.columns = ["chr", "start", "end"]
    counts = [(reads.count(row.chr, row.start, row.end) if (row.chr in read_chrs) else 0) for _, row in bed_regions.iterrows()]
    bed_regions['count'] = counts
    bed_regions.to_csv(output, header=None, index=None, sep="\t")

# 硬编码参数直接调用
bamfile = "/mnt/d/Axiotl/ABC/Input/bam/ENCFF070PWH.bam"
bed_file = "/mnt/d/R/abc_3a_scaffold_hg19.bed"
output = "/mnt/d/Axiotl/ABC/Output/bam" 
genome_sizes = "/mnt/d/Axiotl/ABC/Working_Copy/ABC-Enhancer-Gene-Prediction-master/reference/chr_sizes"

count_bam(bamfile, bed_file, output, genome_sizes)

方案2：使用argparse命令行传参（无需硬编码路径） #

调用parseargs()函数触发参数解析，删除硬编码的路径变量，同时增加参数必填校验，代码如下：

import pysam
import pandas as pd
import argparse

def parseargs():
    class formatter(argparse.ArgumentDefaultsHelpFormatter, argparse.RawTextHelpFormatter):
        pass
    parser = argparse.ArgumentParser(description='统计BAM文件在指定bed区域的reads数',
                                     formatter_class=formatter)
    # 增加required=True标记参数为必填项
    parser.add_argument('--bamfile', required=True, help="输入BAM文件路径")
    parser.add_argument('--bed_file', required=True, help="输入bed区域文件路径")
    parser.add_argument('--genome_sizes', required=True, help="基因组染色体大小文件路径")
    parser.add_argument('--output', required=True, help="输出结果文件路径")
    args = parser.parse_args()
    print(args)
    return args

def count_bam(bamfile, bed_file, output, genome_sizes, use_fast_count=True, verbose=True):
    reads = pysam.AlignmentFile(bamfile)
    read_chrs = set(reads.references)
    bed_regions = pd.read_table(bed_file, header=None, sep="\t", skip_blank_lines=True)
    if bed_regions.shape[1] < 3:
        raise ValueError(f"输入bed文件列数不足3，当前仅{bed_regions.shape[1]}列，请检查文件格式")
    bed_regions = bed_regions[bed_regions.columns[:3]]
    bed_regions.columns = ["chr", "start", "end"]
    counts = [(reads.count(row.chr, row.start, row.end) if (row.chr in read_chrs) else 0) for _, row in bed_regions.iterrows()]
    bed_regions['count'] = counts
    bed_regions.to_csv(output, header=None, index=None, sep="\t")

if __name__ == "__main__":
    args = parseargs()
    count_bam(args.bamfile, args.bed_file, args.output, args.genome_sizes)

使用时直接在命令行执行以下命令即可：

python 你的脚本文件名.py --bamfile <BAM文件路径> --bed_file <bed文件路径> --genome_sizes <染色体大小文件路径> --output <输出文件路径>

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