利用TSV/CSV表格拆分大型多序列FASTA文件用于比对分析

我来给你两个实用的方案，一个是灵活可定制的Python脚本，另一个是适合命令行老手的awk命令组合，都能完美实现你按TSV基因座分组拆分大型FASTA文件的需求。

方案一：Python脚本（推荐，易扩展） #

这个脚本会先解析你的TSV文件，建立序列头到基因座的映射，然后逐行读取FASTA文件，将每个序列写入对应基因座的独立FASTA文件。脚本处理逻辑清晰，方便你根据后续需求调整（比如处理未匹配的序列）。

import sys
from collections import defaultdict

def main(fasta_file, tsv_file):
    # 解析TSV，构建「基因座→序列头集合」的映射
    locus_to_headers = defaultdict(set)
    with open(tsv_file, 'r') as tsv:
        for line in tsv:
            line = line.strip()
            if not line:
                continue
            parts = line.split('\t')
            locus = parts[0]
            # 遍历该行后续所有序列头（跳过空字段）
            for header in parts[1:]:
                if header:
                    locus_to_headers[locus].add(header)
    
    # 反向构建「序列头→基因座」的映射，方便快速查找
    header_to_locus = {}
    for locus, headers in locus_to_headers.items():
        for header in headers:
            header_to_locus[header] = locus
    
    # 处理FASTA文件，写入对应基因座的文件
    file_handles = {}  # 保存已打开的文件句柄，避免重复IO操作
    current_header = None
    current_seq = []

    with open(fasta_file, 'r') as fasta:
        for line in fasta:
            line = line.strip()
            if not line:
                continue
            # 遇到序列头行，先处理上一个序列
            if line.startswith('>'):
                if current_header is not None:
                    locus = header_to_locus.get(current_header)
                    if locus:
                        # 打开对应基因座的文件（如果还没打开）
                        if locus not in file_handles:
                            file_handles[locus] = open(f"{locus}.fa", 'w')
                        # 写入序列
                        fh = file_handles[locus]
                        fh.write(f">{current_header}\n")
                        fh.write('\n'.join(current_seq) + '\n')
                # 更新当前序列的头和内容
                current_header = line[1:]
                current_seq = []
            else:
                current_seq.append(line)
        # 处理最后一个序列
        if current_header is not None:
            locus = header_to_locus.get(current_header)
            if locus:
                if locus not in file_handles:
                    file_handles[locus] = open(f"{locus}.fa", 'w')
                fh = file_handles[locus]
                fh.write(f">{current_header}\n")
                fh.write('\n'.join(current_seq) + '\n')
    
    # 关闭所有打开的文件
    for fh in file_handles.values():
        fh.close()

if __name__ == "__main__":
    if len(sys.argv) != 3:
        print("用法: python split_fasta_by_locus.py <seq.fa> <基因座TSV文件>")
        sys.exit(1)
    main(sys.argv[1], sys.argv[2])

使用方法 #

1. 将上述代码保存为split_fasta_by_locus.py
2. 在命令行运行：

   python split_fasta_by_locus.py seq.fa your_locus_file.tsv
   

3. 运行完成后，当前目录会生成每个基因座对应的[基因座名称].fa文件，每个文件包含该基因座的所有序列。

方案二：awk命令组合（轻量快速） #

如果你熟悉Linux命令行，这个组合更简洁，不需要编写脚本，直接用awk处理：

步骤1：生成序列头到基因座的映射文件 #

首先把TSV转换成便于awk读取的映射格式（每行是「序列头\t基因座」）：

awk -F'\t' '{for(i=2;i<=NF;i++) if($i!="") print $i "\t" $1}' your_locus_file.tsv > header_locus.map

步骤2：拆分FASTA文件 #

用awk读取映射文件和FASTA文件，将序列写入对应基因座的文件：

awk 'BEGIN{FS="\t"; while((getline < "header_locus.map")>0) map[$1]=$2}
     /^>/ {current_head=substr($0,2); target_locus=map[current_head]; if(target_locus!="") out_file=target_locus".fa"}
     target_locus!="" {print >> out_file; close(out_file)}' seq.fa

说明 #

• close(out_file)是为了避免同时打开过多文件句柄，适合处理超大型FASTA文件
• 不在映射中的序列会被自动忽略

额外注意事项 #

• 如果你的TSV中有重复的序列头，两个方案都会自动去重，不会重复写入序列
• 如果你想保留未匹配的序列，可以修改Python脚本，添加一个unmatched.fa的文件句柄，将不在映射中的序列写入其中
• 两个方案都是逐行读取文件，不会将整个大型FASTA加载到内存，适合处理GB级别的文件

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