GEO-基因表达芯片分析（接上文
不要在同一个文件夹中处理两套数据
实验分组（Group）和探针注释（ids）
分组
library(stringr)
##生成Group的三种方法
if(F){
#1.有现成可以分组的列
Group = pd$group #group是分组列
}else if(F){
#2.第二种方法，眼睛数，自己生成
Group = rep(c("Disease","Normal"),each = 10)
}else if(T){
# 第三种方法，使用字符串处理的函数获取分组
k = str_detect(pd$title,"Normal");table(k)
Group = ifelse(k,"Normal","Disease")
}
##检查一下
data.frame(a=pd$title,
b=Group)
# 需要把Group转换成因子，并设置参考水平，指定levels，对照组在前，处理组在后
Group = factor(Group,levels = c("LGOSC","HGOSC"))
Group
##如果是多分组里面取2个分组(分别含有关键词‘A'和'B')
k1 = str_detect(pd$title,'A');table(k1)
k2 = str_detect(pd$title,'B');table(k2)
pd = pd[k1|k2,] ##变成2分类变量
探针注释的获取
探针注释：探针与基因的对应关系
##捷径
library(tinyarray)
find_anno(gpl_number) #辅助写出找注释的代码
#按照返回的代码提示下载↓
library(hgu133plus2.db);ids <- toTable(hgu133plus2SYMBOL)
##如果这个捷径没有成功，则按以下方法获取
注释来源：
- Bioconductor 的注释包
- GPL页面的表格文件解析：ID和Gene Symbol，下载后提取这两列即可用。
- 官网下载对应产品的注释表格
- 自主注释 教程：https://mp.weixin.qq.com/s/mrtjpN8yDKUdCSvSUuUwcA
不是所有的GPL都能找到注释
去除NA值：na.omit()
绘图
PCA图
-
sthda网站查找PCA图，获得代码及示例数据。
-
将自己的数据制作成示例数据的格式使用。
dat=as.data.frame(t(exp))
library(FactoMineR)
library(factoextra)
dat.pca <- PCA(dat, graph = FALSE)
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"
)
热图
- 看帮助文档的示例可以发现自定义热图的方法
g = names(tail(sort(apply(exp,1,sd)),1000)) ##取方差最大的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之间，超出此范围的数字显示极限颜色
)
LIMMA差异分析
library(limma)
design = model.matrix(~Group)
fit = lmFit(exp,design)
fit = eBayes(fit)
deg = topTable(fit,coef = 2,number = Inf)
- 差异分析后的数据整理
#1.加probe_id列，把行名变成一列
library(dplyr)
deg = mutate(deg,probe_id = rownames(deg))
#2.加上探针注释
ids = distinct(ids,symbol,.keep_all = T) ##去重复(随机法)
deg = inner_join(deg,ids,by="probe_id")
nrow(deg) #如果行数为0就是你找的探针注释是错的。
- 如果发现非特异性探针：1个探针对应多个基因，一般去掉该探针。
- 多个探针对应1个基因
- 随机去重，随机选其中之一，按照symbol去重复
- 保留行和/行平均值最大的探针
- 取多个探针的平均值
#3.加change列,标记上下调基因
logFC_t = 1
p_t = 0.05
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)
#4.加ENTREZID列，用于富集分析（symbol转entrezid，然后inner_join）
library(clusterProfiler)
library(org.Hs.eg.db)
s2e = bitr(deg$symbol, #Symbol转换Entrezid
fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db)#人类,注意物种
#其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
nrow(deg)
deg = inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
#SYMBOL与ENTREZID不是一对一的。
nrow(deg)
火山图
#火山图
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()
差异基因热图
# 表达矩阵行名替换为基因名
exp = exp[deg$probe_id,] ##表达矩阵按照deg的探针取子集
rownames(exp) = deg$symbol
diff_gene = deg$symbol[deg$change !="stable"]
n = exp[diff_gene,] #取出差异基因
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_cols = F
)
富集分析
- 富集分析：要求输入差异基因的Entrezid，任何类型的基因转换都会损失一些基因。
- KEGG数据库：8k+基因收录，系统分析基因功能、基因组信息数据库，有助于把基因及表达信息作为一个整体网络进行研究。
- GO数据库：10k+基因收录，涵盖生物学的三个方面：细胞组分、分子功能、生物过程。
library(clusterProfiler)
library(ggthemes)
library(org.Hs.eg.db)
library(dplyr)
library(ggplot2)
library(stringr)
library(enrichplot)
#(1)输入数据
gene_diff = deg$ENTREZID[deg$change != "stable"]
#(2)富集
ekk <- enrichKEGG(gene = gene_diff,organism = 'hsa')
ekk <- setReadable(ekk,OrgDb = org.Hs.eg.db,keyType = "ENTREZID") ##使结果易读，把Entrezid转成symbol
ego <- enrichGO(gene = gene_diff,OrgDb= org.Hs.eg.db,
ont = "ALL",readable = TRUE)
#setReadable和readable = TRUE都是把富集结果表格里的基因名称转为symbol，ont='ALL'是分析三个方面，如果只需要一种可以修改参数
class(ekk)
#KEGG一般看p.adj，不是p.value
#(3)可视化
barplot(ego, split = "ONTOLOGY") +
facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y")
barplot(ekk)
#或者是dotplot
# 更多资料---
# GSEA：https://www.yuque.com/docs/share/a67a180f-dd2b-4f6f-96c2-68a4b86fe862?#
# Y叔的书：http://yulab-smu.top/clusterProfiler-book/index.html
# GOplot：https://mp.weixin.qq.com/s/LonwdDhDn8iFUfxqSJ2Wew
-
富集分析结果解读：
View(ekk@result)- 衡量每个通路里的基因在差异基因里面是否足够多
- 结果看P.adj，一般不使用P.value
- GeneRatio：差异基因中有n个属于该通路/差异基因中有n个被数据库收录
- BgRatio：该通路有n个基因被收录/一共n个基因被收录
- 可视化：barplot/dotplot都可以，根据自己的数据灵活处理
-
富集不到的补救秘籍：
- 调整logFC、pvalue阈值，以改动差异基因数量。
- 不使用默认的padj（富集的），而用原始p值，在文章中说清。
- 换富集方法，GSEA也可以做KEGG富集。
- 调参数maxGSSize = 500，默认参数，表示500个基因以上的通路不考虑，可以调大。（但基因数量超多的通路，要想好生物学解释
问题数据和常见错误
数据问题
- 表达矩阵为空
- 表达矩阵不完整
- 表达矩阵被标准化过，需要从原始数据开始处理
- 表达矩阵有错误或异常值（未log的数据有负值等）
处理流程问题
- 用芯片流程分析转录组数据
- log过程（忘记log/多余的log
- 分组错误
- 探针注释错误
- id转换用错物种
不可抗力
- 找不到探针注释
- 数据错误以及找不到原始数据
复杂数据分析
多分组数据
本质上是多个两两分组的集合
一个对照组，多个实验组：可以对照组分别对应不同实验组进行二二分析。也可以实验组之间二二分析。PCA、热图都可以三个组放在一起，可以应用韦恩图做基因交集。火山图需要分开成两两分析。
多数据联合
-
分别分析：各自差异分析，差异基因取交集
-
先合并数据，然后差异分析
- 原则上（√）选择来自同一芯片平台的GSE
- 需要处理批次效应 Batch effect
- 不要选择一个全是处理组，一个全是对照组的数据去合并
批次效应↑
校正去批次：(这里也先标记一下，以后遇到了再做详细攻略)
limma::removeBatchEffect() sva::ComBat()
加权共表达网络（WGCNA）
图1：软阈值β的选择（0.9/0.85/0.8）
计算过程中一个可以调整的参数
-
R^2：无标度网络的拟合度/判定系数，评估拟合模型对观测数据的解释能力
-
R2越大，说明越接近无标度网络，选择使R2第一次达到0.8/0.85/0.9的β值。（β：软阈值，相关性矩阵向邻接矩阵转换的参数。
-
connectivity：连接度，反应节点的重要程度
mean connectivity：平均连通性，尽可能大。
二者中和。
无标度网络和随机网络
图2：基因模块化
图3.模块与性状之间的关联
将基因聚类成模块后，关注其与性状的关联，数字越大，越相关，一般选择最大值。
↑两种展示方式
图4：MM&GS
- GS代表模块里的每个基因与性状的相关性
- MM代表每个基因和所在模块之间的相关性，表示是否与模块的趋势一致。
- 希望GS和MM尽可能大
图5：TOM-拓扑重叠矩阵
- 希望对角线上出现红色框框
蛋白蛋白基因互作网络（PPI）
-
蛋白互作网络图-网页工具string；输入差异基因，输出一个ppi图，可以导出数据。
-
放入cytoscape进行网络可视化
-
寻找hub基因-使用cytohHubba插件
-
子网络-使用Mcode插件
Tinyarray包
小洁老师的超方便基因表达芯片处理流程包！
library(tinyarray)
packageVersion("tinyarray")
library(stringr)
geo = geo_download("GSE7305")
exp = geo$exp
exp = log2(exp+1)
boxplot(exp,las = 2)
pd = geo$pd
gpl_number = geo$gpl
# 分组信息
k = str_detect(pd$title,"Normal");table(k)
Group = ifelse(k,"Normal","Disease")
Group = factor(Group,levels = c("Normal","Disease"))
# 探针注释
find_anno(geo$gpl)
library(hgu133plus2.db);ids <- toTable(hgu133plus2SYMBOL)
head(ids)
#差异分析和它的可视化
dcp = get_deg_all(exp,Group,ids,entriz = F)
table(dcp$deg$change)
head(dcp$deg)
dcp$plots
library(ggplot2)
ggsave("deg.png",width = 15,height = 5)
#富集分析
deg = get_deg(exp,Group,ids)
genes = deg$ENTREZID[deg$change!="stable"]
head(genes)
#有可能因为网络问题报错
g = quick_enrich(genes,destdir = tempdir())
names(g)
#g是4个元素组成的列表，kk、go、kk.dot、go.dot，后两个为图
g[[1]][1:4,1:4]
library(patchwork)
g[[3]]+g[[4]]
ggsave("enrich.png",width = 12,height = 7)
最后马克一下
https://www.yuque.com/xiaojiewanglezenmofenshen/kzgwzl/pn0s1mmsaxocfynb?singleDoc# 《又一个有点难的探针注释》
https://www.yuque.com/xiaojiewanglezenmofenshen/kzgwzl/sv262capcgg9o8s5?singleDoc# 《一个有点难的探针注释》
https://www.jianshu.com/p/b94449be175a《分组凌乱怎么办》
引用自生信技能树课程！今天放了很多小洁老师的授课直出！有些东西还需要后续消化，先mark住~ 再次比心小洁老师爱你爱你！
5066
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
