基因的富集分析

通过RNA_seq差异表达分析等获得一组基因，想要进一步了解这些基因可能参与哪些生物学过程，具有哪些功能，一般就会用到基因的富集分析了。现在小麦基因注释信息已经可以获得，目前来看包括两套。一个是TGACv1，另一个就是IWGSC1.0。需要用到一个叫“clusterProfiler”的R包。下面用到的数据集来自TGACv1。

1、安装

source("https://bioconductor.org/biocLite.R")biocLite("clusterProfiler")biocLite("topGO")install.packages("DOSE")

需要的文件有3个，格式见下图。

library("clusterProfiler")library(ggplot2)library(stringr)setwd("/Users/mashengwei/Desktop/TGACv1/WGCNA/")####cellular_component#####gene <- read.csv(gene_list.txt",header=FALSE,sep="\t")gene <- as.factor(gene$V1)term2gene <- read.csv("cellular_component_Go_term_gene.txt",header=TRUE,sep="\t")term2name <- read.csv("GO_name.txt",header=TRUE,sep="\t")x <- enricher(gene,TERM2GENE=term2gene,TERM2NAME=term2name)out_file <- paste("TSG_nearest_CG_cellular_component_enricher.out.txt",sep ="\t")write.csv(x,out_file)dotplot(x)ggsave(filename="dotplot_cellular_component.png",dpi=600)dev.off()####molecular_function#####gene <- read.csv("gene_list2.txt",header=FALSE,sep="\t")gene <- as.factor(gene$V1)term2gene <- read.csv("molecular_function_Go_term_gene.txt",header=TRUE,sep="\t")term2name <- read.csv("GO_name.txt",header=TRUE,sep="\t")x <- enricher(gene,TERM2GENE=term2gene,TERM2NAME=term2name)out_file <- paste("molecular_function_enricher.out.txt",sep ="\t")write.csv(x,out_file)dotplot(x)ggsave(filename="dotplot_molecular_function.png",dpi=600)dev.off()####biological_process#####gene <- read.csv("gene_list3.txt",header=FALSE,sep="\t")gene <- as.factor(gene$V1)term2gene <- read.csv("biological_process_Go_term_gene.txt",header=TRUE,sep="\t")term2name <- read.csv("GO_name.txt",header=TRUE,sep="\t")x <- enricher(gene,TERM2GENE=term2gene,TERM2NAME=term2name)out_file <- paste("process_enricher.out.txt",sep ="\t")write.csv(x,out_file)dotplot(x) + scale_y_discrete(labels=function(y) str_wrap(y, width=10))ggsave(filename="dotplot_biological_process.png",dpi=600)dev.off()

将GO信息换成KEGG pathway信息即可进行KEGG pathway分析。