生信小白菜之GEO芯片数据分析流程--附画图代码

最新推荐文章于 2024-05-17 23:29:48 发布
最新推荐文章于 2024-05-17 23:29:48 发布
阅读量1.1k 收藏 11
点赞数 30
文章标签： 数据分析 数据挖掘 r语言 数据库 笔记
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/yangyulu1998/article/details/136941593
版权

title: “GEO data analysis”
author: “yuluyang”
date: “2024-03-22”

生信技能树数据挖掘课程笔记～小洁老师授课

主要内容：

  • 数据分组的内容
  • 关键词的分组和多分组比较
  • id map报错的原因及解决方法
  • 基因组的注释流程
  • 数据的行列互换及方差数值
  • 画图示例代码

示例数据

library(GEOquery)
eSet = getGEO("GSE7305", destdir = '.', getGPL = F)
library(stringr)
# 紧接数据下载后
# 第一步：提取表达矩阵exp
eSet = eSet[[1]]  # 要先打破列表的框架方便后续提取信息
exp <- exprs(eSet) 
range(exp) # 看数据范围决定是否需要log，是否有负值，异常值
exp = log2(exp+1)
# 第二步：提取临床信息
pd <- pData(eSet)
## 要记得自己命名了哪些内容，'exp'是表达矩阵、'pd'是临床/分组信息

可以自己上GEO筛选一个芯片数据集，然后按以下流程复现一遍
只要命名方式不改变，后续作图代码可以直接copy然后run

提取后如何处理分组信息

if(F){    # 括号里是F，即不运行，这里只是示例有这种方法可供选择
  # 第一种方法，有现成的可以用来分组的列
}else if(F){
  # 第二种方法，眼睛数，自己生成
  Group = rep(c("Disease","Normal"),each = 10) # 这个数据是很整齐排好的，前10个为`disease`，后10个为`normal`
  # 如果两个组个数不一样的话咋整？
}else if(T){
  # 第三种方法，使用字符串处理的函数获取分组 
  k = str_detect(pd$title,"Normal");table(k)  # 选中信息这一列，找一个能区分组的关键词"Normal"，然后`str_detect`分成`true`和`false`两组
  ## 记得检查`table`的结果
  Group = ifelse(k,"Normal","Disease")  # 对`true`、`false`结果进行分组分配
}
## 想起什么名字都可以
## 分组关系一定要对应清楚

# 需要把Group转换成因子，并设置参考水平，指定levels，对照组在前，处理组在后
Group = factor(Group,levels = c("Normal","Disease")) ## `factor()`转化为因子型向量，`levels=`指定因子的水平向量
Group
## `levels=`如果不指定的话，默认是按照字母顺序排列的

推荐使用第三种，优先考虑代码梳理
factor因子型向量是R中一类特殊的向量类别
作用是给其他类别元素进行分组、计数，相当于统计学中分类变量，分可枚举【有序】和类别型【无序】两类

补充：如果两个组样本数量不一致
# 例如 对照组有10个，样本组15个
rep(c("control","disease"),times=c(10,15))
## 要记得检查生成的向量是不是先10个`control`再15个`disease`
rep(c("control","disease"),c(10,15)) 
## `times=`可以省略
## 不可以用`each=`
补充：如何检查分组设置与临床信息是否对应
# 新建一个数据框
data.frame(a=pd$title,
           b=Group)
## 两列对照着看
补充：多分组里面取两个分组
# 举🌰，现编一个数据
pd$test=paste0(rep(c("a","b","c","d"),each=5),rep(c(1:5),times=4)) # 新增一列'test'，分成a、b、c、d 4类，每一类5个
k1=str_detect(pd$test,"a");table(k1) # "a"为分组关键词
k2=str_detect(pd$test,"b");table(k2) # "b"为分组关键词
pd=pd[k1|k2,];pd # 再把这两列取子集出来

# 再上述取一遍
##  k = str_detect(pd$,"");table(k)  
## Group = ifelse(k,"","")

探针注释来源

  • Bioconductor的注释包（通常可以用代码直接下载）
  • GPL页面的表格文件解析（手动下载到工作目录然后按列读取）
  • 官网下载对应产品的注释表（部分芯片是定制的，需要到对应平台官网查询，也有可能找不到）
  • 自主注释（哒咩❌费劲）

用代码获取探针注释

# 推荐：用代码下载
library(tinyarray) 
find_anno(gpl_number)  
library(hgu133plus2.db);ids <- toTable(hgu133plus2SYMBOL) 
## 可能会报错，提示你"use ids <- AnnoProbe::idmap('GPL570')...please try different 'type'"
## '?AnnoProbe::idmap'查看帮助文档，有哪些不同type

## 当你成功获取注释，`environment`会出现`ids`
## 检查`ids`列名是否为'probe_id'和'symbol'，如不是，需手动修改

# 有些探针不对应具体的基因，需要去掉
## 去除空缺值
ids=na.omit(ids)

save(exp,Group,ids,file = "step2output.Rdata")

自己处理数据的时候，按分析的内容分段式以Rdata存放好

另一种可能的报错
# 返回 no annotation available in Bioconductor and AnnoProbe
# 那么就需要在GEO网页上找到对应的GPL平台，打开该平台的表格找到ID、symbol等信息
# 下载文件再导入

遇到这种有难度的芯片就参考一下别人的代码

https://www.yuque.com/xiaojiewanglezenmofenshen/kzgwzl/sv262capcgg9o8s5

补充：其他获取探针注释的方法
# 方法1 BioconductorR包(最常用，已全部收入find_anno里面，不用看啦) 
if(F){
  gpl_number #看看编号是多少
  #http://www.bio-info-trainee.com/1399.html #在这里搜索，找到对应的R包
  library(hgu133plus2.db)
  ls("package:hgu133plus2.db") #列出R包里都有啥
  ids <- toTable(hgu133plus2SYMBOL) #把R包里的注释表格变成数据框
}
# 方法2 读取GPL网页的表格文件，按列取子集
##https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570
# 方法3 官网下载注释文件并读取
# 方法4 自主注释，了解一下
#https://mp.weixin.qq.com/s/mrtjpN8yDKUdCSvSUuUwcA

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"
)

以上仅为作图代码示例
详细教程可参考STHDA网站，我后续也会进一步整理

http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials/

heatmap 画图

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之间，超出此范围的数字显示极限颜色
         )

差异分析——画火山图

library(limma)
design = model.matrix(~Group)
fit = lmFit(exp,design)
fit = eBayes(fit)
deg = topTable(fit,coef = 2,number = Inf)

#为deg数据框添加几列
#1.加probe_id列，把行名变成一列
library(dplyr)
deg = mutate(deg,probe_id = rownames(deg))
#2.加上探针注释
ids = distinct(ids,symbol,.keep_all = T)

#其他去重方式在zz.去重方式.R
deg = inner_join(deg,ids,by="probe_id")
nrow(deg) #如果行数为0就是你找的探针注释是错的

#3.加change列,标记上下调基因
logFC_t = 2
p_t = 0.005
#思考，如何使用padj而非p值
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)

## 画火山图
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,]
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)
)

富集分析

# 加ENTREZID列，用于富集分析（symbol转entrezid，然后inner_join）
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"))
# 多了几行少了几行都正常，SYMBOL与ENTREZID不是一对一的。
nrow(deg)
save(exp,Group,deg,logFC_t,p_t,file = "step4output.Rdata")

library(ggthemes)
library(org.Hs.eg.db)
library(dplyr)
library(enrichplot)

# 第一步：输入数据
gene_diff = deg$ENTREZID[deg$change != "stable"] 

# 第二步：富集
ekk <- enrichKEGG(gene = gene_diff,organism = 'hsa')
ekk <- setReadable(ekk,OrgDb = org.Hs.eg.db,keyType = "ENTREZID")
ego <- enrichGO(gene = gene_diff,OrgDb= org.Hs.eg.db,
                ont = "ALL",readable = TRUE)
# setReadable和readable = TRUE都是把富集结果表格里的基因名称转为symbol
class(ekk)

# 第三步：可视化
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
网上的资料和宝藏无穷无尽，学好R语言慢慢发掘~

羊&鹿
关注
  • 30
    点赞
  • 踩
  • 11
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
零基础入门转录组分析-数据处理（GEO数据库-芯片数据）配套资源
04-24
零基础入门转录组分析——数据处理（GEO数据库——芯片数据） https://blog.csdn.net/weixin_49878699/article/details/135436066 教程对应的原始文件+代码+处理好的数据
R软件--GEO数据分析教程：差异性分析、富集分析（GO\KEGG\DO）
weixin_40270911的博客
03-19 9770
title: “GSE93798” output: word_document knitr::opts_chunk$set(echo = TRUE) R Markdown 关于GSE93798的差异性分析、GO分析、KEGG分析、DO分析 rm(list = ls()) #RStudio控制台(console)中出现少量中文乱码的解决方法 Sys.setlocale("LC_ALL","...
零基础入门转录组数据分析——数据处理（GEO数据库——高通量测序数据）
weixin_49878699的博客
02-11 1万+
GEO数据库中高通量数据处理（结合了官方和自己理解），包括：基因symbol转化，获取count，fpkm处理，设置分组信息表。
GEO数据挖掘（一）基础介绍
热门推荐
doctor_yuxiang的博客
05-07 1万+
GEO数据挖掘，火山图，热图，主成分分析
GEO芯片数据基本分析
Annaaphq的博客
12-23 4493
GEO数据库芯片基本分析——R语言
零基础入门转录组数据分析——数据处理（GEO数据库——芯片数据）
weixin_49878699的博客
01-07 4254
GEO数据库全称GENE EXPRESSION OMNIBUS，是由美国国立生物技术信息中心NCBI创建并维护的基因表达数据库。它创建于2000年，收录了世界各国研究机构提交的高通量基因表达数据，也就是说只要是目前已经发表的论文，论文中涉及到的基因表达检测的数据都可以通过这个数据库中找到。并且这个网站作为各种高通量实验数据的公共存储库。这些数据包括基于单通道和双通道微阵列的实验，检测mRNA，基因组DNA和蛋白质丰度，以及非阵列技术，如基因表达系列分析（SAGE），质谱蛋白质组学数据和高通量测序数据。
学习笔记Day9&10：GEO数据挖掘-基因表达芯片代码
m0_71952366的博客
03-16 1165
不要在同一个文件夹中处理两套数据。
零基础入门转录组数据分析-DESeq2差异分析 配套资源
07-25
零基础入门转录组数据分析——DESeq2差异分析 https://blog.csdn.net/weixin_49878699/article/details/136114233 教程配套的原始数据+代码+处理好的数据文件
零基础入门转录组分析-数据处理（GEO数据库-高通量测序数据）配套资源
04-24
零基础入门转录组分析——数据处理（GEO数据库——高通量测序数据） https://blog.csdn.net/weixin_49878699/article/details/136032388 教程配套的原始数据+代码+处理好的数据文件 已验证可以从头到尾正常运行，可...
零基础入门转录组数据分析-GO+KEGG富集分析 配套资源
07-30
零基础入门转录组数据分析——GO+KEGG富集分析 https://blog.csdn.net/weixin_49878699/article/details/137052592 教程配套的原始数据+代码+处理好的数据文件
零基础入门转录组数据分析-机器学习算法之SVM-RFE（筛选特征基因） 配套资源
最新发布
07-31
零基础入门转录组数据分析——机器学习算法之SVM-RFE（筛选特征基因） https://blog.csdn.net/weixin_49878699/article/details/140748316 教程配套的原始数据+代码+处理好的数据文件
GEO数据分析举例GSE3335-GPL5175
weixin_50614492的博客
04-10 1806
GEO数据分析
GEO芯片数据处理、多芯片数据合并、差异分析、核心基因提取、核心基因与免疫细胞相关性
2302_80650915的博客
01-20 2707
下载探针矩阵probe.txt和平台注释ann.txt查看分组信息，复制对照组样本名称s1.txt和实验组样本名称s2.txt。
GEO官网下载series_matrix读取，GPL探针，GSE临床数据
hx2024的博客
10-08 1493
由于网络原因，常常使用getGEO函数可能会出错，下载不了文件。可以从GEO官网下载矩阵（部分可能没有矩阵）然后进行R读取进行ID转换分析等。
geo生存分析的简易做法
z2498249的博客
01-26 343
GEO简易生存分析
GEO数据挖掘-PCA、差异分析
weixin_57975238的博客
05-17 1267
From 生物技能树 GEO数据挖掘第二节。
【ONE·R || 两次作业（二）：GEO数据处理下载分析】
EPIPHANYECHO
01-09 2782
两次作业汇报·其二：GEO数据处理学习汇报。
OSG学习笔记-Geode（2-2）
liangfei868的博客
03-16 974
osg::Geode是OSG中的叶节点，它用于保存几何信息以便渲染。
基础公共数据库
GeekFocus
03-26 2665
生物医学数据库： 核酸，变异，表达，蛋白结构，功能，通路，表型，疾病，诊疗，药物 实验原始数据(公共数据)： 1KG：目标是发现人群中频率大于1%的变异位点。下载 GEO：GEO数据库全称GENE EXPRESSION OMNIBUS，是由美国国立生物技术信息中心NCBI创建并维护的基因表达数据库。它创建于2000年，收录了世界各国研究机构提交的高通量基因表达数据。 GEO芯片数据分析 BIG DDBJ EMBL-EBI ...
RNA-seq数据分析：生信实战指南
书中详细阐述了RNA-seq数据分析的整个流程，从原始数据的质控（Quality Control, QC），到数据预处理、比对（Alignment）、转录本组装（Transcriptome Assembly）、定量（Quantification）、差异表达分析...
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值