GO注释---Tbtools

第一步：根据叶绿体基因组的genbank注释文件获得蛋白编码基因序列

提取序列的python脚本

import sys
from Bio import SeqIO

input_file = sys.argv[1]
output_file = sys.argv[2]

protein_coding = {}

for rec in SeqIO.parse(input_file,'gb'):
    for feature in rec.features:
        if feature.type == "CDS":
            gene_Name = feature.qualifiers["gene"][0]
            gene_Sequence = feature.extract(rec.seq)
            protein_coding[gene_Name] = gene_Sequence

with open(output_file,"w") as fw:
        for a,b in protein_coding.items():
            fw.write(">%s\n%s\n"%(a,b))

使用方法

python extract_CDS_from_gb.py input.gb output.fasta`

第二步：使用diamond将叶绿体的蛋白编码基因与swissprot数据库比对，获得TBtools做GO注释需要的.xml格式文件

参考文献：DIAMOND: 超快的蛋白序列比对软件

下载swissprot数据

wget ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz
bgzip uniprot_sprot.fasta.gz

下载diamond

wget http://github.com/bbuchfink/diamond/releases/download/v0.9.25/diamond-linux64.tar.gz
tar xzf diamond-linux64.tar.gz

无需安装，解压出来即可使用

构建数据库

~/mingyan/Bioinformatics_tools/Diamond/diamond makedb --in uniprot_sprot.fasta -db uniprot_sprot

运行完目录下多了一个uniprot_sprot.dmnd文件

比对自己的数据，我的是核苷酸序列，使用blastx

第三步：使用TBtools进行GO注释

需要准备的文件

idmapping.tb.gz 文件比较大

这里推荐一个下载器 https://motrix.app/ 界面非常干净清爽

go-basic.obo

cp_Protein_coding.xml

第一步获得文件

cp_Protein_coding.xml.TBtools.GOAnno.xls

第二步获得文件

Bhagwa_gene2Go.txt.Level3.count.xls

TBtools做GO注释如何具体操作大家可以关注TBtools作者在腾讯课堂开设的一系列视频课程。

这样GO注释就做好了，TBtools也会对应有可视化工具，这里我选择使用R语言的ggplot2进行展示

library(ggplot2)
df<-read.csv("Bhagwa_cp_protein_coding.csv",header=T)
dim(df)
colnames(df)
df1<-df[which(df$Class=="biological_process"),]
df2<-df[which(df$Class=="cellular_component"),]
df3<-df[which(df$Class=="molecular_function"),]
df1<-df1[order(df1$Counts,decreasing = T),]
df2<-df2[order(df2$Counts,decreasing = T),]
df3<-df3[order(df3$Counts,decreasing = T),]
df4<-rbind(df1,df2,df3)
dim(df4)
levels<-as.vector(df4$Description)
df4$Description<-factor(df4$Description,
                        levels = levels)
ggplot(df4,aes(x=Description,y=Counts,fill=Class))+
  geom_bar(stat="identity")+
  theme_bw()+
  theme(axis.text.x = element_text(angle=90,
                                   vjust=0.2,hjust=1),
        legend.title=element_blank(),
        legend.position = c(0.8,0.8),
        panel.grid = element_blank())

图片.png

对结果进行可视化遇到的问题

• 数据框如何根据指定列分组排序，比如我的数据

  X Y
1 A 1
2 A 2
3 B 3
4 B 4
5 C 5
6 C 6

ABC分别从大到小排序，如何实现自己还没有想到比较好的办法。

• ggplot2X轴文本对齐方式采用的是vjust和hjust参数，更改这两个参数

library(ggplot2)
df<-read.csv("Bhagwa_cp_protein_coding.csv",header=T)
dim(df)
colnames(df)
df1<-df[which(df$Class=="biological_process"),]
df2<-df[which(df$Class=="cellular_component"),]
df3<-df[which(df$Class=="molecular_function"),]
df1<-df1[order(df1$Counts,decreasing = T),]
df2<-df2[order(df2$Counts,decreasing = T),]
df3<-df3[order(df3$Counts,decreasing = T),]
df4<-rbind(df1,df2,df3)
dim(df4)
levels<-as.vector(df4$Description)
df4$Description<-factor(df4$Description,
                        levels = levels)
ggplot(df4,aes(x=Description,y=Counts,fill=Class))+
  geom_bar(stat="identity")+
  theme_bw()+
  theme(axis.text.x = element_text(angle=90,
                                   vjust=-0.5,hjust=0.5),
        legend.title=element_blank(),
        legend.position = c(0.8,0.8),
        panel.grid = element_blank())

图片.png