eggNOGmapper注释结果
具体如何使用eggNOGmapper进行注释在这里不详细赘述
我们先拿到对应的注释结果
eggmapper注释结果
然后我们要将其中的.emapper2.emapper.annotations 后缀的文件转移出来，将其放入annotation文件夹中。
# 设置源文件夹和目标文件夹的路径
source_folder <- "~\\eggmapperannoteation\\result"
target_folder <- "~\\eggmapperannoteation\\annotation"
#使用时请补全自己的路径
# 获取源文件夹中的所有文件
files <- list.files(source_folder, pattern = "\\.annotations$", full.names = TRUE)
# 创建目标文件夹（如果不存在）
dir.create(target_folder, recursive = TRUE, showWarnings = FALSE)
# 遍历源文件列表，并将文件复制到目标文件夹中
for (file in files) {
file.copy(file, target_folder)
}
# 输出完成提示信息
cat("文件复制完成！")
文件复制完毕
R中读入数据
rm(list=ls())
setwd("C:\\Users\\60522\\OneDrive\\桌面\\eggmapperannoteation\\annotation")
#读入注释结果
eggmapper_annotations<- read.table(file = "Alcanivorax_jadensis_C17.faa.emapper2.emapper.annotations",skip = 4
,sep ="\t")
下载感兴趣的kegg通路信息
进入KEGG官网通路下载页面
下载通路的所有KO号
点击进入，然后点击pathway entry
点击All links中的KEGG ORTHOLOGY
然后得到整个通路涉及的所有酶的对应KO号
KO号
将其复制保存到一个txt文件中，如M00720.txt
M00720
注：直接从网站上复制过去，ko号信息中间只是一堆空格和不是tab，所以建议一个个复制或者在excel中排好后复制到txt中，否者后续R读取时会出错
#读入感兴趣的通路信息
M00720<- read.table(file = "info/M00720.txt",sep = "\t",header = T)
使用stringr包来检测注释结果中是否含有对应通路的基因
result <- c(1:length(M00720$ID))*0
for (j in c(1:length(eggmapper_annotations$V1))){
for(k in c(1:length(M00720$ID))){
if (str_detect(string = str_to_lower(eggmapper_annotations$V9[j]),pattern = str_to_lower(M00720$ID[k]))) {
result[M00720$ID[k]] <- 1
break
}
}
}
批量执行
library(stringr)
rm(list=ls())
setwd("C:\\Users\\60522\\OneDrive\\桌面\\eggmapperannoteation\\annotation")
file_name <- dir(pattern = "annotations")
for (i in c(1:length(file_name))){
file_name[i] <- str_remove(string = file_name[i],pattern = ".faa.emapper2.emapper.annotations")
}
#读入感兴趣的通路信息
M00720<- read.table(file = "info/M00720.txt",sep = "\t",header = T)
result <- data.frame(matrix(data=0,nrow = length(file_name),ncol = length(M00720$ID)))
colnames(result) <- M00720$ID
rownames(result) <- file_name
for(i in c(1:length(file_name))){
#依次读入结果
eggmapper_annotations <- read.table(file = paste(file_name[i],".faa.emapper2.emapper.annotations",sep = "")
,sep ="\t")
for (j in c(1:length(eggmapper_annotations$V1))){
for(k in c(1:length(M00720$ID))){
if (str_detect(string = str_to_lower(eggmapper_annotations$V9[j]),pattern = str_to_lower(M00720$ID[k]))) {
cat(file_name[i]," ",M00720$ID[k],'\n')
result[file_name[i],M00720$ID[k]] <- 1
break
}
}
}
}
pheatmap(result,
color = colorRampPalette(c("white","#00FA9A"))(1000),
cluster_rows = F,clustering_method = "average",
cluster_cols = F,
fontsize = 2.5,
cellwidth = 3.5,
cellheight = 3.5,
display_numbers = F,
fontsize_number = 3,
angle_col=c(45)
)
结果
这时结果的横坐标是对应的KO号，可以使用下述代码将其更换为对应的注释信息
colnames(result) <- M00720$Definition #在这里把列名更改
pheatmap(result,
color = colorRampPalette(c("white","#00FA9A"))(1000),
cluster_rows = F,clustering_method = "average",
cluster_cols = F,
fontsize = 2.5,
cellwidth = 3.5,
cellheight = 3.5,
display_numbers = F,
fontsize_number = 3,
angle_col=c(45)
)
将结果投射到进化树上
使用chiplot的可视化工具TVBOT (chiplot.online)
具体使用方法作者在B站有讲，在该工具网站有跳转链接！
首先我们需要对结果进行输出
write.table(x = result,file = "result_tvBOT.txt")
导出后找到该文件所在地方，然后新建一个xlsx文件选择数据--自文本，进行导入数据。
转化为xlsx step 1
转化为xlsx step 2
转化为xlsx step 3
最后效果
最后参照tvBOT的官方教程微调即可实现数据可视化
效果图（条形图为另外加的）
缺点：tvBOT工具画出来后每个基因的注释无法显示，只有一堆散点图，需要后续一个个添加
解决方法：导出成pdf文件，然后在R中也导出一个无进化树版本的热图，这两个图的横坐标顺序是一致的，因而可以在AI中将注释信息移植过去！
使用pathview包绘制KEGG通路图
library(pathview)
KO_count <- colSums(result)/length(result[,1])#统计所有基因组出现的比例
KO_ID <- colnames(result)
KO_PATH <- data.frame(KO_count=KO_count)
rownames(KO_PATH) <- KO_ID
pathview(gene.data = KO_PATH,pathway.id = "00720",species = "ko")
#通常情况下我们还要关注一些稀缺类群是具有的独特基因，需要计算所有物种中出现的可能性
tmp <- KO_count
tmp[tmp>0] <- 1
KO_PATH_2 <- data.frame(KO_count=KO_count,KO_exite=tmp)
pathview(gene.data = KO_PATH_2,pathway.id = "00720",species = "ko")
注：在species中务必选择 ”ko“ 这个类型，否者无法运行。此外该包会从kegg下载数据，因而必须在有网络的条件使用
仅绘制统计比例的通路图
两组数据同时展示图
左边较暗代表比例，右边代表所有基因组中有无。这样可以快速帮助我们找到我们的基因组中有什么基因，大概比例多少！
补充说明：程序中涉及了三个包，其中的stringr 和 pheatmap两个包直接使用install.packages（）安装即可；但是pathview需要借助bioconductor进行安装。
install.packages("stringr")
install.packages("pheatmap")
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install(version = "3.18")#4.3版本的R对应3.18；4.2版本的R对应的3.16的版本
#安装好后再通过BiocManager安装pathview
BiocManager::install("pathview")
