使用seqkit对原始蛋白组序列进行批量长度筛选，剔除短片段序列，得到标准化高质量蛋白序列

一、说明

蛋白序列长度过滤：

1. 去除无效噪声序列：基因组或转录组预测得到的原始蛋白文件中，会存在大量过短氨基酸序列，这类序列可能为错误ORF，随机翻译片段或残缺编码序列等等，直接参与分析会影响结果
2. 统一内外群筛选阈值：比如之前我在分析内群榧属物种时，在CD-HIT聚类阶段，已经过滤50aa以下的短序列，然而外群的蛋白序列还未经过筛选

一般筛选的阈值：过滤小于50aa的蛋白序列，保留≥50aa，目前在TransDecoder，OrfFinder，CD-HIT等默认为50aa过滤

二、环境配置

conda install -c bioconda seqkit biopython

三、核心代码

from pathlib import Path
import subprocess

SPECIES_LIST = ["specie1", "specie2", "specie3"]

DATA_DIR = Path("/data/data")
RESULT_DIR = Path("/data/result")
FINAL_DIR = RESULT_DIR / "specie_proteome"

# 蛋白序列过滤
def filter_protein(specie):
    raw_pep = DATA_DIR / f"{specie}.chr.pep"
    out_dir = FINAL_DIR / f"{specie}"
    out_dir.mkdir(parents=True, exist_ok=True)
    final_pep = out_dir / f"{specie}.filtered.pep.fa"

    # Seqkit过滤：保留50aa及以上蛋白序列
    print("[INFO]过滤50aa以下短序列...")
    filter_cmd = ["seqkit", "seq", "-m", "50", str(raw_pep)]
    with open(final_pep, "w", encoding="utf-8") as f:
        subprocess.run(filter_cmd, stdout=f, check=True)
    print(f"[DONE]{specie} 蛋白序列已导出：{final_pep}\n")

if __name__ == "__main__":
    for sp in SPECIES_LIST:
        filter_protein(sp)

四、结果输出

/data/result/specie_proteome/{specie}/{specie}.filtered.pep.fa

【金山文档 | WPS云文档】 seqkit蛋白序列长度过滤 https://www.kdocs.cn/l/ckEXaChLtpt4