By liyupeng, 22 September, 2025
Forums
Researches

基因功能注释：eggnog & nr篇

介绍

基因功能注释是基因组研究中的关键环节，有助于识别和解析蛋白质编码基因及其生物学功能。通过将预测得到的基因序列与权威数据库（如 NCBI nr、Swiss-Prot、eggNOG 等）进行同源比对，可以推测基因的潜在功能。功能注释不仅能够揭示基因在生物学过程、分子功能和细胞组分等方面的归属（Gene Ontology，GO），还可借助 KEGG 等数据库将基因与代谢通路相联系。系统化的注释结果为深入理解物种的生理特性、进化机制以及特有的生物学功能奠定了基础。

eggnog和nr数据库已经由喵桑和远兄下载到了92服务器上

eggnog的参考：https://www.kdocs.cn/l/crkLpQvsdFCG

不过，在使用的时候，考虑到92上不多的cpu数量以及需要比对蛋白序列的数量，我是把这个eggnog数据库打包传输到工作站上然后运行的，反正数据量也不是特别大。

至于nr数据库，这个数据库解压后有1T多，虽然不知道为什么这么大，咱们不敢问，咱们也不敢说，只能老老实实地在92服务器上跑了。

这里再补充两个线上工具，eggnog和nr的，前者经常挂掉，所以没法指望

eggnog-mapper：eggNOG-mapper

nr（应该）：Protein BLAST: search protein databases using a protein query

除了eggnog-mapper以外，其他的数据库限制上传蛋白序列的数量和文件大小，大约在2mb，按需使用。

还有就是，nr数据库比较大，所以比较花费时间。关于这个，实际上可以根据索引什么的选择库，什么都不选的话肯定是比对一整个数据库，选择了植物、陆生植物或者其他支选项的话，比对的数据量也会少一些。还有就是用dimond来比对了。

操作流程

eggnog-mapper注释

这里贴一个带注释和不带注释的命令行就可以了，输出文件时一个excel表，很直观

emapper.py \
  --data_dir /database/eggnog_data_dir \         # eggNOG 数据库路径
  -i Rhododendron_simsii_longest_pep.fa \        # 输入文件（蛋白序列 FASTA）
  --cpu 50 \                                     # 使用 50 个线程并行
  --mp_start_method forkserver \                 # 多进程启动方式设为 forkserver
  -o Rhododendron_simsii_longest_pep.eggnog.txt \# 输出文件前缀
  --output_dir ./ \                              # 输出目录
  --override \                                   # 覆盖已存在的输出文件
  -m diamond \                                   # 比对方式：diamond（比 BLAST 快）
  --dmnd_ignore_warnings \                       # 忽略 diamond 的警告
  --evalue 0.001 \                               # E-value 阈值 = 1e-3
  --score 60 \                                   # 最低比对得分阈值 = 60
  --pident 40 \                                  # 要求序列相似度 ≥ 40%
  --query_cover 20 \                             # 查询序列覆盖度 ≥ 20%
  --subject_cover 20 \                           # 目标序列覆盖度 ≥ 20%
  --itype proteins \                             # 输入类型：蛋白质序列
  --tax_scope 33090 \                            # 物种范围限制：植物 (taxid 33090)
  --target_orthologs all \                       # 输出所有目标直系同源基因
  --go_evidence non-electronic \                 # 只保留非电子推断的 GO 注释
  --pfam_realign none \                          # 不重新对齐 Pfam
  --report_orthologs \                           # 在输出里报告直系同源基因
  --decorate_gff yes \                           # 如果有 GFF 文件，写回功能注释
  --excel                                        # 额外输出 Excel 格式结果
  1. 运行命令
nohup emapper.py \
  -i female_hap2_rename_protein.fa \
  -o female_hap2_rename_protein.fa.eggnog.txt \
  --cpu 20 \
  --data_dir /database/eggnog_data_dir --mp_start_method forkserver --output_dir ./ --override -m diamond --dmnd_ignore_warnings --evalue 0.001 --score 60 --pident 40 --query_cover 20 --subject_cover 20 --itype proteins --tax_scope 33090 --target_orthologs all --go_evidence non-electronic --pfam_realign none --report_orthologs --decorate_gff yes --excel &

nr注释

  1. nr注释我是参考了微信公众号上的一片文章，大致的操作都相同，diamond比对，然后用一个脚本整理一下比对的结果；不过，关于结果的输出，也可以通过设置参数：outfmt 来调整。
  2. 创建容器，用远兄创建的nr数据库（1.1T），挂在目录下去
docker run -it --name nr_mark -v /data3/liruiyuan/blast/home/work:/home/server ubuntu_mark5
  1. 库已经建好，所以接下来直接比对应该就可以了
    1. 这里指的是用dimond建立的nr数据库，已经建好了就用这个来
    2. 除了dimond比对，也可以用blastp来，不过这个效率没有前者的高
nohup diamond blastp \
    --query female_hap2_rename_protein.fa \
    --db /home/server/nr2/nr.dmnd \
    --out female_hap2_blastp.outfmt0 \
    --outfmt 0 \
    --evalue 1e-5 \
    --max-target-seqs 1 \
    --threads 10 \
    --salltitles &
  1. --salltitles：在输出表格里保留 subject 的完整注释信息（不仅仅是 accession 号）。
  2. --max-target-seqs（1）：保留一个最佳比对结果
  3. 使用他定义的脚本转换成下面格式
    1. 脚本只能用perl版本，perl需要安装bioperl之类之类的（比较麻烦）；作者也有写过python版本的，但那个没法用，没什么报错信息，输出确实空的。
      1. conda install -c bioconda perl-bioperl
      2. 脚本在下方附件
      3. 最终的输出格式是表格式，和eggnog-mapper的输出差不多
perl blast_parser.pl female_hap1_blastp.outfmt0 female_hap1.txt

附件

blast_parser.pl

#!/usr/bin/perl -w
 
use strict;
use Bio::SearchIO;
 
if(@ARGV != 2)  {
    print "perl blast_parser.pl INPUT  OUTPUT \n";
        exit;
}
 
my($input, $output) = @ARGV;
 
open(OUTPUT, ">$output") || die "could not open.\n";
print OUTPUT "Query_name\tQuery_description\tQuery_length\tQuery_start\tQuery_end\tHit_name\tHit_description\tHit_length\tHit_start\tHit_end\tAln_length\tIdentity\tEvalue\n";
my $searchio = Bio::SearchIO->new(-format=>"blast", -file => "$input");
 
while(my $result = $searchio->next_result)  {
    my $query_name = $result->query_name;
    my $query_description = $result->query_description;
    my $query_length = $result->query_length;
    if($result->num_hits) {
    my $hit = $result->next_hit;
        my $hit_name = $hit->name();
        my $hit_desc = $hit->description();
        my $hit_length = $hit->length;
        my $hsp = $hit->next_hsp;
            my $evalue = $hsp->evalue;
            my $strand = $hsp->strand('hit');
            my $query_start = $hsp->start('query');
            my $query_end = $hsp->end('query');
            my $hit_start = $hsp->start('hit');
            my $hit_end = $hsp->end('hit');
            my $aln_length = $hsp->length('total');
            my $identity = $hsp->percent_identity;
            print OUTPUT $query_name,"\t",$query_description,"\t",$query_length,"\t",$query_start,"\t",$query_end,"\t",$hit_name,"\t",$hit_desc,"\t",$hit_length,"\t",$hit_start,"\t",$hit_end,"\t",$aln_length,"\t",$identity,"\t",$evalue,"\n";            
    }
}
close OUTPUT;
exit;