純搬運(yùn)工调违，不錯(cuò)的教程墅垮，留著自學(xué)
【內(nèi)容出自轉(zhuǎn)錄組入門（8）：差異基因結(jié)果注釋 - 簡書 (jianshu.com)】
作業(yè)要求
我們統(tǒng)一選擇p<0.05而且abs(log2FC)大于1的基因?yàn)轱@著差異表達(dá)基因集西饵，對這個(gè)基因集用R包做KEGG/GO超幾何分布檢驗(yàn)分析。
然后把表達(dá)矩陣和分組信息分別作出cls和gct文件兄旬，導(dǎo)入到GSEA軟件分析捺疼。
來源于生信技能樹：http://www.biotrainee.com/forum.php?mod=viewthread&tid=1750#lastpost
實(shí)驗(yàn)過程
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ì)的說明逞敷，這個(gè)包的功能很強(qiáng)大，我小白一個(gè)真的是整不明白灌侣，大致看了一些推捐，不過還是有學(xué)習(xí)到很多，下面就開始貼代碼侧啼。
2.1 安裝clusterProfiler
安裝clusterProfiler以及依賴的包牛柒，因?yàn)閭€(gè)人的電腦都是有差別的愕秫，所以我也不好說，這樣的代碼就一定適合你焰络，因?yàn)樵谖覅⒖紕e人的時(shí)候，就是出現(xiàn)了很多問題符喝，沒法安裝和載入這個(gè)包闪彼。具體問題還是要具體分析，也不要那么容易放棄协饲，稍微折騰一些畏腕，說不定就能解決。
# Bioconductor的包茉稠，安裝都是一個(gè)套路描馅，source一下，bioLite一下而线，就差不多了铭污。
source("https://bioconductor.org/biocLite.R")
biocLite("clusterProfiler")
library(clusterProfiler)
# DOSE和DO.db這兩個(gè)包在我安裝的時(shí)候提示需要安裝，才能載入clusterProfiler膀篮，所以就直接安裝嘹狞。
# 問題是在我安裝的過程中，又提示好多依賴包沒法安裝誓竿，出現(xiàn)了權(quán)限的問題磅网，說是目錄NOT PERMISSION。
# 所以一氣之下筷屡，我就直接修改了R包的讀寫權(quán)限涧偷，因?yàn)閭€(gè)人電腦，也沒有什么特別重要的資料毙死，
# 所以我就直接將相關(guān)的R包的目錄遞歸修改成777燎潮，這可是相當(dāng)危險(xiǎn)的操作，可不要隨意在服務(wù)器上進(jìn)行扼倘，后果自負(fù)哈跟啤。
# 平時(shí)個(gè)人電腦我都是以root身份進(jìn)行操作，一下在Ubuntu上以普通用戶的身份
# 經(jīng)常出現(xiàn)權(quán)限不足的提示唉锌，沒有辦法進(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ù)庫秃症，包括了小鼠候址，人類和擬南芥等常用注釋數(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è)包來構(gòu)建自己的OrgDb驶赏。代碼如下:
# 這個(gè)包應(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這兩大類五续，這里我在分析的時(shí)候發(fā)現(xiàn)，ENTREZID=kegg=ncbi-geneid龄恋，這三者有事相同的ID號返帕。jimmy博客有詳細(xì)的介紹，我進(jìn)行了適當(dāng)?shù)膮⒖迹?a target="_blank">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:主要的分為三種，三個(gè)層面來闡述基因功能煮落，生物學(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ù)定義的基因集（通常來自功能注釋或先前實(shí)驗(yàn)的結(jié)果），將基因按照在兩類樣本中的差異表達(dá)程度排序害驹，然后檢驗(yàn)預(yù)先設(shè)定的基因集合是否在這個(gè)排序表的頂端或者底端富集鞭呕。基因集合富集分析檢測基因集合而不是單個(gè)基因的表達(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ì)算機(jī)化處理和將現(xiàn)有的基因功能解釋標(biāo)準(zhǔn)化咕娄，對基因的功能進(jìn)行系統(tǒng)化的分析。
KEGG的另一個(gè)任務(wù)是一個(gè)將基因組中的一系列基因用一個(gè)細(xì)胞內(nèi)的分子相互作用的網(wǎng)絡(luò)連接起來的過程珊擂，如一個(gè)通路或是一個(gè)復(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)入組入門昼激，從小白慢慢的開始入門庇绽，確實(shí)不容易，中間有想過要放棄橙困，真的太難瞧掺，沒有完全一樣的流程可以讓我參考，只能一遍看別人的凡傅，一遍自己摸索辟狈，慢慢的學(xué)習(xí)，我很慶幸自己堅(jiān)持了來了夏跷，走了一遍流程哼转，雖然沒有那樣的熟悉，但是卻是個(gè)極大的進(jìn)步槽华。這里必須要感謝幾個(gè)大牛的幫助：Jimmy大神壹蔓，徐洲更同學(xué)，還有Y叔猫态，主要參考了他們幾個(gè)人的博客佣蓉，邊學(xué)習(xí)，邊進(jìn)步亲雪，著實(shí)不易勇凭。
作者：lxmic
鏈接：http://www.reibang.com/p/4910d7cec5c8
來源：簡書
著作權(quán)歸作者所有。商業(yè)轉(zhuǎn)載請聯(lián)系作者獲得授權(quán)义辕，非商業(yè)轉(zhuǎn)載請注明出處套像。
