作業(yè)要求

我們統(tǒng)一選擇p<0.05而且abs(log2FC)大于1的基因?yàn)轱@著差異表達(dá)基因集恬砂，對(duì)這個(gè)基因集用R包做KEGG/GO超幾何分布檢驗(yàn)分析。
然后把表達(dá)矩陣和分組信息分別作出cls和gct文件蓬痒，導(dǎo)入到GSEA軟件分析泻骤。
來(lái)源于生信技能樹(shù)：http://www.biotrainee.com/forum.php?mod=viewthread&tid=1750#lastpost

實(shí)驗(yàn)過(guò)程

1.差異基因篩選

我在轉(zhuǎn)錄組入門（7）：差異基因分析已經(jīng)完成了差異基因篩選，為了更好的銜接，我將上一步的代碼也一并寫入狱掂，完整流暢一些演痒，最后我們得到的是數(shù)據(jù)diff_gene_deseq2，包含了差異表達(dá)基因趋惨。(這里就不在詳細(xì)注釋這些代碼鸟顺，請(qǐng)看上一篇文章)

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上有詳細(xì)的說(shuō)明器虾，這個(gè)包的功能很強(qiáng)大讯嫂，我小白一個(gè)真的是整不明白，大致看了一些兆沙，不過(guò)還是有學(xué)習(xí)到很多欧芽，下面就開(kāi)始貼代碼。

2.1 安裝clusterProfiler

安裝clusterProfiler以及依賴的包挤悉，因?yàn)閭€(gè)人的電腦都是有差別的渐裸，所以我也不好說(shuō)，這樣的代碼就一定適合你装悲，因?yàn)樵谖覅⒖紕e人的時(shí)候昏鹃，就是出現(xiàn)了很多問(wèn)題，沒(méi)法安裝和載入這個(gè)包诀诊。具體問(wèn)題還是要具體分析洞渤，也不要那么容易放棄，稍微折騰一些属瓣，說(shuō)不定就能解決载迄。

# Bioconductor的包，安裝都是一個(gè)套路抡蛙，source一下护昧，bioLite一下，就差不多了粗截。
source("https://bioconductor.org/biocLite.R")
biocLite("clusterProfiler")
library(clusterProfiler)
# DOSE和DO.db這兩個(gè)包在我安裝的時(shí)候提示需要安裝惋耙，才能載入clusterProfiler，所以就直接安裝熊昌。
# 問(wèn)題是在我安裝的過(guò)程中绽榛，又提示好多依賴包沒(méi)法安裝，出現(xiàn)了權(quán)限的問(wèn)題婿屹，說(shuō)是目錄NOT PERMISSION灭美。
# 所以一氣之下，我就直接修改了R包的讀寫權(quán)限昂利，因?yàn)閭€(gè)人電腦届腐，也沒(méi)有什么特別重要的資料铁坎，
# 所以我就直接將相關(guān)的R包的目錄遞歸修改成777，這可是相當(dāng)危險(xiǎn)的操作犁苏，可不要隨意在服務(wù)器上進(jìn)行厢呵，后果自負(fù)哈。
# 平時(shí)個(gè)人電腦我都是以root身份進(jìn)行操作傀顾，一下在Ubuntu上以普通用戶的身份
# 經(jīng)常出現(xiàn)權(quán)限不足的提示，沒(méi)有辦法進(jìn)行操作碌奉，尤其是R短曾，非常的麻煩。
# 我主要修改了/usr/lib/R/library 和 /usr/local/lib這兩個(gè)目錄赐劣，全部遞歸修改權(quán)限為777嫉拐，折后貌似可以安裝成功。
biocLite("DOSE")
require(DOSE)
library(DO.db)

2.2 安裝構(gòu)建自己的基因組注釋數(shù)據(jù)

Biocouductor官網(wǎng)已經(jīng)擁有了構(gòu)建好的常用的19個(gè)注釋數(shù)據(jù)庫(kù)魁兼，包括了小鼠婉徘，人類和擬南芥等常用注釋數(shù)據(jù)，可以用安裝bioconductor包的方法來(lái)直接安裝和載入注釋數(shù)據(jù)咐汞，直接使用盖呼。

19個(gè)注釋數(shù)據(jù)庫(kù)

# 我們是小鼠數(shù)據(jù)，所以直接安裝載入就可以了化撕，當(dāng)然人類的也是一樣几晤。
# 人類的注釋數(shù)據(jù)
biocLite("org.Hs.eg.db")
library(org.Hs.eg.db)
# 小鼠的注釋數(shù)據(jù)
biocLite("org.Mm.eg.db")
library(org.Mm.eg.db)

• 如果沒(méi)有包括在這些注釋數(shù)據(jù)庫(kù)里面，那么就需要使用AnnotationHub這個(gè)包來(lái)構(gòu)建自己的OrgDb植阴。代碼如下:

# 這個(gè)包應(yīng)該是clusterProfiler自帶的蟹瘾，可以直接載入
library(AnnotationHub)
hub <- AnnotationHub()
# 可以用query（）函數(shù)來(lái)查找你要的物種注釋信息，這里我參考官網(wǎng)的內(nèi)容掠手，我查找的是番茄的注釋
# 選擇的格式是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 
# 下載注釋數(shù)據(jù)
sl <- hub[["AH55774"]] 

2.3 GO（Gene Ontology）分析

這里涉及到多種類型的ID轉(zhuǎn)換，我們常見(jiàn)的ENSEMBL喷鸽，ENTREZID這兩大類众雷，這里我在分析的時(shí)候發(fā)現(xiàn)，ENTREZID=kegg=ncbi-geneid魁衙，這三者有事相同的ID號(hào)报腔。jimmy博客有詳細(xì)的介紹，我進(jìn)行了適當(dāng)?shù)膮⒖迹?a target="_blank" rel="nofollow">http://www.bio-info-trainee.com/710.html剖淀。

• ID轉(zhuǎn)換函數(shù)介紹

# 看一下數(shù)據(jù)庫(kù)的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()函數(shù)進(jìn)行ID的轉(zhuǎn)換
# keys是原始的ID纯蛾，columns是轉(zhuǎn)換之后的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()函數(shù)可以轉(zhuǎn)換ID纵隔，得到的結(jié)果都是一樣的
anyid <- bitr(gene,fromType = "ENSEMBL",toType = c("GENENAME","SYMBOL","ENTREZID"),OrgDb = org.Mm.eg.db)

• enrichGO()函數(shù)進(jìn)行GO分析及畫圖

主要函數(shù)及參數(shù)：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:主要的分為三種炮姨，三個(gè)層面來(lái)闡述基因功能，生物學(xué)過(guò)程（BP）碰煌，細(xì)胞組分（CC）舒岸，分子功能（MF）；OrgDb:指定物種注釋數(shù)據(jù)芦圾；keytype:ID類型蛾派；pAdjustMethod:P值校正方法。

# 進(jìn)行g(shù)o分析
ego <- enrichGO(
  gene = row.names(diff_gene_deseq2),
  OrgDb = org.Mm.eg.db,
  keytype = "ENSEMBL",
  ont = "MF"
)
# 氣泡圖
dotplot(ego, font.size=5)
# 網(wǎng)絡(luò)圖
enrichMap(ego, vertex.label.cex=1.2, layout=igraph::layout.kamada.kawai)
# GO圖需要安裝額外的包
biocLite("topGO")
biocLite("Rgraphviz")
require(Rgraphviz)
plotGOgraph(ego)

氣泡圖

網(wǎng)絡(luò)圖

GO圖

關(guān)于這些圖的說(shuō)明个少，可以參考諾禾致源的微信文章

2.4 GSEA分析

基因集富集分析 (Gene Set Enrichment Analysis, GSEA) 的基本思想是使用預(yù)定義的基因集（通常來(lái)自功能注釋或先前實(shí)驗(yàn)的結(jié)果）洪乍，將基因按照在兩類樣本中的差異表達(dá)程度排序，然后檢驗(yàn)預(yù)先設(shè)定的基因集合是否在這個(gè)排序表的頂端或者底端富集夜焦】前模基因集合富集分析檢測(cè)基因集合而不是單個(gè)基因的表達(dá)變化，因此可以包含這些細(xì)微的表達(dá)變化茫经，預(yù)期得到更為理想的結(jié)果巷波。
參考資料：GSEA分析是什么鬼(上)和GSEA分析是什么鬼(下)。

# Gene Set Enrichment Analysis（GSEA）
# 獲取按照l(shuí)og2FC大小來(lái)排序的基因列表
genelist <- diff_gene_deseq2$log2FoldChange
names(genelist) <- rownames(diff_gene_deseq2)
genelist <- sort(genelist, decreasing = TRUE)
# GSEA分析（具體參數(shù)參考：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結(jié)果分析圖

2.5 KEGG（pathway）分析

KEGG將基因組信息和高一級(jí)的功能信息有機(jī)地結(jié)合起來(lái)卸伞，通過(guò)對(duì)細(xì)胞內(nèi)已知生物學(xué)過(guò)程的計(jì)算機(jī)化處理和將現(xiàn)有的基因功能解釋標(biāo)準(zhǔn)化抹镊，對(duì)基因的功能進(jìn)行系統(tǒng)化的分析。
KEGG的另一個(gè)任務(wù)是一個(gè)將基因組中的一系列基因用一個(gè)細(xì)胞內(nèi)的分子相互作用的網(wǎng)絡(luò)連接起來(lái)的過(guò)程瞪慧，如一個(gè)通路或是一個(gè)復(fù)合物髓考，通過(guò)它們來(lái)展現(xiàn)更高一級(jí)的生物學(xué)功能。
參考文章：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

# 轉(zhuǎn)換ID適合KEGG
x=bitr(rownames(diff_gene_deseq2),fromType = "ENSEMBL",toType = "ENTREZID", OrgDb = "org.Mm.eg.db")
# 獲取keggID
kegg <- x[,2]
# KEGG分析弃酌，在KEGG官網(wǎng)中氨菇，物種都有對(duì)應(yīng)的縮寫，小鼠mmu妓湘，其他的縮寫看官網(wǎng)：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結(jié)果轉(zhuǎn)換成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：最后查蓉，終于完成了轉(zhuǎn)入組入門，從小白慢慢的開(kāi)始入門榜贴，確實(shí)不容易豌研，中間有想過(guò)要放棄，真的太難唬党，沒(méi)有完全一樣的流程可以讓我參考鹃共，只能一遍看別人的，一遍自己摸索驶拱，慢慢的學(xué)習(xí)霜浴，我很慶幸自己堅(jiān)持了來(lái)了，走了一遍流程蓝纲，雖然沒(méi)有那樣的熟悉阴孟，但是卻是個(gè)極大的進(jìn)步晌纫。這里必須要感謝幾個(gè)大牛的幫助：Jimmy大神，徐洲更同學(xué)永丝，還有Y叔锹漱，主要參考了他們幾個(gè)人的博客，邊學(xué)習(xí)慕嚷，邊進(jìn)步哥牍，著實(shí)不易。