转录组RNAseq基因注释/富集分析与功能分类

3.转录组 | 富集分析(clusterProfiler)

2019-04-25  本文已影响96人  pomela
基因富集分析 (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 转录组学习八(功能富集分析)


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

上一篇下一篇

猜你喜欢

热点阅读