基因功能注释

2020-05-11 本文已影响0人斩毛毛

根据已有的蛋白库，对从基因组上提取到的蛋白序列进行比对，从而获得相应的信息。

常用的数据库：

Nr：NCBI官方非冗余蛋白数据库，包括PDB, Swiss-Prot, PIR, PRF; 如果要用DNA序列，就是nt库
Pfam: 蛋白结构域注释的分类系统
Swiss-Prot: 高质量的蛋白数据库，蛋白序列得到实验的验证
KEGG: 代谢通路注释数据库.
GO: 基因本体论注释数据库

$\color{red}{后续分析，你的蛋白序列中不能有氨基酸以外的字符，比如*}$

可以通过以下集中方式进行注释

blastp

Nr的NCBI收集的最全的蛋白序列数据库，但是无论是用NCBI的BLAST还是用速度比较快DIAMOND对nr进行搜索，其实都没有利用好物种本身的信息。因此在RefSeq上下载对应物种的蛋白序列, 用BLASTP进行注释即可。

# download
wget -4  -nd -np -r 1 -A *.faa.gz ftp://ftp.ncbi.nlm.nih.gov/refseq/release/plant/
mkdir -p ~/db/RefSeq
zcat *.gz > ~/db/RefSeq/plant.protein.faa
# build index
~/opt/biosoft/ncbi-blast-2.7.1+/bin/makeblastdb -in plant.protein.faa -dbtype prot -parse_seqids -title RefSeq_plant -out plant
# search
~/opt/biosoft/ncbi-blast-2.7.1+/bin/blastp -query protein.fa -out RefSeq_plant_blastp.xml -db ~/db/RefSeq/uniprot_sprot.fasta -evalue 1e-5 -outfmt 5 -num_threads 50 &

Swiss-Prot里收集了目前可信度最高的蛋白序列，一共有55w条记录，数据量比较小，

# download
wget -4 -q ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz
gzip -d uniprot_sprot.fasta.gz
## 或者在这里下载植物数据
https://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/taxonomic_divisions/
# builid index
~/opt/biosoft/ncbi-blast-2.7.1+/bin/makeblastdb -in uniprot_sprot.fasta -dbtype prot -title swiss_prot -parse_seqids
# search
~/opt/biosoft/ncbi-blast-2.7.1+/bin/blastp -query protein.fa -out swiss_prot.xml -db ~/db/swiss_prot/uniprot_sprot.fasta -evalue 1e-5 -outfmt 5 -num_threads 50 &

InterProScan

下面介绍的工具是InterProScan, 从它的9G的体量就可以感受它的强大之处，一次运行同时实现多个信息注释。

InterPro注释
Pfam数据库注释(可以通过hmmscan搜索pfam数据库完成)
GO注释(可以基于NR和Pfam等数据库，然后BLAST2GO完成,)
Reactome通路注释，不同于KEGG
一般使用Pfam，PRINTS，PROSITE, ProDom, SMART就可以了

安装

若安装该软件需要以下需求：

64-bit Linux
Perl 5 (default on most Linux distributions)
Python 3 (InterProScan 5.30-69.0 onwards)
Java JDK/JRE version 11 (InterProScan 5.37-76.0 onwards)
Environment variables set
$JAVA_HOME should point to the location of the JVM
$JAVA_HOME/bin should be added to the$ PATH

$\color{red}{具体要求请查看说明文档}$ , 这里

下载并安装interProScan

##下载
wget ftp://ftp.ebi.ac.uk/pub/software/unix/iprscan/5/5.44-79.0/interproscan-5.44-79.0-64-bit.tar.gz
wget ftp://ftp.ebi.ac.uk/pub/software/unix/iprscan/5/5.44-79.0/interproscan-5.44-79.0-64-bit.tar.gz.md5
##文件较大，下载md5并进行检测，若ok则没毛病
md5sum -c interproscan-5.44-79.0-64-bit.tar.gz.md5
##安装
tar -pxzf interproscan-5.44-79.0-64-bit.tar.gz

# where:
#     p = preserve the file permissions
#     x = extract files from an archive
#     z = filter the archive through gzip
#     f = use archive file

安装Panther Models

cd ~/interproscan-5.44-79.0/data
wget ftp://ftp.ebi.ac.uk/pub/software/unix/iprscan/5/data/panther-data-14.1.tar.gz
wget ftp://ftp.ebi.ac.uk/pub/software/unix/iprscan/5/data/panther-data-14.1.tar.gz.md5
#文件较大，下载md5并进行检测，若ok则没毛病
md5sum -c panther-data-14.1.tar.gz.md5
#安装
tar -pxvzf panther-data-14.1.tar.gz

使用Pre-calculated Match Lookup（可选）
电脑能连接 http://www.ebi.ac.uk
即可。
也可以下载本地版；
或者不想要直接在interproscan.properties中注释掉即可
'#precalculated.match.lookup.service.url=https://www.ebi.ac.uk/interpro/match-lookup'

InterProscan的使用

测试

./interproscan.sh -i test_proteins.fasta -f tsv
./interproscan.sh -i test_proteins.fasta -f tsv -dp

若运行上述命令，得到两个文件，则表明可行。

./interproscan-5.29-68.0/interproscan.sh -appl Pfam  \
  -f TSV -i sample.fa -cpu 50 -b sample \
              -goterms -iprlookup -pa

## 参数
-appl告诉软件要执行哪些数据分析，勾选的越多，分析速度越慢，Pfam就行
-I: fasta 序列
-iprlookup:  Option that provides mappings from matched member database: signatures to the InterPro entries that they are integrated into
-goterms:  Option that provides mappings to the Gene Ontology (GO).
-b: output-file-base file_name 
-f: 输出格式
-pa: pathways

输入格式

TSV的格式详解

P51587  14086411a2cdf1c4cba63020e1622579    3418    Pfam    PF09103 BRCA2, oligonucleotide/oligosaccharide-binding, domain 1    2670    2799    7.9E-43 T   15-03-2013
P51587  14086411a2cdf1c4cba63020e1622579    3418    ProSiteProfiles PS50138 BRCA2 repeat profile.   1002    1036    0.0 T   18-03-2013  IPR002093   BRCA2 repeat    GO:0005515|GO:0006302

每一列说明：

蛋白名字
序列MD5 disest
序列长度
所用的库（Pfam/RPINTS et.）
编号
描述
起始位置
终止位置
e-value得分
匹配的状态T: true
日期
interPro 注释编号
interPro 注释描述
GO注释
Pathways 注释

KEGG

可以选择在在线工具KASS (KEGG Automatic Annotation Server)进行KEGG注释。

第一步：选择双端的方法
第二步：设置参数，并上传数据

image.png

其中 $\color{red}{BLAST}$ 可以用于氨基酸和核酸比对，若是核酸，勾选相应选框即可，提交还需在邮箱中确定提交，而后耐心等待结果。

结果分析

image.png

注释文件下载

我看网上大多数关于KAAS的教程，基本都只介绍到上面两步，从上面的界面中下载KO注释以及代谢通路图。

但是这样就有两个小问题，一是KO注释只有ID，无详细名称，还要后续自己去库中对应，麻烦；二是手动一个个下载代谢通路图，也是很费劲，但不下载在线看的话，过一段时间服务器中就没了。

事实上，无需在这两个链接中下载数据，我们有更好的方法。KAAS注释后，会针对咱们这个注释结果，给你提供一个KO注释列表。依次在结果链接的主界面点击“html”>“[BRITE hierarchies]”>“KEGG Orthology (KO)”，可下载该注释表。

默认保存名称“q00001.keg”，这个我也放在的网盘附件中，见“result/q00001.keg”，这个“q00001.keg”就是常规的KEGG通路列表的样式了。随便找个文本编辑器（如Notepad++）打开它。

本次KEGG（KAAS）注释，所有能注释到KO的基因以及其参与的代谢通路，都能在列表中找到。这时候你就可以根据所上传基因（或蛋白）序列ID，或者特定的代谢通路名称，去搜索你想要得知的信息了。当然了，并非所有基因都能匹配到KO（主要是这些基因的功能目前尚位完全知晓），所以未注释到的基因，就不在这个列表中了。

对于这个KO注释列表，若觉得这样看起来还是不方便，我们可以再给它处理一下，比方说排个版（整理成表格样式便于操作），或者统计各代谢通路中所含基因数量，等等。

哦对了，我好像还没说怎么把代谢通路图搞下来是吧。其实KEGG的代谢通路图链接，有种固定格式，我们完全可以仿照它的链接样式，在本地生成链接，据此映射到它的具体通路中，并将所涉及的基因以特定颜色展示出来。

此外，在上文展示的KAAS在线注释中，我们推荐上传编码蛋白的氨基酸序列去作功能注释，并不推荐使用基因的核酸序列。因此在这种情况下，序列比对方式是氨基酸序列-氨基酸序列的比对，同时最终得到的KO注释列表中，所谓的“基因id”其实是“编码蛋白id”。我们这里用的细菌，一个基因只对应一条编码蛋白，倒是还好说。若是真核生物，由于可变剪切等的存在，一个基因会对应多个编码蛋白。因此最好再根据基因id和蛋白id的对应关系（可从基因组注释gff文件中获取），将对应的基因id也和KO注释结果联系起来。

这些操作，需要有编程基础才方便完成。其实也不是很难的说，稍微懂一点就可以实现了。网盘链接中提供了示例脚本，这些都是按个人的习惯随便写的，想着给刚入行的同学们作个参考示例，提供下帮助。

示例：KO列表样式转换，对应基因id
网盘附件python脚本“script/kegg_trans.py”，可以实现KO列表的样式转换以及根据蛋白id进一步对应基因id。

#请使用python3.6或以上版本运行，例如
#python3 kegg_trans.py -h
python3 kegg_trans.py q00001.keg gene.list.txt gene_anno.txt pathway_anno.txt

#q00001.keg，KO注释列表
#gene.list.txt，基因列表，共两列，左侧为基因id，右侧为其对应的蛋白id
#gene_anno.txt，输出文件名称，该文件包中含蛋白id、基因id、KO注释等的对应关系
#pathway_anno.txt，输出文件名称，该文件中给出了注释到的通路图链接等

对于“gene.list.txt”，支持一个基因id对应多个蛋白id的情形（考虑到真核生物一个基因通常有多种编码产物），也支持只提供部分基因或蛋白id的情形（考虑到有时我们只想关注其中的部分特定基因）。

上述转化后得到了两个文件“gene_anno.txt”和“pathway_anno.txt”，可分别见网盘附件“result/gene_anno.txt”和“result/pathway_anno.txt”。这两个文件中的内容观察起来会更方便。

例如打开“gene_anno.txt”，对每个基因归了类，需注意一个基因能够对应注释到多种KO：

再看“pathway_anno.txt”，其中生成了可直达通路图的链接（和上文中注释结果的链接不同，这个链接无时间限制，可即点即看；它的原理是通过这个链接接入后台，每次连入时，服务器那边会根据链接结构当场绘图作为临时展示）：

示例：代谢通路图抓取
如果你想将pathway图片批量下载下来，这儿也提供了一个爬虫参考，见网盘附件“script/pathway_download.py”，可通过第一步转化后的文件中的链接获取通路图。

#python3 pathway_download.py -h
python3 pathway_download.py pathway_anno.txt ./map_download
 
#pathway_anno.txt，上一步得到的转化后的文件，我们将根据其中的通路图链接抓取图片
#map_download，图片存放路径

示例：统计某分类所含基因数量
例如我们可在第二层级作个统计，统计各分类水平下所包含的基因数量，并绘制柱状图，以简要展示本次注释的基本情况。这个在测序公司的结题报告中经常看到，我们仿照他们的样式来一个，见网盘附件“script/pathway_stat.r”和“result/kegg_stat”。

#Rscript pathway_stat.r -h
Rscript pathway_stat.r gene_anno.txt ./kegg_stat
 
#gene_anno.txt，第一步转化后的基因注释文件，该文件包中含蛋白id、基因id、KO注释等的对应关系
#kegg_stat，输出结果的路径

同上文，我在这儿统计的是各分类下对应的注释到的“基因数量”（在第一步转化过程中，通过基因-蛋白id列表，对应到了二者的关系），而不是“编码蛋白数量”。

关于KEGG结果部分来源* KEGG自动注释工具KAAS简介及KO注释表的简单处理，该文提到的网盘地址为：https://pan.baidu.com/s/1gOA08CEMFe7jj9AbYUYhbw

非编码RNA注释

a. rRNA预测
--软件：mammer
--用法：

mammer –S euk –gff rRNA.predict.gff –m Lsu,ssu,tsu –multi < $genome

b. tRNA预测
--软件：tRNAscan-SE
--用法：

tRANscan-SE –o tRNA.predict.gff $genome

c. 其他小RNA预测
--软件：cmscan, 利用rfam数据库
--用法:

cescan –o test.rfam –cpu 20  --tblout $rfam.tab --acc /lustre/Work/software/alignment/Rfrn.12.0.CMs/Rfam.cm $genome