如何从BAM文件中筛选插入长度大于指定阈值的测序reads

核心筛选逻辑：CIGAR字段中I标记代表read相对于参考基因组存在的插入序列，本次筛选目标为携带单段长度大于阈值的I片段的reads。注意不要和BAM中TLEN字段标记的双端插入片段长度混淆，也不要将S（软剪切）、H（硬剪切）的未比对序列误计为插入。

方案1：Pysam脚本（推荐，效率最高，适合大BAM文件） #

依赖提前安装pysam与samtools，脚本直接读写BAM二进制文件，不需要生成中间文本，处理速度快。

import pysam

# 可自行修改配置参数
input_bam_path = "input.bam"
output_bam_path = "insert_over_50bp.bam"
min_insert_threshold = 50
CIGAR_I_OP = 1  # pysam中CIGAR插入操作I对应的固定编码

with pysam.AlignmentFile(input_bam_path, "rb") as in_bam, \
     pysam.AlignmentFile(output_bam_path, "wb", header=in_bam.header) as out_bam:
    for read in in_bam:
        # 跳过未比对、无有效CIGAR的read
        if read.is_unmapped or read.cigartuples is None:
            continue
        keep_read = False
        # 遍历CIGAR操作，存在单段插入长度超阈值直接标记保留
        for op, op_len in read.cigartuples:
            if op == CIGAR_I_OP and op_len > min_insert_threshold:
                keep_read = True
                break
        if keep_read:
            out_bam.write(read)

# 输出文件建索引，方便后续下游分析
pysam.index(output_bam_path)

脚本直接执行python filter_insert.py即可，对你给出的测试用例处理结果完全符合预期：

• Read1（CIGAR: 2M1I89M53I2M）：解析到两段插入，长度分别为1bp、53bp，53bp超过50bp阈值，被保留
• Read2（CIGAR: 2M1I144M）：仅存在1段1bp的插入，小于阈值，被过滤

方案2：Samtools+Awk一行流（适合临时快速处理，无需Python依赖） #

不需要额外写脚本，直接通过命令行管道处理，缺点是逐行解析文本，大BAM处理速度低于方案1：

samtools view -h input.bam | awk '
BEGIN {OFS="\t"; min_len=50}
# BAM头行直接输出
/^@/ {print; next}
{
    cigar_str = $6
    keep = 0
    # 逐段匹配CIGAR中的插入片段，判断长度
    while(match(cigar_str, /([0-9]+)I/, match_res)) {
        if (match_res[1] > min_len) {
            keep = 1
            break
        }
        cigar_str = substr(cigar_str, RSTART+RLENGTH)
    }
    if (keep) print
}' | samtools sort -o insert_over_50bp.bam - && samtools index insert_over_50bp.bam

• 不要混淆插入定义：CIGAR中的I是read相对参考的碱基插入，和双端测序的文库插入片段大小（SAM第9列TLEN字段）完全是两个概念，不要用TLEN做筛选。
• 不要误算其他CIGAR操作：S软剪切、H硬剪切是read两端未比对上的序列，不属于参考基因组上的插入片段，筛选时不要计入长度。
• 如果需求是read上所有插入片段的总长度大于阈值，只要把脚本/命令里的判断逻辑从「单段I长度超阈值即保留」改成「累加所有I片段长度，总和超阈值即保留」即可。

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