基因功能注释-SwissProt和Gene Ontology（四

写在前面

前面的流程，我们已经可以

1. 得到原始数据（下载或者自己测序）
2. 质控和比对
3. 定量
4. 得到差异表达分析

差异表达分析，得到了一堆的基因。但是这些基因都是干啥的，我们并不知道。所以这个时候，如果没有功能注释的话，就需要进行功能注释。既然要注释，那么后续还会需要做功能富集分析。所以一般不会只针对差异表达基因注释。而是对物种的所有蛋白或者所有转录本进行功能注释。

很久以前，一般是比对到NCBI的NR注释，现在往往是直接比对到Swissprot。
又写了一个小脚本，课题组的人，可以一个命令完成注释。
（当然，脚本会自动识别输入的是 ，蛋白序列，还是 核酸序列）

使用脚本

需要准备的文件只有一个
genome.cds.fa # 即 所有CDS序列
或者是
genome.pep.fa # 即 所有PEP，即蛋白序列

需要运行的只有一个命令

perl /tools/XiaLabRNAseqPipe/GOanno.pl genome.pep.fa genome.query2go.map

其中
genome.pep.fa 即输入的文件
genome.query2go.map # 即输出的GO注释文件

此外，运行过程中会生成有用的中间文件，包括

• sample.fa.2.swissprot.xml 即 blast输出的XML比对结果文件，可用于可视化
• sample.fa.2.swissprot.xml.TBtools.table 即 XML转换成TBtools Table，可用于逐个查看基因的功能
• sample.fa.2.swissprot.xml.xls 即 原始的GO注释结果

写在后面

使用这个脚本，目前可以完成GO注释。
而或许后面还会做KEGG注释，那个再算。
而有了GO注释，那么就可以使用TBtools进行GO富集分析（无论是命令行，还是界面版）