分析记录|细菌基因组功能注释

2017-11-11 本文已影响0人琼脂糖

功能注释

需要faa蛋白序列，prodigal的pep就是蛋白序列

1.BLAST数据库注释

1. blast

需要blast的分析项目：nr,trembl,pfam,cog,cazy,swissprot,vfdba,vfdbb,ardb,phi
大小和耗时
nr：67G 7天
trembl: 35G 5天
pfam去重复：3.8G 1天
cog：792M
cazy：442M
swissprot: 256M
vfdb-a,vfdb-b: 10M
ardb：3.6M
phi: 2.5M

2.bestHit：

blast直接参数指定只保留最优结果 -maxtargetseqs 或者按score排序geneid去重，只保留最优score
文件夹bestHit，文件作为第3步的输入文件。

3. 添加注释信息

3.1 统一处理（7/10)

nr,trembl,cazy,swissprot,vfdba,vfdbb,ardb
fasta文件格式见下

swissprot

从fasta文件取所有>列，得到anno文件

image.png
python db_anno_extract.py。输入bestHits和anno文件，得到根据ID匹配的后面所有内容。得到anno.out文件

image.png

3.2 需要特殊处理的(3/10)

phi的fasta文件格式不同，anno的提取
cog的大类和term信息需要另外获取。
pfam的注释需要另外获取。

A.phi

blast结果
注释处理：包含全部信息的第2列是井号分隔的。第2列仅输出其中的1，2；添加第13列输出空格连接的3，4，5，6列

B.pfam

比对的输出结果

image.png

ID并不是pfam ID，而是uniprot ID. 而且只需下划线之前的部分
根据结果得到 gene-uniprot ID（第1列，第二列部分）
利用swiss-ID得到Pfam ID （Pfam-A.regions.uniprot.tsv）

image.png
从pdb_pfam_mapping.txt 提取pfam（只要点号前）和description，去重复。

image.png
pfam ID 得到term
pfam-besthit 加上term

C.COG

1.3个database

cognames2003-2014.tab

cog2003-2014.csv

fun2003-2014.tab

blast结果有注释，但需要更大的类的注释和分组可以统计和作图。

image.png

3.根据gi编号查询cog2003-2014.csv.tab得到COG编号；COG编号查询cognames2003-2014.tab 得到大的类别分类和注释

D.CAzy

数据库(link)没有提供下载。只有分物种的online信息
网上很多parsing和抓去的工具 link link
一个online数据库 link，得到fasta文件进行blast。
结果中没有注释，id是geneID+分类

AFY26876.1|AA0

各个分类CBM，CE，GT等等需要计数统计，填写报告。

数据库中有CAZyDB-ec-info.txt.07-20-2017

genebankID-cazyClass-st-term

没有进一步注释，注释信息不全（缺少AA分类）

2.GO注释

gene ID 和GO的对应不是一对一，而是一对多。

通过swissprot-bestHit得到gene和unirpotID对应关系
goa_uniprot_all.gaf 第二列（uniprotID）和第五列（GO编号）对应去重复

image.png
按照swissprotID 得到GO分类
awk 'BEGIN{FS=OFS="\t"}FNR==NR{a[$2]=$1;next}($2 in a){print a[$2],$0}' gene_swiss.id gaf > go.temp
ID可能对应多个GO分类，因此gaf作为逐行读取。只要我们需要的gene ID有匹配，就输出。
按照GO分类和go.obo得到term和分类。
GO计数用于作图。
将同一gene ID的GO整合。

image.png
存在问题 GO的backgroud 如何计算；作图只能选择count

3.KEGG

在线分析
(1) 进入kaas网上在线服务网站：http://www.genome.jp/kaas-bin/kaas_main
(2) 提交faa蛋白序列，点击compute，预留邮箱，进入邮箱点击相关链接提交任务
(3) 下载ko注释信息：http://rest.kegg.jp/list/ko
(4) 给注释到的基因进行ko注释
将KO编号和term连接。

4.软件分析注释

1.分泌蛋白signalP

signalP 同 rRNAmmer
都要vim signalIP，找到安装目录，修改脚本的目录位置。

输出文件为gff和summary-out
根据gff的YES计数，可得到分泌蛋白个数。
summary-out 计数YES，但需要除2.每个gene出现两次。

2.分泌蛋白T3SS EffectiveT3

link

###Usage: 
java -jar TTSS_GUI.jar -f <inputfile> -m <modulename> -t <restriction> -o <outputfile> [-q] 
##or
java -jar TTSS_GUI.jar -f <inputfile> -c <configurationfile> [-q]

java -jar TTSS_GUI-latest.jar -f test.fasta -m TTSS_ANIMAL-1.0.1.jar -t selective -o test.out -q
java -jar TTSS_GUI-latest.jar -f test.fasta -m TTSS_STD-latest.jar -t cutoff=0.995 -o test.out -q

关于module

Download the Effective classification modules:
TTSS_STD-2.0.1.jar     (based on all training data; 09/2015)
TTSS_STD-1.0.1.jar     (based on all training data; 08/2009)
TTSS_ANIMAL-1.0.1.jar  (based on training data from human/animal-associated bacteria; 08/2009)
TTSS_PLANT-1.0.1.jar   (based on training data from plant-associated bacteria; 08/2009)

3.运行注意
要将要用的module即-m参数指定的文件，复制到工作目录的module文件夹下

3. 次级代谢基因簇 antismash

本地版本

安装版本报错

web版本

分为细菌和真菌 link
选择了clusterblast 和subclusterblast ，其他参数默认
运行完，下载全部result.zip

html的表格直接复制下来，作为最终结果。
统计基因个数clusterblastoutput.txt
cat clusterblastoutput.txt |sort -k1,1 -u|awk '{print $1}' > temp
wc -l temp
554
cat temp |grep "|c1|"
第1簇有9个基因
html可以直观的看出
html结果
从起始和终止位置：71278-81315 ，结合gff也可以数出来
image.png

其他折腾

Ardb

1.网站提供了ARdb注释工具（包含蛋白database）link
2.按照说明文档
进入blastdb中，建立index。
当ARdb文件目录下，修改输入文件路径
在ARdb文件目录下，运行perl

makeblastdb -in resisGenes.pfasta -dbtype prot -out T
vim genomeList.tab
perl ardbAnno.pl

3.运行输出

######################################
BLASTing /home/liyb/data/assembly.faa against ARDB
######################################
######################################
Annotating BLAST file test.blastp
######################################
Create Excel Table output.xls
4.没有结果。。

方法二
还是用blast