前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >使用TBtools对叶绿体蛋白编码基因进行GO注释

使用TBtools对叶绿体蛋白编码基因进行GO注释

作者头像
用户7010445
发布2020-03-03 14:55:50
5.2K0
发布2020-03-03 14:55:50
举报
文章被收录于专栏:小明的数据分析笔记本
第一步:根据叶绿体基因组的genbank注释文件获得蛋白编码基因序列

提取序列的python脚本

代码语言:javascript
复制
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数据

代码语言:javascript
复制
wget ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz
bgzip uniprot_sprot.fasta.gz

下载diamond

代码语言:javascript
复制
wget http://github.com/bbuchfink/diamond/releases/download/v0.9.25/diamond-linux64.tar.gz
tar xzf diamond-linux64.tar.gz

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

构建数据库

代码语言:javascript
复制
~/mingyan/Bioinformatics_tools/Diamond/diamond makedb --in uniprot_sprot.fasta -db uniprot_sprot

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

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

代码语言:javascript
复制
~/mingyan/Bioinformatics_tools/Diamond/diamond blastx --db uniprot_sprot -q output.fasta -o cp_Protein_coding.xml --outfmt 5
第三步:使用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进行展示

代码语言:javascript
复制
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())

image.png

对结果进行可视化遇到的问题
  • 数据框如何根据指定列分组排序,比如我的数据
代码语言:javascript
复制
  X Y
1 A 1
2 A 2
3 B 3
4 B 4
5 C 5
6 C 6

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

  • ggplot2X轴文本对齐方式采用的是vjust和hjust参数,更改这两个参数
代码语言:javascript
复制
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())
本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2019-09-09,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 小明的数据分析笔记本 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 第一步:根据叶绿体基因组的genbank注释文件获得蛋白编码基因序列
  • 第二步:使用diamond将叶绿体的蛋白编码基因与swissprot数据库比对,获得TBtools做GO注释需要的.xml格式文件
  • 第三步:使用TBtools进行GO注释
  • 对结果进行可视化遇到的问题
相关产品与服务
数据库
云数据库为企业提供了完善的关系型数据库、非关系型数据库、分析型数据库和数据库生态工具。您可以通过产品选择和组合搭建,轻松实现高可靠、高可用性、高性能等数据库需求。云数据库服务也可大幅减少您的运维工作量,更专注于业务发展,让企业一站式享受数据上云及分布式架构的技术红利!
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档