求助：基于origin.bed生成nucleotide.fa及protein.fa的可行方法

核心处理逻辑 #

bedtools getfasta无法直接实现按TranscriptID分组拼接的需求，因此需要自定义流程：

• 严格按TranscriptID对bed条目分组，同一TranscriptID的片段才进行拼接（即使属于同一GeneID，不同TranscriptID也需分开处理）
• 根据链方向排序片段：
  • 正链(+)：按起始坐标从小到大排序，匹配转录方向
  • 负链(-)：按起始坐标从大到小排序，符合负链转录的顺序
• 提取每个片段的基因组序列，负链片段需做反向互补处理
• 拼接同一TranscriptID的所有序列段，生成指定格式的FASTA条目（header包含染色体、合并后的坐标范围、TranscriptID）

具体实现（Python脚本） #

用Python结合pysam、pandas、biopython可以灵活实现需求，脚本如下：

import pandas as pd
import pysam
from Bio.Seq import Seq

# 读取origin.bed文件，指定列名
bed_df = pd.read_csv('origin.bed', sep='\t', header=None, 
                     names=['chr', 'start', 'end', 'transcript_id', 'strand', 'gene_id'])

# 加载参考基因组FASTA（需提前用samtools faidx建立索引）
genome_fasta = pysam.FastaFile('genome.fa')
output_fa = open('nucleotide.fa', 'w')

# 遍历每个TranscriptID分组
for transcript_id, group in bed_df.groupby('transcript_id'):
    strand = group.iloc[0]['strand']
    # 按链方向排序片段
    if strand == '+':
        sorted_group = group.sort_values('start', ascending=True)
    else:
        sorted_group = group.sort_values('start', ascending=False)
    
    # 生成FASTA header所需的坐标信息
    chr_name = sorted_group.iloc[0]['chr']
    min_start = sorted_group['start'].min()
    max_end = sorted_group['end'].max()
    
    # 拼接序列
    full_seq = ''
    for _, row in sorted_group.iterrows():
        # pysam采用0-based闭区间，bed的end是开区间，因此end需减1
        seq = genome_fasta.fetch(row['chr'], row['start'], row['end'])
        if strand == '-':
            # 负链序列做反向互补
            seq = str(Seq(seq).reverse_complement())
        full_seq += seq
    
    # 写入FASTA文件，按每行80碱基换行符合规范
    output_fa.write(f">{chr_name}:{min_start}-{max_end}({transcript_id})\n")
    for i in range(0, len(full_seq), 80):
        output_fa.write(f"{full_seq[i:i+80]}\n")

# 关闭文件句柄
genome_fasta.close()
output_fa.close()

脚本说明 #

1. 依赖安装：执行以下命令安装所需库

   pip install pandas pysam biopython
   

2. 预处理要求：参考基因组FASTA需提前用samtools faidx genome.fa建立索引
3. 关键细节：
   • 负链片段先按start从大到小排序，再做反向互补，保证拼接后的序列符合转录逻辑
   • FASTA header严格匹配需求格式，序列自动按80碱基换行

替代命令行组合方案（可选） #

如果偏好命令行工具，可通过awk+bedtools+samtools组合实现，但步骤较繁琐：

1. 用awk按TranscriptID分组，输出每个分组的排序后bed条目
2. 对每个分组的bed条目，用bedtools getfasta -s（-s参数自动处理链特异性，负链反向互补）提取序列
3. 用awk拼接同一TranscriptID的序列段，生成目标FASTA

该方案需编写循环处理每个分组，灵活性不如Python脚本，仅适合熟悉命令行的用户。

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