生信新手实战指南从PASA注释结果高效提取gff与fasta文件刚接触基因组注释的研究者常会遇到这样的困境费尽周折获得的PASA注释结果却不知如何提取所需的gff和fasta文件。本文将手把手教你如何从原始数据中提取关键信息并转化为下游分析所需的标准化格式。1. 理解PASA输出结构与文件用途PASAProgram to Assemble Spliced Alignments是基因组注释流程中的关键环节其输出包含多种文件类型每种都有特定用途gff3文件记录基因结构注释信息包含基因、外显子、CDS等特征的位置fasta文件包含基因组序列、蛋白质序列或CDS序列lens文件记录染色体或contig的长度信息关键点不同版本的PASA可能生成略有差异的文件命名格式但核心内容结构相似。典型的PASA输出文件名可能包含类似gene_structures_post_PASA_updates.[数字].gff3的模式。2. 准备工作与环境配置在开始提取前需要确保工作环境准备就绪# 检查Perl环境许多PASA工具依赖Perl perl -v # 安装必要的生物信息学工具 conda install -c bioconda evm-utils bedtools提示建议使用conda或mamba管理生物信息学软件环境避免依赖冲突需要准备的基础文件包括PASA生成的gff3注释文件对应的基因组fasta文件可选EVMEvidenceModeler工具包3. 使用EVM工具提取序列文件EVM工具包中的gff3_file_to_proteins.pl脚本是处理PASA输出的利器可提取多种序列格式# 提取蛋白质序列 evidencemodeler-1.1.1/EvmUtils/gff3_file_to_proteins.pl \ input.gene_structures_post_PASA_updates.gff3 \ genome.fasta prot output.prot.fasta # 提取CDS序列 evidencemodeler-1.1.1/EvmUtils/gff3_file_to_proteins.pl \ input.gene_structures_post_PASA_updates.gff3 \ genome.fasta CDS output.cds.fasta # 提取基因序列包含内含子 evidencemodeler-1.1.1/EvmUtils/gff3_file_to_proteins.pl \ input.gene_structures_post_PASA_updates.gff3 \ genome.fasta gene output.gene.fasta参数说明第一个参数输入的gff3文件第二个参数基因组fasta文件第三个参数序列类型prot/CDS/gene输出重定向到指定文件4. 处理与验证提取结果提取完成后应对文件进行基本验证# 检查fasta文件格式 head -n 4 output.prot.fasta # 统计序列数量 grep -c ^ output.prot.fasta # 检查gff3文件完整性 grep -v ^# input.gff3 | head常见问题处理问题现象可能原因解决方案空输出文件gff3与基因组序列ID不匹配检查并统一序列ID命名部分序列缺失注释不完整检查PASA运行日志序列异常短提取参数错误确认使用了正确的序列类型参数5. 生成染色体长度文件lens染色体长度信息对许多分析至关重要可通过多种方式生成方法一使用基因组fasta文件# 使用seqkit工具 seqkit fx2tab -n -l genome.fasta | awk {print $1\t$2} genome.lens # 使用bioawk bioawk -c fastx {print $name\tlength($seq)} genome.fasta genome.lens方法二基于gff3注释# 提取每个序列区域的最大位置 awk $3 gene {if($5 len[$1]) len[$1]$5} END {for(id in len) print id\tlen[id]} input.gff3 estimated.lens注意基于gff3估计的长度可能略小于实际基因组长度建议优先使用fasta直接计算6. 进阶处理与格式转换对于需要进一步清洗或转换格式的情况可以考虑以下工具使用GeneClean工具包# 安装 pip install GeneClear # 准备配置文件 GeneClear -getpasa ? run.conf # 编辑run.conf指定输入文件路径 vi run.conf # 运行处理 GeneClear -getpasa run.conf手动处理常见格式问题序列ID统一化# 修改蛋白质fasta头部的ID格式 sed s/\.CDS$// input.prot.fasta clean.prot.fasta染色体命名标准化# 将EVM_01改为chr01 sed s/EVM_/chr/g input.gff3 standard.gff37. 下游分析准备检查清单在进入后续分析前建议确认以下文件已正确生成并验证[ ] 标准化的gff3注释文件[ ] 基因组fasta文件[ ] 蛋白质序列fasta[ ] CDS序列fasta[ ] 基因序列fasta可选[ ] 染色体长度文件lens文件组织结构示例project/ ├── annotations/ │ ├── genome.gff3 │ ├── genome.fasta │ ├── proteins.fasta │ ├── cds.fasta │ └── genome.lens ├── scripts/ └── results/8. 实际案例分析拟南芥注释处理以拟南芥数据为例演示完整处理流程# 下载示例数据 wget ftp://ftp.ensemblgenomes.org/pub/plants/release-55/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.55.gff3.gz wget ftp://ftp.ensemblgenomes.org/pub/plants/release-55/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz # 解压 gunzip *.gz # 提取蛋白质序列 gff3_file_to_proteins.pl \ Arabidopsis_thaliana.TAIR10.55.gff3 \ Arabidopsis_thaliana.TAIR10.dna.toplevel.fa prot athaliana.prot.fasta # 生成长度文件 seqkit fx2tab -n -l Arabidopsis_thaliana.TAIR10.dna.toplevel.fa | awk {print $1\t$2} athaliana.lens处理过程中的典型挑战大型gff3文件处理时内存不足 → 使用流式处理或分块处理ID命名不一致 → 提前统一命名规范特殊字符问题 → 清洗非ASCII字符9. 自动化脚本示例为提高效率可创建自动化处理脚本extract_pasa.sh#!/bin/bash # 参数检查 if [ $# -ne 2 ]; then echo Usage: $0 input.gff3 genome.fasta exit 1 fi # 定义输出文件名前缀 prefix$(basename $1 .gff3) # 提取各类序列 gff3_file_to_proteins.pl $1 $2 prot ${prefix}.prot.fasta gff3_file_to_proteins.pl $1 $2 CDS ${prefix}.cds.fasta gff3_file_to_proteins.pl $1 $2 gene ${prefix}.gene.fasta # 生成长度文件 seqkit fx2tab -n -l $2 | awk {print $1\t$2} ${prefix}.lens # 验证输出 echo 处理完成生成文件 ls -lh ${prefix}.* echo 蛋白质序列数量: $(grep -c ^ ${prefix}.prot.fasta)使用方式chmod x extract_pasa.sh ./extract_pasa.sh annotation.gff3 genome.fasta10. 常见问题排查指南遇到问题时可按照以下步骤排查文件无法生成检查输入文件路径是否正确验证文件权限检查磁盘空间输出文件内容异常# 检查gff3与fasta的序列ID是否匹配 head -n 1 input.gff3 grep ^ genome.fasta | head -n 1性能优化建议对大基因组使用pigz替代gzip考虑使用awk或bioawk处理大型文本对重复操作使用并行处理性能对比表工具/方法适用场景内存占用处理速度标准Perl脚本小中型数据集中等中等awk/bioawk大型数据集低快Python脚本复杂转换高取决于实现11. 扩展应用多物种处理流程当需要处理多个物种时可采用批量化处理# 创建物种列表 cat species.list EOF species1 species2 species3 EOF # 批量处理 while read sp; do gff3_file_to_proteins.pl \ ${sp}.annotation.gff3 \ ${sp}.genome.fasta prot ${sp}.prot.fasta done species.list对于超大规模分析建议使用工作流管理系统如Nextflow或Snakemake。以下是一个简单的Snakemake规则示例rule extract_sequences: input: gff3 {species}.annotation.gff3, fasta {species}.genome.fasta output: prot {species}.prot.fasta, cds {species}.cds.fasta shell: gff3_file_to_proteins.pl {input.gff3} {input.fasta} prot {output.prot} gff3_file_to_proteins.pl {input.gff3} {input.fasta} CDS {output.cds} 12. 质量评估与结果验证为确保提取结果的可靠性建议进行以下检查序列完整性验证# 检查蛋白质序列是否包含终止符 grep -v athaliana.prot.fasta | grep -c *$ # 检查CDS长度是否为3的倍数 bioawk -c fastx {if(length($seq)%3 ! 0) print $name} athaliana.cds.fastagff3文件验证# 使用AGAT工具验证gff3结构 agat_sp_validate_gff.pl --gff input.gff3一致性检查# 比较gff3中的基因数与fasta中的序列数 grep -c $\tgene\t input.gff3 grep -c ^ output.prot.fasta差异分析表检查项预期结果异常可能原因基因数 vs 蛋白数相等注释不完整或提取错误CDS长度3的倍数注释错误或序列问题染色体长度与参考一致版本不匹配或组装问题