转录组入门（8）：差异基因结果注释

2017-08-16 本文已影响2642人 lxmic

作业要求

我们统一选择p<0.05而且abs(log2FC)大于1的基因为显著差异表达基因集，对这个基因集用R包做KEGG/GO超几何分布检验分析。
然后把表达矩阵和分组信息分别作出cls和gct文件，导入到GSEA软件分析。
来源于生信技能树：http://www.biotrainee.com/forum.php?mod=viewthread&tid=1750#lastpost

实验过程

1.差异基因筛选

我在转录组入门（7）：差异基因分析已经完成了差异基因筛选，为了更好的衔接，我将上一步的代码也一并写入，完整流畅一些，最后我们得到的是数据diff_gene_deseq2，包含了差异表达基因。(这里就不在详细注释这些代码，请看上一篇文章)

require(DESeq2)
control1 <- read.table("~/disk2/data/rna-seq/matrix/SRR3589959.count", sep="\t", col.names = c("gene_id","control1"))
control2 <- read.table("~/disk2/data/rna-seq/matrix/SRR3589961.count", sep="\t", col.names = c("gene_id","control2")) 
rep1 <- read.table("~/disk2/data/rna-seq/matrix/SRR3589960.count", sep="\t", col.names = c("gene_id","akap951")) 
rep2 <- read.table("~/disk2/data/rna-seq/matrix/SRR3589962.count", sep="\t",col.names = c("gene_id","akap952"))
raw_count <- merge(merge(control1, control2,by="gene_id"),merge(rep1,rep2, by="gene_id"))
raw_count_filt <- raw_count[-48823:-48825,]
raw_count_filter <- raw_count_filt[-1:-2,]
ENSEMBL <- gsub("\\.\\d*", "", raw_count_filter$gene_id) 
row.names(raw_count_filter) <- ENSEMBL
raw_count_filter <- raw_count_filter[ ,-1]
condition <- factor(c(rep("control",2),rep("akap95",2)), levels = c("control","akap95"))
countData <- raw_count_filter[,1:4]
colData <- data.frame(row.names=colnames(raw_count_filter), condition)
dds <- DESeqDataSetFromMatrix(countData, colData, design= ~ condition)
head(dds)
dds2 <- DESeq(dds)
resultsNames(dds2)
res <- results(dds2)
summary(res)
table(res$padj<0.05)
res <- res[order(res$padj),]
diff_gene_deseq2 <-subset(res,padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1))

2.GO/KEGG分析及GSEA

我们主要用到的就是Y叔的R包：clusterProfiler包，github上有详细的说明，这个包的功能很强大，我小白一个真的是整不明白，大致看了一些，不过还是有学习到很多，下面就开始贴代码。

2.1 安装clusterProfiler

安装clusterProfiler以及依赖的包，因为个人的电脑都是有差别的，所以我也不好说，这样的代码就一定适合你，因为在我参考别人的时候，就是出现了很多问题，没法安装和载入这个包。具体问题还是要具体分析，也不要那么容易放弃，稍微折腾一些，说不定就能解决。

# Bioconductor的包，安装都是一个套路，source一下，bioLite一下，就差不多了。
source("https://bioconductor.org/biocLite.R")
biocLite("clusterProfiler")
library(clusterProfiler)
# DOSE和DO.db这两个包在我安装的时候提示需要安装，才能载入clusterProfiler，所以就直接安装。
# 问题是在我安装的过程中，又提示好多依赖包没法安装，出现了权限的问题，说是目录NOT PERMISSION。
# 所以一气之下，我就直接修改了R包的读写权限，因为个人电脑，也没有什么特别重要的资料，
# 所以我就直接将相关的R包的目录递归修改成777，这可是相当危险的操作，可不要随意在服务器上进行，后果自负哈。
# 平时个人电脑我都是以root身份进行操作，一下在Ubuntu上以普通用户的身份
# 经常出现权限不足的提示，没有办法进行操作，尤其是R，非常的麻烦。
# 我主要修改了/usr/lib/R/library 和 /usr/local/lib这两个目录，全部递归修改权限为777，折后貌似可以安装成功。
biocLite("DOSE")
require(DOSE)
library(DO.db)

2.2 安装构建自己的基因组注释数据

Biocouductor官网已经拥有了构建好的常用的19个注释数据库，包括了小鼠，人类和拟南芥等常用注释数据，可以用安装bioconductor包的方法来直接安装和载入注释数据，直接使用。

19个注释数据库

# 我们是小鼠数据，所以直接安装载入就可以了，当然人类的也是一样。
# 人类的注释数据
biocLite("org.Hs.eg.db")
library(org.Hs.eg.db)
# 小鼠的注释数据
biocLite("org.Mm.eg.db")
library(org.Mm.eg.db)

如果没有包括在这些注释数据库里面，那么就需要使用AnnotationHub这个包来构建自己的OrgDb。代码如下:

# 这个包应该是clusterProfiler自带的，可以直接载入
library(AnnotationHub)
hub <- AnnotationHub()
# 可以用query（）函数来查找你要的物种注释信息，这里我参考官网的内容，我查找的是番茄的注释
# 选择的格式是OrgDb，所以我们选择AH55774
query（hub, "Solanum lycopersicum"）
## AnnotationHub with 2 records
## # snapshotDate(): 2017-04-25 
## # $dataprovider: Inparanoid8, ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/
## # $species: Solanum lycopersicum
## # $rdataclass: Inparanoid8Db, OrgDb
## # additional mcols(): taxonomyid, genome, description,
## #   coordinate_1_based, maintainer, rdatadateadded, preparerclass, tags,
## #   rdatapath, sourceurl, sourcetype 
## # retrieve records with, e.g., 'object[["AH10593"]]' 

##          title                               
##   AH10593 | hom.Solanum_lycopersicum.inp8.sqlite
##  AH55774 | org.Solanum_lycopersicum.eg.sqlite 
# 下载注释数据
sl <- hub[["AH55774"]]

2.3 GO（Gene Ontology）分析

这里涉及到多种类型的ID转换，我们常见的ENSEMBL，ENTREZID这两大类，这里我在分析的时候发现，ENTREZID=kegg=ncbi-geneid，这三者有事相同的ID号。jimmy博客有详细的介绍，我进行了适当的参考：http://www.bio-info-trainee.com/710.html。

ID转换函数介绍

# 看一下数据库的ID类型
keytype(org.Mm.eg.db)
##  [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS"
##  [6] "ENTREZID"     "ENZYME"       "EVIDENCE"     "EVIDENCEALL"  "GENENAME"    
##  [11] "GO"           "GOALL"        "IPI"          "MGI"          "ONTOLOGY"    
##  [16] "ONTOLOGYALL"  "PATH"         "PFAM"         "PMID"         "PROSITE"     
##  [21] "REFSEQ"       "SYMBOL"       "UNIGENE"      "UNIPROT" 
# Jimmy推荐的是使用select()函数进行ID的转换
# keys是原始的ID，columns是转换之后的ID，keytype是要指定的原始ID类型
gene <- row.names(diff_gene_deseq2)
tansid <- select(org.Mm.eg.db,keys = gene,columns = c("GENENAME","SYMBOL","ENTREZID"),keytype = "ENSEMBL")
## ENSEMBL                                                GENENAME SYMBOL ENTREZID
## 1 ENSMUSG00000003309 adaptor protein complex AP-1, mu 2 subunit  Ap1m2    11768
## 2 ENSMUSG00000046323    developmental pluripotency-associated 3  Dppa3    73708
## 3 ENSMUSG00000001123       lectin, galactose binding, soluble 9 Lgals9    16859
## 4 ENSMUSG00000018569                                  claudin 7  Cldn7    53624
## 5 ENSMUSG00000023906                                  claudin 6  Cldn6    54419
## 6 ENSMUSG00000000184                                  cyclin D2  Ccnd2    12444
# 此外还有bitr()函数可以转换ID，得到的结果都是一样的
anyid <- bitr(gene,fromType = "ENSEMBL",toType = c("GENENAME","SYMBOL","ENTREZID"),OrgDb = org.Mm.eg.db)

enrichGO()函数进行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:差异基因ID；ont:主要的分为三种，三个层面来阐述基因功能，生物学过程（BP），细胞组分（CC），分子功能（MF）；OrgDb:指定物种注释数据；keytype:ID类型；pAdjustMethod:P值校正方法。

# 进行go分析
ego <- enrichGO(
  gene = row.names(diff_gene_deseq2),
  OrgDb = org.Mm.eg.db,
  keytype = "ENSEMBL",
  ont = "MF"
)
# 气泡图
dotplot(ego, font.size=5)
# 网络图
enrichMap(ego, vertex.label.cex=1.2, layout=igraph::layout.kamada.kawai)
# GO图需要安装额外的包
biocLite("topGO")
biocLite("Rgraphviz")
require(Rgraphviz)
plotGOgraph(ego)

气泡图

网络图

GO图

关于这些图的说明，可以参考诺禾致源的微信文章

2.4 GSEA分析

基因集富集分析 (Gene Set Enrichment Analysis, GSEA) 的基本思想是使用预定义的基因集（通常来自功能注释或先前实验的结果），将基因按照在两类样本中的差异表达程度排序，然后检验预先设定的基因集合是否在这个排序表的顶端或者底端富集。基因集合富集分析检测基因集合而不是单个基因的表达变化，因此可以包含这些细微的表达变化，预期得到更为理想的结果。
参考资料：GSEA分析是什么鬼(上)和GSEA分析是什么鬼(下)。

# Gene Set Enrichment Analysis（GSEA）
# 获取按照log2FC大小来排序的基因列表
genelist <- diff_gene_deseq2$log2FoldChange
names(genelist) <- rownames(diff_gene_deseq2)
genelist <- sort(genelist, decreasing = TRUE)
# GSEA分析（具体参数参考：https://mp.weixin.qq.com/s/p-n5jq5Rx2TqDBStS2nzoQ）
gsemf <- gseGO(genelist,
               OrgDb = org.Mm.eg.db,
               keyType = "ENSEMBL",
               ont="MF"
)
# 查看大致信息
head(gsemf)
# 画出GSEA图
gseaplot(gsemf, geneSetID="GO:0000977")

GSEA结果分析图

2.5 KEGG（pathway）分析

KEGG将基因组信息和高一级的功能信息有机地结合起来，通过对细胞内已知生物学过程的计算机化处理和将现有的基因功能解释标准化，对基因的功能进行系统化的分析。
KEGG的另一个任务是一个将基因组中的一系列基因用一个细胞内的分子相互作用的网络连接起来的过程，如一个通路或是一个复合物，通过它们来展现更高一级的生物学功能。
参考文章：http://blog.sciencenet.cn/blog-364884-779116.html
KEGG物种缩写：http://www.genome.jp/kegg/catalog/org_list.html
GO和KEGG输出文件解读：http://www.bio-info-trainee.com/370.html

# 转换ID适合KEGG
x=bitr(rownames(diff_gene_deseq2),fromType = "ENSEMBL",toType = "ENTREZID", OrgDb = "org.Mm.eg.db")
# 获取keggID
kegg <- x[,2]
# KEGG分析，在KEGG官网中，物种都有对应的缩写，小鼠mmu，其他的缩写看官网：http://www.genome.jp/kegg/catalog/org_list.html
ekk <- enrichKEGG(kegg, keyType = "kegg",organism = "mmu", pvalueCutoff = 0.05, pAdjustMethod = "BH", qvalueCutoff = 0.1)
head(summary(ekk))
# 气泡图
dotplot(ekk, font.size=5)
# 将GO/KEGG结果转换成CSV格式输出
write.csv(as.data.frame(ekk),"KEGG-enrich.csv",row.names =F)
write.csv(as.data.frame(ego),"GO-enrich.csv",row.names =F)

KEGG分析可视化

PS：最后，终于完成了转入组入门，从小白慢慢的开始入门，确实不容易，中间有想过要放弃，真的太难，没有完全一样的流程可以让我参考，只能一遍看别人的，一遍自己摸索，慢慢的学习，我很庆幸自己坚持了来了，走了一遍流程，虽然没有那样的熟悉，但是却是个极大的进步。这里必须要感谢几个大牛的帮助：Jimmy大神，徐洲更同学，还有Y叔，主要参考了他们几个人的博客，边学习，边进步，着实不易。