基因富集分析 (gene set enrichment analysis)：在一组基因或蛋白中找到过表达的基因或蛋白。

常用的注释数据库：
①GO( The Gene Ontology Consortium): 描述基因的层级关系
②KEGG( Kyoto Encyclopedia of Genes and Genomes)：提供pathway的数据库
富集分析的方法：
①ORA(over-repressentation analysis):目前用得最多
②FCS(funcyional class scoring)：有一个可视化软件GSEA
③PT(pathway topology)：难度比较大

本次富集分析的工具是clusterProfiler，支持上面提到的三类方法。

1.准备数据集和软件

①数据准备

rm(list = ls())
options(stringsAsFactors = T)
library(DESeq2)
control <- read.table("SRR3589956_matrix.count",
                      sep="\t", col.names = c("gene_id","control"))
rep1 <- read.table("SRR3589957_matrix.count",
                   sep="\t", col.names = c("gene_id","rep1"))
rep2 <- read.table("SRR3589958_matrix.count",
                   sep="\t",col.names = c("gene_id","rep2"))
raw_count <- merge(merge(control, rep1, by="gene_id"), rep2, by="gene_id")
raw_count_filt <- raw_count[-1:-5,]
ENSEMBL <- gsub("(.*?)\\.\\d*?_\\d", "\\1", raw_count_filt$gene_id)
row.names(raw_count_filt) <- ENSEMBL

delta_mean <- abs(mean(raw_count_filt$rep1) - mean(raw_count_filt$rep2))

sampleNum <- length(raw_count_filt$control)
sampleMean <- mean(raw_count_filt$control)
control2 <- integer(sampleNum)

for (i in 1:sampleNum){
  if(raw_count_filt$control[i] < sampleMean){
    control2[i] <- raw_count_filt$control[i] + abs(raw_count_filt$rep1[i] - raw_count_filt$rep2[i])
  }
  else{
    control2[i] <- raw_count_filt$control[i] + rpois(1,delta_mean)
  }
}

raw_count_filt$control2 <- control2

library(DESeq2)
countData <- raw_count_filt[,2:5]
condition <- factor(c("control","KD","KD","control"))
dds <- DESeqDataSetFromMatrix(countData, DataFrame(condition), design= ~ condition )
nrow(dds)
dds <- dds[ rowSums(counts(dds)) > 1, ]
nrow(dds)

##使用DESeq进行差异表达分析
dds <- DESeq(dds)
res <- results(dds)
mcols(res, use.names = TRUE)
summary(res)

###准备富集分析的数据集，根据padj < 0.05 且Log2FoldChange的绝对值大于1的标准。
deseq2.sig <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)

②安装包下载注释数据 （rog.HS.eg.db)

source("https://bioconductor.org/biocLite.R")
biocLite("clusterProfiler")
biocLite("AnnotationHub")
library(AnnotationHub)
ah <- AnnotationHub()
head(unique(ah$species))
ah[ah$species == 'Homo sapiens' & ah$rdataclass == 'OrgDb']

目标数据库

org.hs <- ah[['AH61777']]
##若是在rstudio中下载不下来，可以将网址复制到网页中下载，之后再加载进来
loadDb("org.Hs.eg.db.sqlite")
org.hs <- loadDb("org.Hs.eg.db.sqlite")

##对数据库进行简单的探索
unique(ah$dataprovider)    #查看数据来源 
unique(ah$species)         #查看数据库中有多少个物种的注释信息
unique(ah$rdataclass)      #查看注释信息的数据存放格式，FaFile表示fasta文件； TxDb表示转录组数据库； OrgDb,用于基因ID、基因名、GO、KEGG ID之间的相互映射
ah$species[which(ah$species == "Arabidopsis thaliana")]   #筛选目标物种，查找目标物种是否存在，例如“拟南芥”
ah[ah$species == 'Arabidopsis thaliana' & ah$rdataclass == 'OrgDb']   #筛选目标物种，查找目标物种、目标数据库是否存在，例如“拟南芥的OrgDb库”

2.GO富集分析

GO分析主要的函数是：

enrichGO(gene, OrgDb, keyType = "ENTREZID", ont = "MF",
  pvalueCutoff = 0.05, pAdjustMethod = "BH", universe, qvalueCutoff = 0.2,
  minGSSize = 10, maxGSSize = 500, readable = FALSE, pool = FALSE)

##主要关注的参数如下：
gene: 差异表达的基因。不推荐使用"A1BG"这种类型的命名，请用setReadable转换。
OrgDb： 物种注释数据库，一般是org开头。
keytpe: ID类型。
ont： 有三个层次，分别是BP（Biological Process）, CC（Cellular Component）, MF（Molecular Function）。一个基因的功能可以从生物学过程、所属细胞部分和分子功能三个角度定义。
pAdjustMethod：p值矫正的方法。

提供相应的数据即可：

ego <- enrichGO(
  gene = row.names(deseq2.sig),
  OrgDb = org.hs,
  keyType = "ENSEMBL",
  ont = "MF"
)

head(ego)

可视化

dotplot(ego,font.size=10)
barplot(ego,showCategory=20,drop=T)
ego.MF <- enrichGO(gene = row.names(deseq2.sig), OrgDb = org.hs, ont='MF',pAdjustMethod = 'BH',pvalueCutoff = 0.05, qvalueCutoff = 0.2,keyType = 'ENSEMBL')
plotGOgraph(ego.MF)

dotplot

barplot

plotGOgraph

3.KEGG富集分析

KEGG分析主要的函数是：

enrichKEGG(gene, organism = "hsa", keyType = "kegg", pvalueCutoff = 0.05,
  pAdjustMethod = "BH", universe, minGSSize = 10, maxGSSize = 500,
  qvalueCutoff = 0.2, use_internal_data = FALSE)

主要参数：
gene： 基因名，要和keyType对应
organism: 需要参考 [http://www.genome.jp/kegg/catalog/org_list.html](https://link.jianshu.com/?t=http://www.genome.jp/kegg/catalog/org_list.html)， 人类是hsa
keyType: 基因的命名方式， "kegg", 'ncbi-geneid', 'ncib-proteinid' and 'uniprot'选择其一

library(clusterProfiler)
gene_list <- mapIds(org.hs, keys = row.names(deseq2.sig),
                       column = "ENTREZID", keytype = "ENSEMBL" )


kk <- enrichKEGG(gene_list, organism="hsa",
                 keyType = "ncbi-geneid",
                 pvalueCutoff=0.05, pAdjustMethod="BH",
                 qvalueCutoff=0.1)
head(summary(kk))

红色字符串表示 一个基因可以有多个功能（KEGG注释），当然一个KEGG注释下也有可以多个基因

可视化

dotplot(kk)

dotplot

4.GSEA分析

用clusterProfiler的gseGO或GSEA函数分析，后者可以自定义输入数据：

gseGO(geneList, ont = "BP", OrgDb, keyType = "ENTREZID", exponent = 1,
  nPerm = 1000, minGSSize = 10, maxGSSize = 500, pvalueCutoff = 0.05,
  pAdjustMethod = "BH", verbose = TRUE, seed = FALSE, by = "fgsea")
主要参数：
geneList： 排序数据， 可以根据log2foldchange, 也可以是pvalues
nPerm： 重抽取次数
minGSSize： 每个基因集的最小数目maxGSSize： 用于测试的基因注释最大数目

genelist <- sig.deseq2$log2FoldChange
names(genelist) <- rownames(sig.deseq2)
genelist <- sort(genelist, decreasing = TRUE)
gsemf <- gseGO(genelist,
      OrgDb = org.hs,
      keyType = "ENSEMBL",
      ont="MF"
      )
head(gsemf)

可视化

gseaplot(gsemf, geneSetID="GO:0004871")

GSEA

参考：
01 转录组入门(8): 富集分析
02 GEO数据挖掘小尝试：（三）利用clusterProfiler进行富集分析
03 转录组学习八（功能富集分析）

写在最后 | 不管是这一部分中的富集分析还是之前的差异表达分析，虽然参考前辈们写出的攻略便捷了很多，参看这些攻略会认为这部分的分析是满满的套路，但不跑一遍是不知道其中隐藏着哪些细节上的问题。同时这些也只是很粗浅的一些分析，仍然需要更进一步的学习。