基因功能注释

基因功能注释，简单来说，即是根据已有的蛋白库，对从基因组上提取到的蛋白序列进行比对，从而获得相应的信息。

这里整理了拟南芥、Uniprot、Swissprot和Interproscan数据库比对，从而得到基因的信息。

根据数据库中已知编码基因的注释信息（包括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 ： 代谢通路注释数据库。

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

---

1.拟南芥基因比对

这里需要先从TAIR（https://www.arabidopsis.org/download/index-auto.jsp?dir=%2Fdownload_files%2FGenes%2FTAIR10_genome_release）下载拟南芥蛋白序列（pep.fa）及描述文件（description）（如图所示）。安装DIAMOND进行本地blastp。

下载拟南芥描述文件

diamond makedb --in Ara.fasta -d Ara#构建blastp索引
diamond blastp -d Ara.dmnd -q proteins.fasta --evalue 1e-5 > blastp.outfmt6#blastp得到结果

#比对结果中筛选每个query的最佳subject（jcvi可用conda安装）
python -m jcvi.formats.blast best -n 1 blastp.outfmt6
生成blastp.outfmt6.best

2.Uniprot数据库比对

这里以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;
zcat uniprot_sprot_plants.dat.gz uniprot_trembl_plants.dat.gz > uniprot_plants.dat
#转换格式（用python的SeqIO模块转换）
from Bio import SeqIO
count = SeqIO.convert("uniprot_plants.dat","swiss","uniprot_plants.fa","fasta")
print("Converted %i records" % count)

如此便得到了Uniprot数据库，本地blastp即可得到结果（方法同上）。得到结果后，需要进行ID转换，这里使用Uniprot官网转换（https://www.uniprot.org/id-mapping/），输入ID号后即可下载基因名称等相关信息。

3.Swissprot数据库比对

#数据库下载与解压
wget https://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz
gunzip -d uniprot_sprot.fasta.gz

解压完成后，同样需要进行基因ID转换后得到对应的基因名称。

4.Intrproscan数据库比对

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

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

• InterPro注释
• Pfam数据库注释(可以通过hmmscan搜索pfam数据库完成)
• PANTHER数据库注释
• GO注释(可以基于NR和Pfam等数据库，然后BLAST2GO完成)
• Reactome通路注释，不同于KEGG还有很多数据库。

interproscan的安装

软件本身自带了很多数据库，不需要安装，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的一部分，interproscan 5.47-82.0软件自带PANHTER数据库，所以不用单独下载。

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文件，都下下来。推荐用md5sum检查下载是否成功md5sum -c interproscan-5.47-82.0-64-bit.tar.gz.md5。若成功继续下一步。解压文件tar -pxzf interproscan-5.47-82.0-64-bit.tar.gz。

检查

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

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

输入

fasta 格式的蛋白或核酸序列，序列中不能含有-或 * 等非法字符，若出现会报错，这里需要注意。

Note

• 如果本地化运行，可以添加-dp 参数或把interproscan.properties文件里的#precalculated.match.lookup.service.url=https://www.ebi.ac.uk/interpro/match-lookup注释掉，就不进行match.lookup检查，这个应该是更新数据库的检查。
• 如果不注释，在运行时存储interproscan数据库端网络慢时（大部分情况下）会输出一些Connection timed out和ERROR - Lookup version check failed的报错信息，但the analysis will continue to run locally。本地运行仍在继续。

运行

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

参数

• -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值来快速检验这上传的数据是否已经被注释了，如果是已经注释了就直接出结果。

Interproscan结果解读

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

tsv格式文件每一列含义：

蛋白名字
序列MD5 disest
序列长度
所用的库（Pfam/PANTHER/RPINTS etc）
匹配的数据编号（PF00001/PTHR00001 etc)
功能描述
起始位置
终止位置
e-value得分
匹配的状态T: true
日期
interPro 注释编号
interPro 注释描述
GO注释
Pathways 注释
每行是一条蛋白序列注释到的每个数据库的每一个匹配，所以每个基因会有多行。
根据第四列-数据库名称来提取不同数据库的注释结果。这里提取pfam和PANTHER两个数据库的结果。

Interproscan结果整理

• 1.提取interpro数据库注释内容并转化成单基因单行格式
  cat *.iprscan.tsv |cut -f1,12,13|grep "IPR"|awk '{print $1"\t"$2": "$3}'|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\tinterpro' >iprs.interpro.anno
• Note
  注意排序和去重。
  因为我的geneID是tg00001这种格式，只有两个字母，所以sort排序按数字从第三个sort -k 1.3n开始排序，可以根据不同geneID更改。
  其中awk -F"\t" '{a[1]//g" |sort -k 1.3n部分是单基因多行注释转化成单基因单行注释的命令。
  用cut截取1，12，13三列是因为有些行没有12列及之后数据，直接用awk会失效。
  用grep搜索IPR是因为直接用awk匹配第12列不能匹配完全，比如SUPERFAMILY数据库的比对就匹配不了。
• 2.提取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
• 3.提取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

至此，Interiorscan注释完成。

5.ITAK预测TF/TR/PK

iTAK 来自康奈尔大学，FeiLab的一个预测工具。是依赖于数据库的用于从蛋白质或核苷酸序列中识别植物转录因子 (TF)、转录调节因子 (TR) 和蛋白激酶 (PK)，然后将单个 TF、TR 和 PK 分类为不同的基因家族的工具。

• ITAK在线网站http://itak.feilab.net/cgi-bin/itak/index.cgi

使用

进入网站后点击ITAK Online，如图。

ITAK

提交序列，进行比对。

提交序列

下载结果

下载结果

结果解读

• tf_sequence.fasta：所有鉴定的TF/TR序列
• tf_classification.txt：所有TF/TR的分类，tab制表符分割，包含序列的ID和各自的家族。
• tf_alignment.txt：制表符分割的txt文档，包含所有鉴定到的TF/TR比对蛋白结构域数据库。
• pk_sequence.fasta：所有鉴定到的PK蛋白序列。
• Shiu_classification.txt：所有鉴定到的PK蛋白分类。制表符分割的txt文件，包含序列ID和相应的蛋白激酶家族。
• Shiu_alignment.txt：制表符分割的txt文档，包含所有鉴定到的PK比对蛋白结构域数据库。

结果整合

通过不同数据库的比对，得到了注释结果。通过R语言dplyr包中left_join根据基因ID将多份结果整合在一起，基因注释完成。