GEO數(shù)據(jù)庫視頻學(xué)習(xí)筆記（差異分析、可視化李命、GSEA）

筆記回顧：
1.GEO數(shù)據(jù)庫視頻學(xué)習(xí)筆記（芯片數(shù)據(jù)下載和數(shù)據(jù)讀取）
2.GEO數(shù)據(jù)庫視頻學(xué)習(xí)筆記（ID轉(zhuǎn)換）
3.GEO數(shù)據(jù)庫視頻學(xué)習(xí)筆記（了解你的表達矩陣）

這一講的視頻地址是：[https://www.bilibili.com/video/BV1is411H7Hq?p=7]

（一）差異分析

上面三個視頻筆記松邪，記錄了如何獲得矩陣坞琴、如何過濾探針、如何獲得分組信息逗抑，具備了這些剧辐，就可以進行差異分析了归形。這里Jimmy大神用的是limma包（當(dāng)然差異分析還有很多其他的包可以用）鞭缭，參考文章：用limma對芯片數(shù)據(jù)做差異分析

首先做分組矩陣：

> dim(exprSet)
[1] 18821     6
> library(limma)
#這一步就是告訴design，哪幾個樣品是對照压语，哪幾個樣品是處理褂傀，1代表是忍啤，0代表不是
> design <- model.matrix(~0+factor(group_list))
> colnames(design)=levels(factor(group_list))
> rownames(design)=colnames(exprSet)
> design
             control Vemurafenib
control1           1           0
control2           1           0
control3           1           0
Vemurafenib1       0           1
Vemurafenib2       0           1
Vemurafenib3       0           1
attr(,"assign")
[1] 1 1
attr(,"contrasts")
attr(,"contrasts")$`factor(group_list)`
[1] "contr.treatment"

再做比較矩陣：

#這個矩陣聲明，我們要把Vemurafenib組跟control組進行差異分析比較
> contrast_matrix <- makeContrasts(Vemurafenib-control,levels=design)
> contrast_matrix
             Contrasts
Levels        Vemurafenib - control
  control                        -1
  Vemurafenib                     1
# 對照是用來被比的

當(dāng)你做比較矩陣的時候仙辟，如果你搞不清到底是誰和誰在比同波，請看文章：
【r<-差異分析】當(dāng)使用limma時，它在比較什么

然后就是差異分析了：

> library(limma)
# step1
> fit <- lmFit(exprSet,design)
# step2
> fit2 <- contrasts.fit(fit,contrast_matrix)
> fit2 <- eBayes(fit2)
# step3
> tempOutput = topTable(fit2,coef=1,n=Inf)
> nrDEG = na.omit(tempOutput)
> head(nrDEG)
          logFC   AveExpr         t      P.Value    adj.P.Val        B
CD36   5.780170  7.370282  79.74674 1.224250e-16 2.304160e-12 26.74942
DUSP6 -4.212683  9.106625 -62.45965 1.822336e-15 1.339567e-11 25.00170
DCT    5.633027  8.763220  61.57004 2.135222e-15 1.339567e-11 24.88864
SPRY2 -3.801663  9.726468 -53.97076 9.143161e-15 4.302086e-11 23.80028
MOXD1  3.263063 10.171635  47.09632 4.109399e-14 1.318367e-10 22.58677
ETV4  -3.843247  9.667077 -47.00035 4.202860e-14 1.318367e-10 22.56798
> dif <- tempOutput[tempOutput[,"P.Value"]<0.01,] #你可以只保存p<0.01的基因叠国，也可以不篩選未檩，這里不是強制的
> write.csv(dif,file="Differential_limma_matrix.csv")

重點到了！Ｋ诤浮冤狡！如何檢查你的差異分析有沒有問題？项棠？悲雳？

先看一下剛得到的差異分析矩陣：

我們就以第一個基因為檢查例子：CD36

> exprSet[rownames(exprSet)=='CD36',]
    control1     control2     control3 Vemurafenib1 Vemurafenib2 Vemurafenib3 
     4.54610      4.40210      4.49239     10.25060     10.21480    10.31570

在表達矩陣里，CD36這個基因在對照組里大概是4沾乘，在處理組里是10左右怜奖。那么它的logFC應(yīng)該是個正數(shù)，因為表達量經(jīng)過處理后上升了翅阵⊥崃幔可以看到limma差異分析矩陣里logFC是5.78017，是個正數(shù)掷匠，那么就應(yīng)該沒什么問題滥崩。如果你發(fā)現(xiàn)表達上調(diào)的基因logFC是個負值，那么可能是你之前的分組讹语，或者比較矩陣沒有做對钙皮。

簡單的畫個熱圖吧，選差異分析列表里的前25個基因：

> library(pheatmap)
> getgene = head(rownames(dif),25)
> getgene_matrix=exprSet[getgene,]
> getgene_matrix=t(scale(t(getgene_matrix)))
> pheatmap(getgene_matrix)

（二）火山圖

先畫一個最丑火山圖：

> plot(nrDEG$logFC,-log10(nrDEG$P.Value))

太丑了，可以給它美化一下短条，但是Jimmy在這個視頻里并沒有提到如何美化导匣，所以我用了自己之前用過的代碼：

> library(EnhancedVolcano)
> library(airway)
> EnhancedVolcano(nrDEG,
                 lab = rownames(nrDEG),
                 x = 'logFC',
                 y = 'P.Value',
                 xlim = c(-6, 6),
                 title = 'Vemurafenib versus control',
                 pCutoff = 0.01,
                 FCcutoff = 2,
                 transcriptPointSize = 1.5,
                 transcriptLabSize = 3.0,
                 col=c('black', 'blue', 'green', 'red1'),
                 colAlpha = 1,
                 legend=c('NS','logFC','P.value',
                          'P.value & logFC'),
                 legendPosition = 'right',
                 legendLabSize = 14,
                 legendIconSize = 5.0,
                 gridlines.major = FALSE,
                 gridlines.minor = FALSE
)

是不是比前面那個黑白火山圖強多了？

（三）GO,KEGG,GSEA分析

在視頻里茸时，關(guān)于GO和KEGG這兩個部分講的很快贡定，沒有涉及到具體的代碼，需要的童鞋可以參考我之前寫的文章：
RNA-seq練習(xí) 第三部分（DEseq2篩選差異表達基因,可視化）

做完差異分析后可都，實際上你得到的差異基因遠遠小于總基因數(shù)缓待，你只是知道它上調(diào)下調(diào)了多少實際上是不夠的，因為有些基因雖然在處理后有變化渠牲，但你還需要知道這些基因哪些是非常重要的旋炒。這時就需要GSEA分析。關(guān)于這一塊签杈，視頻講的有些凌亂瘫镇。。芹壕。對于GESA軟件汇四，大神也木有細講。可以參考我之前的筆記：GSEA學(xué)習(xí)筆記踢涌，里面都是非常非常詳細的代碼通孽，也詳細的講了怎么看GSEA的圖以及簡單的原理，還有如何使用軟件進行GSEA分析睁壁。這里只根據(jù)GEO數(shù)據(jù)貼上代碼：

> geneList <- nrDEG$logFC
> names(geneList) <- gene.df$ENTREZID
> geneList_new <- sort(geneList,decreasing = T)
#go_result <- gseGO(geneList  = geneList_new,
#                   ont = "BP", 
#                   OrgDb = org.Hs.eg.db,
#                   minGSSize    = 10,
#                   maxGSSize = 500,
#                   pvalueCutoff = 1)
#go_result <- setReadable(go_result,OrgDb = org.Hs.eg.db)
#go_result_df <- as.data.frame(go_result)
#gseaplot(go_result,geneSetID = c("GO:0000727"))
> kegg <-gseKEGG(geneList  = geneList_new,
               organism     = 'hsa',
               nPerm        = 1000,
               minGSSize    = 10,
               maxGSSize = 500,
               pvalueCutoff = 1,
               verbose      = FALSE)
> kegg_result <- as.data.frame(kegg)
> gseaplot(kegg,geneSetID = 'hsa03030') #這里你可以根據(jù)上一行代碼得到的data.frame里任何一個通路代碼來進行展示
#圖就不放了

上面你會發(fā)現(xiàn)我分別用了GO和KEGG的結(jié)果來展示GESA背苦，這就是根據(jù)富集結(jié)果來展示哪些基因（通路）是重要的一種手段。

參考文章潘明；生信技能樹：差異分析得到的結(jié)果注釋一文就夠

GEO數(shù)據(jù)庫視頻學(xué)習(xí)筆記（差異分析潘酗、可視化然眼、GSEA）

GEO數(shù)據(jù)庫視頻學(xué)習(xí)筆記（差異分析、可視化琅攘、GSEA）

（一）差異分析

重點到了！Ｋ诤浮冤狡！如何檢查你的差異分析有沒有問題？项棠？悲雳？

（二）火山圖

（三）GO,KEGG,GSEA分析