作業(yè)要求
我們統(tǒng)一選擇p<0.05而且abs(log2FC)大于1的基因為顯著差異表達(dá)基因集已烤，對這個基因集用R包做KEGG/GO超幾何分布檢驗分析敞嗡。
然后把表達(dá)矩陣和分組信息分別作出cls和gct文件密任，導(dǎo)入到GSEA軟件分析扼仲。
來源于生信技能樹：http://www.biotrainee.com/forum.php?mod=viewthread&tid=1750#lastpost
實驗過程
1.差異基因篩選
我在轉(zhuǎn)錄組入門（7）：差異基因分析已經(jīng)完成了差異基因篩選虎锚，為了更好的銜接壁查，我將上一步的代碼也一并寫入，完整流暢一些璧尸，最后我們得到的是數(shù)據(jù)diff_gene_deseq2咒林，包含了差異表達(dá)基因。(這里就不在詳細(xì)注釋這些代碼爷光，請看上一篇文章)
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ì)的說明，這個包的功能很強(qiáng)大蛀序，我小白一個真的是整不明白欢瞪，大致看了一些，不過還是有學(xué)習(xí)到很多徐裸，下面就開始貼代碼遣鼓。
2.1 安裝clusterProfiler
安裝clusterProfiler以及依賴的包，因為個人的電腦都是有差別的重贺，所以我也不好說骑祟，這樣的代碼就一定適合你，因為在我參考別人的時候气笙，就是出現(xiàn)了很多問題次企，沒法安裝和載入這個包。具體問題還是要具體分析潜圃，也不要那么容易放棄缸棵，稍微折騰一些，說不定就能解決秉犹。
# Bioconductor的包蛉谜，安裝都是一個套路，source一下崇堵，bioLite一下型诚，就差不多了。
source("https://bioconductor.org/biocLite.R")
biocLite("clusterProfiler")
library(clusterProfiler)
# DOSE和DO.db這兩個包在我安裝的時候提示需要安裝鸳劳，才能載入clusterProfiler狰贯，所以就直接安裝。
# 問題是在我安裝的過程中赏廓，又提示好多依賴包沒法安裝涵紊，出現(xiàn)了權(quán)限的問題，說是目錄NOT PERMISSION幔摸。
# 所以一氣之下摸柄，我就直接修改了R包的讀寫權(quán)限，因為個人電腦既忆，也沒有什么特別重要的資料驱负，
# 所以我就直接將相關(guān)的R包的目錄遞歸修改成777嗦玖，這可是相當(dāng)危險的操作，可不要隨意在服務(wù)器上進(jìn)行跃脊，后果自負(fù)哈宇挫。
# 平時個人電腦我都是以root身份進(jìn)行操作，一下在Ubuntu上以普通用戶的身份
# 經(jīng)常出現(xiàn)權(quán)限不足的提示酪术，沒有辦法進(jìn)行操作器瘪，尤其是R，非常的麻煩绘雁。
# 我主要修改了/usr/lib/R/library 和 /usr/local/lib這兩個目錄橡疼，全部遞歸修改權(quán)限為777，折后貌似可以安裝成功咧七。
biocLite("DOSE")
require(DOSE)
library(DO.db)
2.2 安裝構(gòu)建自己的基因組注釋數(shù)據(jù)
Biocouductor官網(wǎng)已經(jīng)擁有了構(gòu)建好的常用的19個注釋數(shù)據(jù)庫衰齐，包括了小鼠，人類和擬南芥等常用注釋數(shù)據(jù)继阻，可以用安裝bioconductor包的方法來直接安裝和載入注釋數(shù)據(jù)，直接使用废酷。
# 我們是小鼠數(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)
- 如果沒有包括在這些注釋數(shù)據(jù)庫里面墨辛，那么就需要使用AnnotationHub這個包來構(gòu)建自己的OrgDb。代碼如下:
# 這個包應(yīng)該是clusterProfiler自帶的趴俘，可以直接載入
library(AnnotationHub)
hub <- AnnotationHub()
# 可以用query（）函數(shù)來查找你要的物種注釋信息睹簇，這里我參考官網(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)換太惠，我們常見的ENSEMBL，ENTREZID這兩大類疲憋，這里我在分析的時候發(fā)現(xiàn)凿渊，ENTREZID=kegg=ncbi-geneid，這三者有事相同的ID號缚柳。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ù)庫的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:主要的分為三種堵幽，三個層面來闡述基因功能狗超，生物學(xué)過程（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)
關(guān)于這些圖的說明，可以參考諾禾致源的微信文章
2.4 GSEA分析
基因集富集分析 (Gene Set Enrichment Analysis, GSEA) 的基本思想是使用預(yù)定義的基因集（通常來自功能注釋或先前實驗的結(jié)果）溃肪，將基因按照在兩類樣本中的差異表達(dá)程度排序免胃，然后檢驗預(yù)先設(shè)定的基因集合是否在這個排序表的頂端或者底端富集”棺基因集合富集分析檢測基因集合而不是單個基因的表達(dá)變化羔沙，因此可以包含這些細(xì)微的表達(dá)變化，預(yù)期得到更為理想的結(jié)果厨钻。
參考資料：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分析（具體參數(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")
2.5 KEGG（pathway）分析
KEGG將基因組信息和高一級的功能信息有機(jī)地結(jié)合起來，通過對細(xì)胞內(nèi)已知生物學(xué)過程的計算機(jī)化處理和將現(xiàn)有的基因功能解釋標(biāo)準(zhǔn)化夯膀，對基因的功能進(jìn)行系統(tǒng)化的分析诗充。
KEGG的另一個任務(wù)是一個將基因組中的一系列基因用一個細(xì)胞內(nèi)的分子相互作用的網(wǎng)絡(luò)連接起來的過程，如一個通路或是一個復(fù)合物诱建，通過它們來展現(xiàn)更高一級的生物學(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)中俺猿，物種都有對應(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)
PS：最后押袍，終于完成了轉(zhuǎn)入組入門诵冒，從小白慢慢的開始入門，確實不容易伯病，中間有想過要放棄造烁，真的太難，沒有完全一樣的流程可以讓我參考午笛，只能一遍看別人的惭蟋，一遍自己摸索，慢慢的學(xué)習(xí)药磺，我很慶幸自己堅持了來了告组，走了一遍流程，雖然沒有那樣的熟悉癌佩，但是卻是個極大的進(jìn)步木缝。這里必須要感謝幾個大牛的幫助：Jimmy大神便锨，徐洲更同學(xué)，還有Y叔我碟，主要參考了他們幾個人的博客放案，邊學(xué)習(xí)，邊進(jìn)步矫俺，著實不易吱殉。
