对序列长度进行过滤

D:\1CAAS\Lab\songqianlin\Cas新蛋白\Cas12\cas12_lmnopq_fasta\重复项去除\序列长度.py

点击查看代码

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

"""
功能：
1. 统计FASTA文件中序列总数、最短、最长、平均长度
2. 过滤掉长度 <300 aa 的序列（默认300，可自定义）
3. 输出保留序列FASTA（单行序列）
4. 输出被过滤序列FASTA（单行序列）
5. 输出被过滤序列Excel表（ID、描述、长度、序列）
6. 直接打印统计结果
"""

from openpyxl import Workbook

# =========================
# 输入输出路径（直接修改这里）
# =========================
input_file = r"D:\1CAAS\tree\20260510_PLV_CM_HG_5AHG_MDM_TnpB\PLV_CM_HG_5AHG_MDM_TnpB.fasta"

# 保留序列
output_file = r"D:\1CAAS\tree\20260510_PLV_CM_HG_5AHG_MDM_TnpB\PLV_CM_HG_5AHG_MDM_TnpB-t300.fasta"

# 被过滤序列FASTA
dropped_fasta = r"D:\1CAAS\tree\20260510_PLV_CM_HG_5AHG_MDM_TnpB\PLV_CM_HG_5AHG_MDM_TnpB-below300.fasta"

# 被过滤序列Excel
dropped_excel = r"D:\1CAAS\tree\20260510_PLV_CM_HG_5AHG_MDM_TnpB\PLV_CM_HG_5AHG_MDM_TnpB-below300.xlsx"

min_length = 300  # 最小保留长度

# =========================
# 初始化
# =========================
lengths = []
records = []

header = None
seq_parts = []

# =========================
# 读取FASTA
# =========================
with open(input_file, "r") as f:
    for line in f:
        line = line.strip()

        if not line:
            continue

        if line.startswith(">"):
            if header is not None:
                sequence = "".join(seq_parts)
                seq_len = len(sequence)

                lengths.append(seq_len)
                records.append((header, sequence, seq_len))

            header = line
            seq_parts = []

        else:
            seq_parts.append(line)

    # 最后一条
    if header is not None:
        sequence = "".join(seq_parts)
        seq_len = len(sequence)

        lengths.append(seq_len)
        records.append((header, sequence, seq_len))

# =========================
# 基本统计
# =========================
total_count = len(lengths)

if total_count == 0:
    print("输入文件中没有检测到序列。")
    exit()

min_len = min(lengths)
max_len = max(lengths)
avg_len = sum(lengths) / total_count

# =========================
# 过滤输出
# =========================

keep_count = 0
drop_count = 0
dropped_records = []

with open(output_file, "w") as keep_f, open(dropped_fasta, "w") as drop_f:

    for header, sequence, seq_len in records:

        if seq_len >= min_length:
            # 保留序列
            keep_f.write(f"{header}\n")
            keep_f.write(f"{sequence}\n")

            keep_count += 1

        else:
            # 删除序列
            drop_f.write(f"{header}\n")
            drop_f.write(f"{sequence}\n")

            dropped_records.append((header, sequence, seq_len))
            drop_count += 1

# =========================
# 输出被过滤序列Excel
# =========================
wb = Workbook()
ws = wb.active
ws.title = "Below_Min_Length"

# 表头
ws.append(["Sequence_ID", "Description", "Length", "Sequence"])

for header, sequence, seq_len in dropped_records:
    full_header = header[1:]  # 去掉 >
    seq_id = full_header.split()[0]

    # 描述（如果有）
    description = full_header[len(seq_id):].strip() if len(full_header.split()) > 1 else ""

    ws.append([seq_id, description, seq_len, sequence])

# 自动调整列宽
for col in ws.columns:
    max_length = 0
    col_letter = col[0].column_letter

    for cell in col:
        try:
            if cell.value:
                max_length = max(max_length, len(str(cell.value)))
        except:
            pass

    ws.column_dimensions[col_letter].width = min(max_length + 2, 100)

wb.save(dropped_excel)

# =========================
# 输出结果
# =========================
print("====== FASTA统计结果 ======")
print(f"输入文件: {input_file}")
print(f"总序列数: {total_count}")
print(f"最短序列长度: {min_len} aa")
print(f"最长序列长度: {max_len} aa")
print(f"平均序列长度: {avg_len:.2f} aa")

print(f"\n====== 过滤结果（>={min_length} aa） ======")
print(f"保留序列数: {keep_count}")
print(f"删除序列数: {drop_count}")

print("\n====== 输出文件 ======")
print(f"保留序列FASTA: {output_file}")
print(f"删除序列FASTA: {dropped_fasta}")
print(f"删除序列Excel: {dropped_excel}")