Linux系统下用Python3过滤遗传学.psl格式文件的技术求助

嘿，我来帮你搞定这个PSL文件过滤的问题！在Linux上用Python3处理遗传学的PSL文件其实挺直接的，尤其是处理那些带逗号分隔值的单元格。先给你一套通用的解决方案，你可以根据自己的具体筛选条件调整。

首先，先明确你的示例文件结构：这是制表符分隔的文件，部分列（比如第6、7列，索引从0开始算的话是第5、6列）包含逗号分隔的数值，有的还带末尾的逗号（比如8,9,）。

示例输入文件 (input.psl) #

1 2 3 x read1 8,9, 2001,2002,
1 2 3 mt read2 8,9,10 3001,3002,3003
1 2 3 9 read3 8,9,10,11 4001,4002,4003,4004
1 2 3 9 read4 8,9,10,11 4001,4002,4003,4004

通用过滤脚本 #

下面的脚本会逐行读取PSL文件，处理逗号分隔的字段，并根据你定义的条件筛选行：

#!/usr/bin/env python3

def meets_target_criteria(row_fields):
    """定义你的筛选条件，返回True则保留该行"""
    # 示例条件1: 第4列（参考序列ID，索引3）为"9"
    ref_id = row_fields[3]
    if ref_id != "9":
        return False
    
    # 处理第6列（read位置，索引5）的逗号分隔值，过滤空字符串
    read_positions = [pos.strip() for pos in row_fields[5].split(',') if pos.strip()]
    # 示例条件2: 位置数量不少于4个
    if len(read_positions) < 4:
        return False
    
    # 可以添加更多自定义条件，比如检查坐标范围：
    # coords = [int(c) for c in row_fields[6].split(',') if c.strip()]
    # if not any(4000 <= c <= 5000 for c in coords):
    #     return False
    
    return True

# 输入输出文件路径
input_path = "input.psl"
output_path = "filtered_targets.psl"

# 处理文件
with open(input_path, 'r') as infile, open(output_path, 'w') as outfile:
    for line in infile:
        cleaned_line = line.strip()
        if not cleaned_line:
            continue  # 跳过空行
        
        # 按制表符分割行数据
        fields = cleaned_line.split('\t')
        # 检查是否符合条件，符合则写入输出文件
        if meets_target_criteria(fields):
            outfile.write(cleaned_line + '\n')

print(f"筛选完成！结果已保存到 {output_path}")

使用说明 #

• 修改筛选条件：你只需要修改meets_target_criteria函数里的逻辑，比如：
  • 筛选特定的read名称（第5列，索引4）
  • 检查坐标列（第7列，索引6）是否落在你的目标区间内
  • 验证逗号分隔值的数量或具体数值
• 运行脚本：在Linux终端里，先给脚本加执行权限：

  chmod +x filter_psl.py
  

  然后直接运行：

  ./filter_psl.py
  

• 处理特殊情况：脚本已经处理了末尾带逗号的字段（比如8,9,会被转成['8','9']），避免空元素干扰判断。

示例调整：筛选坐标在4000-5000之间的行 #

如果你的目标是保留第7列（索引6）中存在坐标在4000到5000之间的行，可以把条件函数改成这样：

def meets_target_criteria(row_fields):
    # 处理坐标列，转成整数
    coords = [int(c.strip()) for c in row_fields[6].split(',') if c.strip()]
    # 检查是否有坐标在目标区间内
    return any(4000 <= c <= 5000 for c in coords)

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