基因功能的注释依赖于上一步的基因结构预测，根据预测结果从基因组上提取翻译后的蛋白序列和主流的数据库进行比对，完成功能注释。(有些也支持用DNA序列做功能注释的，但比较少。)

比较推荐的方案是：eggNOG-mapper网页注释和PANNZER2网页注释，如果有余力再加上interpro的本地注释。

功能注释的一般流程：1.选择数据库并下载，2.构建blastp索引，3.用blastp比对蛋白序列到数据库，4.结果整理。

有些软件集成了数据库和功能注释的多个流程，甚至提供服务器用于网页版注释，例如eggNOG-mapper和interproscan分别集成了eggNOG和interpro数据库，PANNZER2使用最新的UniProt数据库实现网页注释。

根据数据库中已知编码基因的注释信息（包括motif、domain），基于同源比对，对基因中的模序和结构域、新基因编码的蛋白质功能、所参与的信号传导通路和代谢途径等的预测。基因组注释内容还可涉及蛋白激酶、病原与宿主互作、致病毒力因子预测、抗性基因等。常用数据库：

Nr：NCBI官方非冗余蛋白数据库，包括PDB, Swiss-Prot, PIR, PRF; 如果要用DNA序列，就是nt库。

UniProt：分为两部分。Swiss-Prot是UniProt数据库中经过人工校正的高质量蛋白数据库，蛋白序列得到实验的验证。TrEBL是UniProt数据库中未经过人工校正、机器自动注释的蛋白数据库。

eggNOG：EMBL创建，是对NCBI的COG数据库的拓展，提供了不同分类水平蛋白的直系同源分组（Orthologous Groups，OG）。

InterPro：EBI开发的一个整合的蛋白家族功能注释数据库，包括Gene3D、CDD、Pfam等10几个数据库。

Pfam: 蛋白结构域注释的分类系统。

GO: 基因本体论注释数据库。

KEGG: 代谢通路注释数据库。

数据库的下载非常耗时，不管是eggnog还是interproscan，linux下几十KB/s，window下几百KB/s(需科学上网)，耗时长度以天计算，可以尝试先在window系统下载，再转移到linux下。

nr数据库；nt数据库；

UniProt数据库包括人工校正的高质量数据库Swiss-Prot和软件自动注释的TrEMBL数据库。

Swiss-Prot数据库：高质量，手工注释，非冗余的蛋白质序列数据库。非常可靠。

TrEMBL数据库：自动注释的蛋白质序列数据库。

UniProt将完整数据库分类拆分成几个子库。包括古菌，真菌，植物，人，哺乳动物，脊椎动物，无脊椎动物，啮齿动物等。

[GO数据库](http://geneontology.org/docs/download-ontology/）

KEGG(Kyoto Encyclopedia of Genes and Genomes)京都基因和基因组百科全书，是研究Pathway代谢通路最主要的数据库，整合了基因组、化学、系统及疾病和健康的信息。

COG数据库:NCBI开发的用于同源蛋白注释的数据库。fun2003-2014.tab是COG的分类信息，将COG的功能分为26个类别，每个类别用一个字母表示。cognames2003-2014.tab是COG的详细信息，包括编号，对应的分类，功能描述等。cog2003-2014.csv是蛋白和COG对应关系。prot2003-2014.fa.gz是fa格式的蛋白序列和注释信息。

KOG（Clusters of orthologous groups for eukaryotic complete genomes）真核生物蛋白相邻类的聚簇。构成每个KOG的蛋白都是被假定来自一个祖先蛋白，要么是orthologs，要么是paralogs。orthologs是指来自不同物种的由垂直家系（物种形成）进化而来的蛋白，并且保留与原始蛋白相同的功能；paralogs是指在一定物种中来源于基因复制的蛋白，可能会进化成新的与原来有关的功能。KOG数据库

eggNOG（evolutionary genealogy of genes：non-supervised orthologous groups）由EMBL创建，是对NCBI的COG数据库的拓展，提供了不同分类水平蛋白的直系同源分组（Orthologous Groups，OG），包括真核生物，原核生物，病毒的信息。拓展了COG的分类，采用无监督聚类算法在全基因组范围内推导基因功能，更适合谱系特征基因的分析。

eggNOG数据库中，e5.proteomes.faa是所有的蛋白组序列，e5.taxid_info.tsv是taxid对应的物种名及完整的谱系信息，e5.og_annotations.tsv是所有NOG group信息（第一列为Taxid，第二列为NOG groups，第三列为COG归属，第四列为Function）。per_tax_level下不同taxonomy level的members.tsv文件分别是相应level的蛋白序列id和NOG group的对应关系（第一列taxid，第二列NOG groups，第三列为该NOG group所包含的蛋白序列数目，第四列为该NOG group所包含的物种数目，第五列是该NOG group所包含的蛋白序列id，第六列是该NOG group所包含的物种taxid。

pfam数据库是蛋白质家族数据库，根据多序列比对结果和隐马尔科夫模型，将蛋白分成不同家族。

string数据库。搜寻已知蛋白质和预测蛋白质之间的相关关系的系统，包括蛋白质之间的物理作用和间接的功能相关性。基于染色体临近，系统进化谱，基因融合和基因芯片数据等计算基因或蛋白的共表达。包括动物转录因子数据AnimalTFDB3.0，植物转录因子数据PlnTFDB，真菌转录因子数据库Fungal TFDB。

植物转录因子数据库收录了大部分植物模式物种共20个种的84个转录因子家族，包含28193 protein models，26184 distinct protein sequences，支持在线blast和本地blast。分析工具：iTAK，软件内置了PlantTFDB的数据库，可直接用于预测植物转录因子。结果在classification.txt文件中。

notes: 注释分析中要保证蛋白序列中没有代表氨基酸字符以外的字符，比如说有些软件会把最后一个终止密码子翻译成”.”或者”*”，需要删除。

这里以Uniprot为例，下载其他数据库后可用同样方法进行手动注释。

可以uniprot_sprot和uniprot_trembl的植物数据库合并在一起。wget ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/taxonomic_divisions/uniprot_sprot_plants.dat.gz #;wget ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/taxonomic_divisions/uniprot_trembl_plants.dat.gz # 11.5G;zcat uniprot_sprot_plants.dat.gz uniprot_trembl_plants.dat.gz > uniprot_plants.dat;

转换格式：用python的SeqIO模块转换

makeblastdb -in uniprot_plants.fasta -dbtype prot -title uniprot_plants -parse_seqids或者更快的diamond makedb --in uniprot_plants.fasta -d uniprot_plants

blastp -query protein.fa -out uniprot_plants.xml -db uniprot_plants.fasta -evalue 1e-5 -outfmt 6 -num_threads 24 &或者更快的diamond blastp -d uniprot_plants.dmnd -q proteins.fasta --evalue 1e-5 > blastp.outfmt6

jcvi可用conda安装python -m jcvi.formats.blast best -n 1 blastp.outfmt6生成blastp.outfmt6.best

python add_annotation_from_dat.py blastp.outfmt6.best uniprot_plants.dat生成swiss_annotation.tsv文件；包含以下几列。

参考脚本add_annotation_from_dat.py

blast2go的功能是将blast注释结果转换成GO注释，但使用较为复杂，下载数据库耗时长，不推荐使用。blast2go教程；blast2go教程下载数据库wget ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/gene_info.gz #650Mwget ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/gene2accession.gz #2.3Gwget ftp://ftp.pir.georgetown.edu/databases/idmapping/idmapping.tb.gz #10Gwget http://release.geneontology.org/2017-01-01/mysql_dumps/go_monthly-assocdb-data.gz #6.3G

eggNOG(evolutionary genealogy of genes: Non-supervised Orthologous Groups)。eggNOG数据库是数据库，eggnog-mapper是使用eggnog数据库作为参考的注释软件或网页，eggnog-mapper网页版。

下载数据库本地注释的方法：在eggnog网站下载对应使用的数据库。http://eggnog5.embl.de/download/emapperdb-5.0.0/eggnog.db.gzhttp://eggnogdb.embl.de/download/emapperdb-5.0.0/eggnog_proteins.dmnd.gz注释命令：python2 /home/leon/software/3.Function/eggnog-mapper.orginal/emapper.py -i Tov.final.all.maker.transcripts.fasta -o run_egg –override -m diamond –cpu 80

注释分为三步：

推荐使用eggnog-mapper网页版注释，避免下载数据库的繁琐。限制10万条蛋白序列，数量超过可以切分后分开注释。

只需要上传氨基酸序列，点击提交；然后在收到的邮件里点击链接，开始任务即可。

eggnog-mapper会生成五个文件:

其中[project_name].emapper.annotations.tsv每一列对应的记录如下：

interpro是数据库，interproscan是可以使用interpro的注释软件。如果序列少，可以使用interproscan网页版注释，避免下载数据库的繁琐。限制一次注释一条，长度小于40000个蛋白序列，超过可以切分后分开注释。

interproscan可实现多个数据库同时注释，包括：

如果注释的文件比较大或者比较多，可以下载本地版注释，下载过程非常慢。根据你的linux版本和发布日期来选择最适版本，软件很大最近版大约9.1G。建议下载对应的md5文件，用 md5sum -c *.md5来检查下载的是否完全。

interproscan软件地址：ftp://ftp.ebi.ac.uk/pub/software/unix/iprscan/5/, 选择最新版本，比如5.47-82.0，则wget ftp://ftp.ebi.ac.uk/pub/software/unix/iprscan/5/5.47-82.0/*，版本号目录下有两个文件，一个软件的tar.gz压缩文件，一个对应的md5文件，都下下来。

软件本身自带了很多数据库，不需要安装，5.47-82.0版本带有CDD-3.17,Coils-2.2.1,Gene3D-4.2.0,Hamap-2020_01,MobiDBLite-2.0,PANTHER-15.0,Pfam-33.1,PIRSF-3.10,PRINTS-42.0,ProSitePatterns-2019_11,ProSiteProfiles-2019_11,SFLD-4,SMART-7.1,SUPERFAMILY-1.75,TIGRFAM-15.0数据库。

PANHTER(Protein ANalysis THrough Evolutionary Relationships) 数据库是Gene Ontology Phylogenetic Annotation Project的一部分。

软件的运行依赖java11和python3，java -version和python --version检查版本是否正确。

解压缩后运行 interproscan.sh 测试是否成功安装，弹出help界面就是成功安装，用测试文件运行下，interproscan.sh -i test_proteins.fasta，没报错则软件可正常运行。

常用运行参数interproscan -i pep.fa -b out.iprscan -goterms -iprlookup -pa -dp -cpu 24

interproscan的参数：

-i,–input 输入文件，一般要为fasta格式，注意不要带有除氨基酸符号的其他特殊符号（比如代表终止密码子的*），gene/氨基酸的名称不能为空。

-b,–output-file-base 指定输出文件的路径和文件名，默认是输入文件的路径

-appl,–applications 指定使用Interpro中哪些数据库，默认使用全部数据库

-f,–formats 用于指定输出文件的后缀，蛋白序列默认输出TSV, XML and GFF3

-t,–seqtype 输入文件的序列类型，p为protein，n为dna/rna，默认是p

-cpu 指定使用线程数

-goterms 打开GO注释，依赖iprlookup参数，switch on lookup of corresponding Gene Ontology annotation (IMPLIES -iprlookup option)

-iprlookup Also include lookup of corresponding InterPro annotation in the TSV and GFF3 output formats.

-pa Optional, switch on lookup of corresponding Pathway annotation (IMPLIES -iprlookup option)

-dp 关闭precalculated match lookup service，默认的是开启。根据md5值来快速检验这上传的数据是否已经被注释了，如果是已经注释了就直接出结果。

生成的out.tsv和out.gff3等是不同格式的注释文件，这里使用tsv格式。

tsv格式文件每一列含义：

每行是一条蛋白序列注释到的每个数据库的每一个匹配，所以每个基因会有多行。根据第四列-数据库名称来提取不同数据库的注释结果。这里提取pfam和PANTHER两个数据库的结果。

notes：

提取Pfam数据库注释结果并转化成单基因单行格式cat *.iprscan.tsv |awk '$4 == "Pfam" {print $1"\t"$5": "$6}' |sort -k 1.3n |uniq|awk -F"\t" '{a[$1]=a[$1]$2"; "}END{for( i in a){print i"\t"a[i]}}'|sed "s/; $//g" |sort -k 1.3n |sed '1i\gene\tinterproscan.pfam'>iprs.pfam.anno

提取PANTHER数据库注释结果并转化成单基因单行格式cat *.iprscan.tsv |awk '$4 == "PANTHER" {print $1"\t"$5": "$6}'|sort -k 1.3n|uniq |awk -F"\t" '{a[$1]=a[$1]$2"; "}END{for( i in a){print i"\t"a[i]}}'|sed "s/; $//g" |sort -k 1.3n |sed '1i\gene\tinterproscan.panther'>iprs.panther.anno

网站注释：http://ekhidna2.biocenter.helsinki.fi/sanspanz/

注释速度：7万条序列4个小时。

notes:

每个基因有只有一行第二列type为DE的数据，此时第六列的描述是对基因的功能描述，即我们想要的基因功能注释信息。

cat pannzer2/anno.out.txt |awk '$2 == "DE" {print $1"\t"$6}' |sort -k 1.3n |uniq |sed '1i\gene\tpannzer2.uniprot' >pannzer2.uniprot.anno把Uniprot数据库功能注释信息提取出来。

cat pannzer2/GO.out.txt |awk '{print $1"\tGO:"$3,$4}' |awk -F"\t" '{a[$1]=a[$1]$2"; "}END{for( i in a){print i"\t"a[i]}}'|sed "s/; $//g" |sort -k 1.3n |uniq |sed "1s/.*/gene\tpannzer2.go/" >pannzer2.go.anno把Uniprot数据库的GO注释提取出来并转化成单基因单行信息。

网站注释：https://www.plabipd.de/portal/web/guest/mercator4/-/wiki/Mercator4/FrontPage

注释速度：7万条蛋白3min注释完，神速！

两个文件，点击即可下载。

cat mercator4v3/*.annotations.results.txt |awk -F"\t" '$3 ~ "geneIDprefix" {print $3"\t"$2": "$4}' |sed "s/'//g" |sort -k 1.3n |uniq |grep -v "no hits"|awk -F"\t" '{a[$1]=a[$1]$2"; "}END{for( i in a){print i"\t"a[i]}}'|sed "s/; $//g" |sort -k 1.3n |sed "s/not assigned.annotated/-/g" |sed "1i\gene\tmercator">mercator.anno

注意检查是否还有no hits或者not assigned.annotated，有的话都替换掉。

GFAP是windows下的一款基于python开发的软件，用于为植物基因组做功能注释。2022年初发表在bioRxiv上，特地是运行非常快。

提交功能注释文件有两种方案：

功能注释以product属性值的形式提供，添加到结构注释gff文件的第九列。

不同软件的基因功能注释结果都预先进行了整理（参考各个软件后的注释整理部分），整理成两列的格式，首行为表头，第一列为基因ID（gene），第二列为基因功能注释的描述信息（标明注释来源），每个注释软件/数据库的描述信息单独为一列。

目前共有9个文件被整理出来：

合并步骤：

整理注释结果先把所有功能注释转换成单基因单注释的格式(单基因多行注释)，确保数据只有两列(基因ID列和功能注释列)，并删除除了product值以外的内容，比如iprs.panther.anno文件中的pantherid号，整理成pannzer2.uniprot.anno文件的形式，两列数据，第一列基因ID，第二列注释到的产物product值。

合并多个注释结果cat *anno|sort -f|uniq -i|awk ‘$2 != “” {print $0}’ >merge.anno合并多个注释结果，排序sort，去重uniq，和去除空值行等处理；按照基因ID排序，生成单基因多行注释的merge.anno文件

过滤注释最好的办法是在这一步做过滤filter。

sed -i -E -e "s/[,|()[]{}]//g" -e "s/^[-:/]//g" -e "s/[-:/]$//g" -e "s/[^a-zA-Z]+//g" merge.anno

product值的不规范情况：

单注释转多注释用cat merge.anno |awk -F"\t" '{a[$1]=a[$1]$2","}END{for( i in a){print i"\t"a[i]}}' |sed "s/,$/;/g"|sort -gk 1.3n|awk '{print $1"\tproduct="$2}' >product.list命令从merge.anno生成product.list，每个基因单行的格式，还是两列，第一列基因ID，第二列若有多个注释用逗号隔开，用分号结尾。

添加product属性值sample.gff与product.list合并，把结构注释和功能注释文件合并，功能注释以product属性值的形式加入gff文件的第九列内容。

不同软件的基因功能注释结果都预先进行了整理（参考各个软件后的注释整理部分），整理成两列的格式，首行为表头，第一列为基因ID（gene），第二列为基因功能注释的描述信息（标明注释来源），每个注释软件/数据库的描述信息单独为一列。

目前共有9个文件被整理出来：

需要整合所有注释软件/数据库的信息，只需根据第一列进行文件的合并。

实现合并a.anno和b.anno的方法：a.anno和b.anno都只有两列的情况

实现合并ab.anno和c.anno的方法：ab.anno有三列，c.anno有两列的情况

如果多个文件，则两两依次合并。还没想到什么更简单的方法，就这样用join手动合并所有注释软件/数据库的功能注释到一个文件吧。

请大佬写了个合并功能注释的程序，merge_file_key_property.title.cpp，可以一次合并任意多个文件。

./merge_file_key_property.o -a .anno -n 7 -p mc -q gene -o merged 使用示范：-a指定合并的文件类型（比如.anno）；-n指定第一列字段长度（比如sp00001长度为7）；-p指定第一列关键词的前缀（基因名前缀，比如sp）；-q指定第一列关键词的标题（比如gene），作为表头；-o指定合并文件名的前缀。

如果文件多且特征复杂，还是琢磨一下python或者R的merge功能吧。