文章目录
- 一、准备材料
- 二、提取外显子区间为BED文件
- 1. 提取GTF中exon为BED
- 三、用bedtools提取外显子fasta
- 四、后续拼接外显子为ORF序列
- 五、流程总结
一、准备材料
- 基因组fasta（如：genome.fa）
- RiboCode生成的GTF文件（如：ribocode.gtf，含ORF的exon注释）
二、提取外显子区间为BED文件
RiboCode的GTF文件每个ORF有多条exon记录。我们需要把这些exon条目提取出来，转为BED格式。
1. 提取GTF中exon为BED
这里我们使用RiboCode生成的GTF格式：
NC_000001.11 RiboCode exon 17233 17310 . - . orf_id "WASH7P_17310_16745_98"; ...
Python脚本示例：
# gtf_to_bed.py
with open("ribocode.gtf") as gtf, open("orf_exons.bed", "w") as bed:for line in gtf:if line.startswith("#"): continuefields = line.strip().split("\t")if fields[2] != "exon": continuechrom = fields[0]start = str(int(fields[3]) - 1) # GTF是1-based闭区间，BED是0-based左闭右开end = fields[4]strand = fields[6]attrs = {kv.split(' ')[0]: kv.split(' ')[1].replace('"','').replace(';','') for kv in fields[8].split('; ') if kv}orf_id = attrs.get("orf_id", "NA")bed.write(f"{chrom}\t{start}\t{end}\t{orf_id}\t.\t{strand}\n")
运行：
python gtf_to_bed.py
三、用bedtools提取外显子fasta
确保你已安装bedtools。
bedtools getfasta -fi genome.fa -bed orf_exons.bed -s -name -fo orf_exons.fa
-fi genome.fa：输入基因组fasta-bed orf_exons.bed：外显子区间-s：考虑链向-name：fasta头用bed第四列（即orf_id）
这样会得到每个外显子的fasta序列，头部格式如：
>WASH7P_17310_16745_98::NC_000001.11:17232-17310(-)
ATGCTG...
四、后续拼接外显子为ORF序列
用你之前优化的merge_exon_seq.py脚本即可，将同一orf_id的外显子按顺序拼接，得到完整的ORF核酸序列。
五、流程总结
- GTF转BED：提取exon条目，转为bed格式。
- bedtools getfasta：用bed文件和基因组fasta提取外显子序列。
- 拼接外显子：按orf_id分组拼接，得到完整ORF序列。
