测序了，然后呢（二） | 基因功能注释

豆豆写于19.1.16 

上次介绍了关于基因预测的知识（为什么要搞全基因组测序（一））
这次是第二部——基因功能注释

背景

得到一个基因集以后，需要知道基因有哪些功能，参与哪些生物过程，只有理解了基因的功能以后，才能联系起来基因型与表型。

继人类基因组计划（HGP）完成以后，找到了2万多个编码蛋白基因，但是基因功能还未知，然后又有了ENCODE（Encyclopedia of DNA Elements）计划，找到了400万个基因开关，平均一个基因有约200个开关，这是人类既有共性又有个性的基础。ENCODE计划使人类基因组不再是个”空壳“
https://www.encodeproject.org/

关于ENCODE计划：https://kknews.cc/science/pq2nx8.html

简而言之，功能注释就是预测蛋白序列的功能，是最基本的分析之一

功能注释应用场景

• 从头拼接的基因组并做了结构注释，知道了哪些地方是外显子，哪些是内含子，接下来就是功能注释，预测每一条基因编码什么蛋白，并且蛋白是什么功能

• 无参转录组需要从头拼接转录本，拼接的转录本功能需要做注释

• 得到了差异表达基因，想做下富集分析，就必须要了解每个基因对应哪个GO分类，也是需要进行功能注释

原理

基因不同于蛋白，不能通过结构来预测功能，只能通过与已知基因功能数据库的比对去推测。一般的数据库包括了两部分内容：一是基因序列（核酸+氨基酸）FASTA格式；二是基因功能信息（可以写到FASTA的ID行中或者单独放在一个文件中）

一般采用氨基酸序列与数据库进行相似性比对，比对结果去数据库中进行过滤

这里看到，基因功能注释主要依赖数据库，如果数据库中没有这个基因，那么就无法注释。更可怕的是，数据库中有错误，就会进行错误注释

比对的结果并不是百分之百完全比对的，那么怎么判断氨基酸序列和数据库的关系呢？比对到多少才能被接受？这里需要考虑比对长度、比对分值、identity值等，过滤掉一部分人为认定不满足同源关系的序列。但是又有一个问题，不同区域的基因会发生不同程度的突变，如果仅设置一个值进行过滤——”一刀切“，这个结果还是有待优化

另外，如果结果提示：Selenocysteine (U) at position ** replaced by X 说明U氨基酸被替代成了X（当然并不是错误，可以忽略），因为在blastp/tblastn的打分矩阵中不存在U和- 这两个字符，替换成任意字符X就可以任意打分
https://www.biostars.org/p/111143/

基本流程

如果手头仅仅有几条蛋白序列想做下功能注释，那么直接甩给uniprot/ncbi在线blast比对就可以了，但是我们这里说的情况是成千上万条基因，肯定不能在线提交，那么怎么办？

要进行大量蛋白序列的功能注释，需要包括：同源注释、功能分类

同源注释

• 基于相似性的注释：就是将要研究的序列与蛋白数据库进行比对，将数据库中比对相似性高的蛋白序列可以作为研究序列的功能，常用的是Nr、Uniprot数据库 ，常用软件是blast和diamond 【其中，blast速度很慢，比对几万条序列可能好几天甚至一周；diamond也是基于blast但速度最快达到blast的两万倍，准确性差不多，因此一般就使用diamond就好】

blast是基于动态规划算法，就是将每个位点都进行比对，比对上就得分，比对失败就罚分。从准确性讲是不错的，但是这个方法对于背后的生物学特性欠缺考虑。因为不是每个氨基酸都是一样重要的，对于某些抗性基因或者转录因子，真正起作用的往往是一些保守结构域

• 基于结构域的注释：Pfam（https://pfam.xfam.org/）数据库中有各种基因家族的保守域模型，可以用HMMER软件将研究序列与数据库中的模型进行比对，如果序列上存在某个结构域，那么推测序列含有该结构域功能；另外Interpro（https://www.ebi.ac.uk/interpro/）是一个综合数据库，使用interproscan软件比对

做完同源注释，就知道了研究的序列和数据库中的哪个蛋白最相似，我们主要利用了nr、uniprot、pfam、interpro这些蛋白数据库，它们又和下游的GO、KEGG、COG等分类数据库有关联，然后就能知道研究的蛋白属于哪个GO分类，哪个Pathway，哪个基因家族，就是功能分类

功能分类

只了解单个基因的功能是不够的，因为基因间是相互作用、协同完成生物功能的，所有需要进行分类，这就是在RNA-seq中得到差异表达基因后做的功能分类（GO）和富集分析（KEGG）过程，看看基因是不是协同完成某一个生物过程，它的原理与功能注释相似，也是利用已有的分类去推测未知的分类

小Tip：功能注释相当于一个过滤筛。GO 注释=》粗筛；KEGG=》细筛，例如：某一个蛋白，GO 只能将它注释到与细胞凋亡有关；而 KEGG 则可以将它注释到细胞凋亡通路中的某一个环节

例如COG数据库（Cluster of Orthologous Groups of proteins， https://www.ncbi.nlm.nih.gov/COG/）是细菌、藻类和真核生物的21个完整基因组的编码蛋白，根据系统进化关系分类得来的。

ftp://ftp.ncbi.nih.gov/pub/COG/COG 数据库还是2003年的，所以做出来的东西，看看就好了

# 下载数据库数据
$wget ftp://ftp.ncbi.nih.gov/pub/COG/COG/myva
$wget ftp://ftp.ncbi.nih.gov/pub/COG/COG/fun.txt
$wget ftp://ftp.ncbi.nih.gov/pub/COG/COG/whog

# 清洗COG数据库(只挑有注释的那些序列)
# https://gist.github.com/Buttonwood/96f9a9ef8159ca111a69
$cog_db_clean.pl -myva myva whog > cog_clean.fa

# blast+比对
$makeblastdb -dbtype prot -in cog_clean.fa
$blastp -query yourdata.fa -db cog_clean.fa -e 1e-4 -out blast.out -outfmt 7 -num_threads 10 -seg no

#整理结果 https://github.com/kodayu/blog_html/blob/master/blast_cog.py
$blast_cog.py blast.out fun.txt whog out

参考：http://yk.cancerbio.info/myblog/cog%E5%88%86%E6%9E%90/

又例如GO数据库 ，其中每个注释都是对基因产物的描述，有特定的分子功能（MF），涉及到特定的生物过程（BP），作用在特定的细胞组分（CC）。它把所有候选的靶基因向GO的各个term进行映射，然后计算映射到每个term的靶基因数量，在整个参考基因背景中利用超几何分布检验，选出候选靶标基因中显著富集的GOterm

再例如KEGG数据库，关于生物化学途径的描述，许多活细胞的功能不能仅仅依赖于单个基因，它将基因组信息与更高一级的功能信息结合；另外它可以将基因组中的许多基因利用细胞内分子互作网络联系起来，通过通路或者复合物来展示更高级的生物学功能

下次介绍整合的流程软件