2023-02-08对候选基因(GWAS\转录组)进行GO注释（

作者: 麦冬花儿 | 来源:发表于2023-02-07 14:31 被阅读0次

2023-02-08对候选基因(GWAS\转录组)进行GO注释（
python从基因组注释文件中提取GO信息
2022-12-14使用bedtools进行gwas基因注释
NGS基础 - 参考基因组和基因注释文件
PheWAS(全表型组关联分析)----GWAS and Phe
eggNOG 对基因组进行GO注释
GO注释结果整理
非模式物种的GO富集分析
基因组重测序与转录组联合分析简介
植物GO注释（在线）

目前网上的教程大多是模式植物，或者是非模式物种的重新构建的教程。
今天分享一个新的方法，对非模式植物研究中获得的基因（例如：GWAS、转录组差异基因等）进行转录因子富集分析气泡图和柱状图可视化。
教程如下：

利用plantTFDB对转录因子进行预测

1、准备差异基因序列

提取差异基因序列
faSomeRecords ../Trinity.gene.fasta DEG_ID_FC1.5FDR0.05.txt leaf_vs_fruit_FC1.5FDR0.05.fasta
faSomeRecords ../Trinity.gene.fasta DEG_ID_FC1.5FDR0.05.txt stem_vs_fruit_FC1.5FDR0.05.fasta
faSomeRecords ../Trinity.gene.fasta DEG_ID_FC1.5FDR0.05.txt stem_vs_leaf_FC1.5FDR0.05.fasta
http://plantregmap.gao-lab.org/index.php

2、打开PlantTFDB(5.0)

数据库链接如下：http://plantregmap.gao-lab.org/index.php
点击导航栏的Prediction按钮，即可打开预测页面。

图片.png

3、TF预测

然后输入你的基因序列，比如差异基因序列，即可开始在线分析：

图片.png

获得TF_and_best1_in_Ath.list文件

图片.png
注意：在这一步点击下载选项不会直接下载文件，而是弹出结果页面。可以用wget命令通过服务器进行下载

图片.png

这个分析的作用就是：判断你的哪些（差异）基因可能是转录因子。

4、go分析预测

将TF_and_best1_in_Ath.list中的第三列粘贴到框内

图片.png

开始submit运行
对结果进行下载

图片.png

5、转录因子富集分析

图片.png

富集分析的原理

网上已有很多人介绍其原理，我们平时最常用的富集分析方法就是Over-Represence Analysis(ORA)，本质上是一次不放回抽样的概率事件。在这里首推CJ大神写的推文：GO富集分析从原理到实践 ~ 零基础掌握。大家可以自行阅读，熟悉一下富集分析的原理。

转录因子的富集分析

当你了解了富集分析背后的原理后，你就能明白：所谓转录因子富集分析，就是以全基因组各个家族的转录因子为背景，通过比较给定的基因中各个家族的转录因子的数量及比例来计算其显著性。为了达到这一目的，我们需要做的无非是三件事：

准备包含待研究物种的基因组内所有转录因子的文件；
准备个记录转录因子基因ID和家族的对应关系的文件；
准备待分析的转录因子的基因ID列表；
当然，我们还需要准备富集分析用到的软件，并根据软件要求将数据调整成恰当的格式。该富集分析需要用的软件包是Y叔开发的clusterProfile包，所以针对该包对数据的要求，上述提高的三个数据文件应遵循如下形式：

TF id和物种内所有转录因子基因ID的对应关系，其中TFid是人为规定的（TF2gene）；

图片.png

TF id和转录因子家族的对应关系（TF2term）；

图片.png
待分析的基因ID，比如某一时期处理A和对照B的所有差异基因ID (genelist)

图片.png

整理文件如上述格式

输入文件：PlantTFDB网站对全基因组蛋白序列的鉴定结果

文件名：TF_and_best1_in_Ath.list（还是上面的那个转录因子预测文件）

文件内容如下：

图片.png

library(tidyverse)
# 1.读取数据
df <- read_tsv('TF_and_best1_in_Ath.list', comment = '#', col_names = F)
# 只保留前两列
df <- select(df, X1:X2)

然后是获得TF id和转录因子家族的对应关系

TF2term <- df
TF2term <- df %>% 
  select(term = X2) %>%
  mutate(TF = paste('TF', 1:nrow(TF2term), sep = '_')) %>%
  select(TF, term)

以该文件为基础，给所有转录因子的基因ID分配TF id

TF2gene <- df %>% 
  left_join(TF2term, by = c('X2' = 'term')) %>%
  select(TF, gene = X1)

富集分析

library(clusterProfiler)
options(stringsAsFactors = F)
genelist <- read.table('genelist.txt')$V1

TFenrich <- enricher(gene = genelist,
                     TERM2GENE = TF2gene, 
                     TERM2NAME = TF2term,
                     pvalueCutoff = 1,
                     qvalueCutoff = 1,
                     pAdjustMethod = 'BH')
# plot
dotplot(TFenrich, showCategory = 20)

图片.png

另外，你可以把富集结果转换成数据框输出或进一步绘图

TF_datafram <- as.data.frame(TFenrich)
write.table(TF_datafram, 'TF_enrich.result', sep = '\t', row.names = F, quota =F)

6、GO富集柱形图绘制

下载go_enrichment_full_all.txt文件

图片.png

将go_enrichment_full_all.txt中的GO.ID、Annotated（Description）、Count（GeneNumber）和Aspect（type）四列手动提取出来，得到文件leaf_vs_fruit.csv

图片.png

绘图

data=read.csv("leaf_vs_fruit.csv",header=T,stringsAsFactors = F)
1.按照qvalue升序排序，分别选出前20个BP,CC,MF的条目，由于enrichGO函数生成的数据框默认是按照qvalue升序排序，所以这里我们只用选取前二十个就行了
go_MF<-data[data$type=="molecular function",][1:20,]
go_CC<-data[data$type=="cellular component",][1:20,]
go_BP<-data[data$type=="biological process",][1:20,]
go_enrich_df<-data.frame(ID=c(go_BP$ID, go_CC$ID, go_MF$ID),
                         Description=c(go_BP$Description, go_CC$Description, go_MF$Description),
                         GeneNumber=c(go_BP$GeneNumber, go_CC$GeneNumber, go_MF$GeneNumber),
                         type=factor(c(rep("biological process", 20), rep("cellular component", 20),rep("molecular function",20)),levels=c("molecular function", "cellular component", "biological process")))
2.将GO_term设定为factor即可按照顺序输出
GO_term_order=factor(as.integer(rownames(go_enrich_df)),labels=go_enrich_df$Description)
library(ggplot2)
ggplot(data=go_enrich_df, aes(x=GO_term_order,y=GeneNumber, fill=type)) + geom_bar(stat="identity", width=0.8) + coord_flip() +  xlab("GO term") + ylab("Num of Genes") + theme_bw() + theme_classic() + theme(panel.border = element_rect(colour = "black", fill=NA, size=1))