如何在Linux服务器用Python脚本提取特定格式的序列ID？

问题分析 #

你需要提取以TRINITY开头、.p1结尾的序列ID，之前的尝试存在两处核心问题：

• awk '{print$1}'会输出每行第一个字段，但如果文件包含非ID行（比如蛋白序列行），就会混入无关内容
• 你的Python脚本存在多处逻辑错误，导致无法匹配目标ID

脚本错误点 #

你的Python脚本有以下关键问题：

1. 判断条件错误："\.p1" in file是检查字符串是否在文件对象中，而非当前读取的行myline，应改为".p1" in myline
2. 文件打开模式错误：每次匹配到就用w模式打开文件，会覆盖之前写入的内容，需用a追加模式或提前打开输出文件
3. 写入内容缺失：new_file.write()未传入要写入的内容，需提取目标ID
4. 提示逻辑错误：循环中每一行不匹配就打印提示，会重复输出大量信息，应在遍历完整个文件后判断是否有匹配结果再输出

修正方案 #

方案1：改进awk命令 #

直接用awk过滤目标ID行，只提取纯ID内容：

awk '/^>TRINITY.*\.p1$/ {print substr($1,2)}' filename.cdhit > seq_id.fasta

• /^>TRINITY.*\.p1$/：匹配以>开头、紧跟TRINITY、最后以.p1结尾的行
• substr($1,2)：去掉开头的>，只保留ID部分

方案2：修正后的Python脚本 #

import re

file_path = '/var2/user/de_novo/data/transdecoder_dir/Trinity.fasta.transdecoder.pep.cdhit'
new_file_path = '/var2/user/de_novo/data/transdecoder_dir/seqID.fasta'
match_count = 0

# 提前打开输出文件，避免重复打开覆盖内容
with open(new_file_path, 'w') as new_file:
    with open(file_path, 'rt') as file:
        for line in file:
            line = line.strip()
            # 正则精确匹配目标ID行
            if re.match(r'^>TRINITY.*\.p1$', line):
                # 提取ID并写入
                seq_id = line[1:]
                new_file.write(f"{seq_id}\n")
                match_count += 1

# 统一输出结果提示
if match_count == 0:
    print('No match found.')
else:
    print(f'成功提取{match_count}个序列ID，已保存到{new_file_path}')

• 用正则表达式精确匹配目标行，避免误匹配
• 提前打开输出文件，提升效率同时防止内容被覆盖
• 统计匹配数量，最后统一输出提示信息

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