笔记回顾：
1.GEO数据库视频学习笔记（芯片数据下载和数据读取）
2.GEO数据库视频学习笔记（ID转换）
3.GEO数据库视频学习笔记（了解你的表达矩阵）
这一讲的视频地址是：[https://www.bilibili.com/video/BV1is411H7Hq?p=7]
（一）差异分析
上面三个视频笔记，记录了如何获得矩阵、如何过滤探针、如何获得分组信息，具备了这些，就可以进行差异分析了。这里Jimmy大神用的是limma包（当然差异分析还有很多其他的包可以用），参考文章：用limma对芯片数据做差异分析
首先做分组矩阵：
> 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
# 对照是用来被比的
当你做比较矩阵的时候，如果你搞不清到底是谁和谁在比，请看文章：
【r<-差异分析】当使用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应该是个正数，因为表达量经过处理后上升了。可以看到limma差异分析矩阵里logFC是5.78017，是个正数，那么就应该没什么问题。如果你发现表达上调的基因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分析
在视频里，关于GO和KEGG这两个部分讲的很快，没有涉及到具体的代码，需要的童鞋可以参考我之前写的文章：
RNA-seq练习 第三部分（DEseq2筛选差异表达基因,可视化）
做完差异分析后，实际上你得到的差异基因远远小于总基因数，你只是知道它上调下调了多少实际上是不够的，因为有些基因虽然在处理后有变化，但你还需要知道这些基因哪些是非常重要的。这时就需要GSEA分析。关于这一块，视频讲的有些凌乱。。。对于GESA软件，大神也木有细讲。可以参考我之前的笔记：GSEA学习笔记，里面都是非常非常详细的代码，也详细的讲了怎么看GSEA的图以及简单的原理，还有如何使用软件进行GSEA分析。这里只根据GEO数据贴上代码：
> 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') #这里你可以根据上一行代码得到的data.frame里任何一个通路代码来进行展示
#图就不放了
上面你会发现我分别用了GO和KEGG的结果来展示GESA，这就是根据富集结果来展示哪些基因（通路）是重要的一种手段。
网友评论