ixxmu
commented
5 months ago 生信技能树-day9 GEO数据挖掘-PCA、差异分析 by 生信菜鸟团
从我们生信技能树历年的几千个马拉松授课学员里面募集了一些优秀的创作者，某种意义来说是传承了我们生信技能树的知识整理和分享的思想！
今天的是三周合计15天的数据挖掘授课学员一点一滴整理的授课知识点笔记哦，还有互动练习题哈，欢迎大家点击文末的阅读原文去关注我们学员的公众号哦！
PCA图、top1000基因热图
rm(list = ls())
load(file = "step2output.Rdata")
#输入数据：exp和Group
PCA 图
dat=as.data.frame(t(exp))
#根据sthda上的示例数据，更改自己的数据，需要转置后转为dataframe
Principal Component Analysis
http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials
library(FactoMineR)
library(factoextra)
dat.pca <- PCA(dat, graph = FALSE)
#对dat做PCA分析
fviz_pca_ind(dat.pca,
geom.ind = "point", # show points only (nbut not "text")
col.ind = Group, # color by groups
palette = c("#00AFBB", "#E7B800"),
addEllipses = TRUE, # Concentration ellipses
legend.title = "Groups"
)
#画PCA图，使用pca分析后的数据dat.pca
top 1000 sd 热图
#离散基因热图，top1000表达基因，只是看一下，不用放文章里
g = names(tail(sort(apply(exp,1,sd)),1000))
#对exp遍历，取每一行的sd，sort排序后取后1000个，再提取基因名字（向量名字）
n = exp[g,]
library(pheatmap)
annotation_col = data.frame(row.names = colnames(n),
Group = Group)
pheatmap(n,
show_colnames =F,
show_rownames = F,
annotation_col=annotation_col,
scale = "row", #按行标准化，只保留行内差别，不保留行间差别，会把数据范围缩放到大概-5~5之间
breaks = seq(-3,3,length.out = 100) #设置色带分布范围为-3~3之间，超出此范围的数字显示极限颜色
)
#pheatmap的帮助文档非常详细，有其他需要可以去查看帮助文档
# 关于scale的进一步学习：zz.scale.R
差异分析
rm(list = ls())
load(file = "step2output.Rdata")
#load上一步做完后得到用于差异分析的结果表格
#需要输入数据：exp、ids、group，使用limma包根据贝叶斯检验原理进行差异分析
library(limma)
design = model.matrix(~Group)
#根据分组信息生成一个模型矩阵
fit = lmFit(exp,design)
#线性拟合函数
fit = eBayes(fit)
#贝叶斯检验
deg = topTable(fit,coef = 2,number = Inf)
#提取差异分析结果：coef = 2指design的第二列，number = infinity指提取全部
得到一个probe_id和对应的logFC、pValue的表格，但还需要和symbol还有entrezid连接在一起
代码如下：
加probe_id列，把行名变成一列，防止行名丢失
library(dplyr)
deg = mutate(deg,probe_id = rownames(deg))
#mutate合并，deg和新的一列，新列名称叫probe_id，内容是rownames.deg
加上探针注释
probe_id和基因symbol不是一对一的关系，是多对一或者一对多的关系，因此需要去重
一个探针对应多个基因：非特异性探针，在表格中可以看到，则需要直接去除
一个基因对应多个探针：相当于一个基因测了好多遍
处理方法：随机去重/保留行和或行平均值最大的探针/取多个探针的平均值
ids = distinct(ids,symbol,.keep_all = T)
#此处是随机去重的方法，其他去重方式在zz.去重方式.R
deg = inner_join(deg,ids,by="probe_id")
#用inner_join取交集并把差异分析的结果deg和之前的id—symbol表格连接在一起
nrow(deg)
## [1] 20824
#检查一下，如果行数为0就是你找的探针注释是错的。
加change列,标记上下调基因
logFC_t = 1
p_t = 0.05
#设置logFC和pValue的阈值
#使用ifelse两次或者casewhen判断down、up、stable并输出成新的一列change
k1 = (deg$P.Value < p_t)&(deg$logFC < -logFC_t)
k2 = (deg$P.Value < p_t)&(deg$logFC > logFC_t)
deg = mutate(deg,change = ifelse(k1,"down",ifelse(k2,"up","stable")))
table(deg$change)
##
## down stable up
## 579 19621 624
火山图
library(ggplot2)
ggplot(data = deg, aes(x = logFC, y = -log10(P.Value))) +
geom_point(alpha=0.4, size=3.5, aes(color=change)) +
scale_color_manual(values=c("blue", "grey","red"))+
geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",linewidth=0.8) +
geom_hline(yintercept = -log10(p_t),lty=4,col="black",linewidth=0.8) +
theme_bw()
#使用ggplot的geom_point画火山图，vline和hline画阈值的线
差异基因热图
# 表达矩阵行名替换为基因名，分为两步：
exp = exp[deg$probe_id,]
#按deg中的symbol列的内容在exp中按行取子集
rownames(exp) = deg$symbol
#把exp中行名改为deg的symbol列，此时已经是一一对应的
diff_gene = deg$symbol[deg$change !="stable"]
#取出change列不是stable的基因symbol
n = exp[diff_gene,]
#按有差异的基因symbol在exp中按行取子集，即为有差异的基因的logFC和pValue，赋值到数据框n中，用于画差异基因热图
library(pheatmap)
annotation_col = data.frame(group = Group)
rownames(annotation_col) = colnames(n)
pheatmap(n,show_colnames =F,
show_rownames = F,
scale = "row",
#cluster_cols = F,
annotation_col=annotation_col,
breaks = seq(-3,3,length.out = 100)
)
#如果差异基因的聚类分组还是错乱的，则加cluster_col = F
#如果加了还是错乱的，去看：小洁老师的语雀/分组聚类的热图
加ENTREZID列，用于后续的富集分析
entrezid是富集分析指定用的，需要symbol转entrezid，然后inner_join
使用clusterProfiler中的bitr函数实现，另外数据库根据物种不同
library(clusterProfiler)
library(org.Hs.eg.db)
s2e = bitr(deg$symbol,
fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db)#人类,注意物种
一部分基因没匹配上是正常的。<30%的失败都没事。
其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb nrow(deg)
看剩下数量，如果只有几十说明有问题
deg = inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
#把差异基因的表和entrezid通过innerjoin连接在一起，用于后面的富集分析
多了几行少了几行都正常，SYMBOL与ENTREZID不是一对一的
nrow(deg)
## [1] 20827
#再看看还有几行，然后保存
save(exp,Group,deg,logFC_t,p_t,file = "step4output.Rdata")
❝from 生信技能树
文末友情宣传
强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：
https://mp.weixin.qq.com/s/jJfwC77CsoPimQJTIbkI2Q