非模式物種

step1:注釋
把本物種的蛋白序列用diamond到eggnog數(shù)據(jù)庫比對

##/home/lx_sky6/software/eggnog-mapper-2.1.8/emapper.py -i GCF_905147795.1_ilPieRapa1.1_protein.faa -m diamond -o eggnog --cpu 60

拿到eggnog的注釋結果記得吧以#開頭的行全部刪掉
step2:構建自己的物種數(shù)據(jù)庫

eggnogmapper<- read_delim(
  file = 'eggnog.emapper.annotations2',
  "\t",
  escape_double =FALSE,
  col_names=FALSE,
  comment = "#",trim_ws = TRUE)%>%
  dplyr::select(GID=X1,
                Gene_Symbol=X9,
                GO=X10,
                KO=X12,
                Pathway=X13,
                OG=X7,
                Gene_Name=X21)
gene_info <- dplyr::select(eggnogmapper,  GID, Gene_Name) %>%
  dplyr::filter(!is.na(Gene_Name))
# 提取GO信息
gene2go <- dplyr::select(eggnogmapper, GID, GO) %>%
  separate_rows(GO, sep = ',', convert = F) %>%
  filter(!is.na(GO)) %>%
  mutate(EVIDENCE = 'IEA')
gene2go <- subset(gene2go,grepl("^.*GO.*$",GO))


# 構建 OrgDB,參數(shù)十分玄學，請按照以下格式
library(AnnotationForge)
AnnotationForge::makeOrgPackage(gene_info=gene_info,
                                               go=gene2go,
                                               maintainer='yuantao <909474045@qq.com>',
                                               author='yuantao',
                                               outputDir="/home/lx_sky6/yt/rworkplace",
                                               tax_id=64459,
                                               genus='Pieris',                     #只改動這里
                                               species='rapae',               #和這里，進行命名
                                               goTable="go",
                                               version="1.0")
pkgbuild::build('/home/lx_sky6/yt/rworkplace/org.Prap.eg.db', dest_path = "/home/lx_sky6/yt/rworkplace/")       #這兩個路徑自己創(chuàng)建

step3:安裝注釋包進行注釋

install.packages("./org.Prapae.eg.db", repos = NULL,type = "source")

library(org.Prapae.eg.db)
library(readr)
fuji1_gene <- read_table2("../anita/fuji1.txt")
##fuji1:254個基因
de_gene <- pull(fuji1_gene, GID)
de_ego <- enrichGO(gene = de_gene,
                   OrgDb = org.Prapae.eg.db,
                   keyType = 'GID',
                   ont = 'ALL',
                   qvalueCutoff = 0.05,
                   pvalueCutoff = 0.05)
de_ego_df <- as.data.frame(de_ego)
enrichplot::dotplot(de_ego,showCategory = 20)


##run__KEGG富集 
pathway2gene <- dplyr::select(eggnogmapper, Pathway, GID) %>%
  separate_rows(Pathway, sep = ',', convert = F) %>%
  filter(str_detect(Pathway, 'ko')) %>%
  mutate(Pathway = str_remove(Pathway, 'ko'))
library(magrittr)
get_path2name <- function(){
  keggpathid2name.df <- clusterProfiler:::kegg_list("pathway")
  keggpathid2name.df[,1] %<>% gsub("path:map", "", .)
  colnames(keggpathid2name.df) <- c("path_id","path_name")
  return(keggpathid2name.df)
}
pathway2name <- get_path2name()
library(clusterProfiler)
de_ekp <- enricher(de_gene,
                   TERM2GENE = pathway2gene,
                   TERM2NAME = pathway2name,
                   pvalueCutoff = 0.05,
                   qvalueCutoff = 0.05)
de_ekp_df <- as.data.frame(de_ekp)
head(de_ekp_df)
enrichplot::dotplot(de_ekp, showCategory = 20)
##導出數(shù)據(jù)
write.table(de_ego_df,file = '/home/lx_sky6/yt/rworkplace/914_GO.csv',sep = ',',quote = FALSE)

write.table(de_ekp_df,file = '/home/lx_sky6/yt/rworkplace/914_Kegg.csv',sep = ',',quote = FALSE)