GEO数据挖掘

原创

可乐同学与生信死磕到底

修改于 2024-04-09 21:45:20

780

修改于 2024-04-09 21:45:20

文章被收录于专栏：今天我入门生信了吗？今天我入门生信了吗？

1 图表介绍

1.1 热图

输入数据：数值型矩阵/数据框

颜色深浅代表数值的大小

1.2 散点图

1.3 箱线图

1.3.1 输入数据

横坐标：一个有重复值的离散型变量

纵坐标：连续型向量

1.3.2 箱线图中五条线的含义

箱线图比较分布情况

箱型图不显示原始数据点，而是采用样本数据，根据四分位数用盒和线来显示值的范围。此外，它们用星号显示落在箱须之外的离群值

箱形图显示五个数据：

1、最小的数字（最小值）

2、第一个四分位数（25%位点值）

3、中间的数字（中位值）

4、第三个四分位数（75%位点值）

5、最大的数字（最大值）

箱线图用于比较单个基因在两组之间（control/treat）之间的表达量差异

在多基因中用于选出分布差异较大的基因

1.4 火山图

1.4.1 火山图的横纵坐标及其含义

1.4.1.1 横坐标：logFC

Foldchange（FC）：处理组平均值/对照组平均值

logFoldchange（logFC）：Foldchange取log2

表达矩阵中的count一般为取过log之后的数值

处理组在前，对照组在后！反了全错!

1.4.1.2 纵坐标：p值

1.4.2 火山图怎么看?

logFC>0,treat>control,基因表达量上升；

logFC<0,treat<control,基因表达量下降。

基因表达上升/下降≠差异基因！还要看阈值的设定（logFC+p值）

通常所说的上调/下调基因：表达量显著上升/下降的基因

上调基因:logFC>1，p<0.01

下调基因:logFC<-1，p<0.01

可以通过“阈值调整”改变差异基因的数量

logFC的常见阈值：1、2、1.2、1.5、2.2

火山图是用limma差异分析结果来做的，limma差异分析结果是一个10列的数据框

1.5 PCA主成分分析图

1.5.1 PCA的原理

主成分分析：旨在利用降维的思想，把多指标转化为少数几个综合指标（即主成分）

根据这些主成分对样本进行聚类，代表样本的点（中心点除外）在坐标轴上的距离越远，说明样本差异越大

1.5.2 PCA的用途

用于“预实验”，简单查看组间是否有差别

同一分组是否聚成一簇（组内重复好）
中心点之间是否有距离（组间差别大）

从这里开始没有课件，以下内容为自己结合课堂视频整理得出~

2 GEO背景知识+表达芯片分析思路

2.1 表达数据实验设计

实验目的：通过基因表达量数据的差异分析和富集分析来解释生物学现象

有差异的材料→差异基因→找功能/找关联→解释差异，缩小基因范

?数据挖掘本质：逐步缩小感兴趣的基因范围！

2.2 GEO数据库介绍

GSM：用户提交给GEO的样本数据（Sample）

GSE：一个完整的研究，提供了整个研究的描述（Series）

GPL：用户测定表达量使用的芯片/平台（Platform）

2.3 基因表达芯片的原理

探针的表达量代表基因的表达量

2.4 分析思路

2.5 表达矩阵

探针id要找到对应的基因sample
样本编号GSM要获取分组信息group

2.6 富集分析

2.6.1 什么是基因的Entrezid？

输入数据：差异基因的Entrezid

Symbol为常说的基因名

并非一一对应，会损失/增加一部分基因

2.6.2 富集分析的数据库

2.6.2.1 KEGG数据库

通路pathway

2.6.2.2 GO数据库

细胞组分（cellular component）

分子功能（molecular function）

生物过程（biological process）

2.6.3 富集分析结果

GeneRatio：差异基因中有哪些属于这条通路/差异基因中有多少个被数据库收录

BgRatio：这条通路一共有多少个基因/数据库中所有通路总共多少个基因

?富集分析本质：衡量某个通路里的基因在差异基因里是否足够多

2.6.4 富集分析的可视化

气泡图、柱状图/条形图

Y叔Clusterprofiler 默认使用p.adjust

可以按照CC、MF、BP对图片进行分面

也可以上、下调基因分开富集，合并画图

3 代码分析流程

3.1 安装需要的R包

options("repos"="https://mirrors.ustc.edu.cn/CRAN/")
if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")

cran_packages <- c('tidyr',
                   'tibble',
                   'dplyr',
                   'stringr',
                   'ggplot2',
                   'ggpubr',
                   'factoextra',
                   'FactoMineR',
                   'devtools',
                   'cowplot',
                   'patchwork',
                   'basetheme',
                   'paletteer',
                   'AnnoProbe',
                   'ggthemes',
                   'VennDiagram',
                   'tinyarray') 
Biocductor_packages <- c('GEOquery',
                         'hgu133plus2.db',
                         'ggnewscale',
                         "limma",
                         "impute",
                         "GSEABase",
                         "GSVA",
                         "clusterProfiler",
                         "org.Hs.eg.db",
                         "preprocessCore",
                         "enrichplot")

for (pkg in cran_packages){
  if (! require(pkg,character.only=T) ) {
    install.packages(pkg,ask = F,update = F)
    require(pkg,character.only=T) 
  }
}


for (pkg in Biocductor_packages){
  if (! require(pkg,character.only=T) ) {
    BiocManager::install(pkg,ask = F,update = F)
    require(pkg,character.only=T) 
  }
}

#前面的所有提示和报错都先不要管。主要看这里
for (pkg in c(Biocductor_packages,cran_packages)){
  require(pkg,character.only=T) 
}

3.2 GEO数据下载并从中提取有用信息

rm(list = ls())
library(GEOquery)
#先去网页确定是否是表达芯片数据，不是的话不能用本流程。
gse_number = "GSE56649"
eSet <- getGEO(gse_number, destdir = '.', getGPL = F) #实现下载并读取

eSet = eSet[[1]] #eSet脱离列表的壳子

R语言中狭义的对象：R包的作者以某种特定的方式组织起来的数据

ExpressionSet对象出自Biobase包

#(1)提取表达矩阵exp
exp <- exprs(eSet)
dim(exp)
exp[1:4,1:4]
#检查矩阵是否正常，如果是空的就会报错，空的和有负值的、有异常值的矩阵需要处理原始数据。
#如果表达矩阵为空，大多数是转录组数据，不能用这个流程(后面另讲)。
#自行判断是否需要log
exp = log2(exp+1)
boxplot(exp)

取过log的数据正常范围在0-20之间

画箱线图看有没有异常数据

#(2)提取临床信息
pd <- pData(eSet)

#(3)让exp列名与pd的行名顺序完全一致 临床信息中的分组信息与表达矩对应
p = identical(rownames(pd),colnames(exp));p
if(!p) exp = exp[,match(rownames(pd),colnames(exp))]

#(4)提取芯片平台编号
gpl_number <- eSet@annotation;gpl_number #@给对象提取子集
save(gse_number,pd,exp,gpl_number,file = "step1output.Rdata")

3.3 ?获取分组信息和探针注释

长脚本管理方式：

2个脚本之间的衔接：清空环境变量+load Rdata

3.3.1 获取分组信息的三种方法：

有现成的可以用来分组的列
自己生成
使用字符串处理的函数获取分组

# Group(实验分组)和ids(探针注释)
rm(list = ls())  
load(file = "step1output.Rdata")
library(stringr)
# 标准流程代码是二分组，多分组数据的分析后面另讲
# 生成Group向量的三种常规方法，三选一，选谁就把第几个逻辑值写成T，另外两个为F。如果三种办法都不适用，可以继续往后写else if
if(F){
  # 1.Group----
  # 第一种方法，有现成的可以用来分组的列
  Group = pd$`disease state:ch1` 
}else if(F){
  # 第二种方法，自己生成
  Group = c(rep("RA",times=13),
            rep("control",times=9))
  Group = rep(c("RA","control"),times = c(13,9))
}else if(T){
  # 第三种方法，使用字符串处理的函数获取分组
  Group=ifelse(str_detect(pd$source_name_ch1,"control"),
               "control",
               "RA")}

因子：带有levels属性的特殊的向量

# 需要把Group转换成因子，并设置参考水平，指定levels，对照组在前，处理组在后
Group = factor(Group,levels = c("control","RA"))
Group

在第一个位置上的为参考水平，做差异分析时作为对照组

3.3.2 探针注释的获取

3.3.2.1 探针注释的定义及来源

探针注释：探针与基因的对应关系不是所有的GPL都能找到注释！

注释来源：

Bioconductor的注释包
GPL的表格文件解析
官网下载对应产品的注释表格
自主注释3.3.2.2 探针注释的代码

library(tinyarray)
find_anno(gpl_number) #打出找注释的代码
ids <- AnnoProbe::idmap('GPL570') #此时已经找到了探针注释，后面的代码不需要再运行

找探针注释的四种方法：（原始、基础）

Bioconductor中的R包（最常用）在> http://www.bio-info-trainee.com/1399.html 网站上找

gpl_number 
#http://www.bio-info-trainee.com/1399.html
if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db") #加后缀.db
library(hgu133plus2.db)
ls("package:hgu133plus2.db") #列出一个R包里有什么函数/数据
ids <- toTable(hgu133plus2SYMBOL)
head(ids)

读取GPL网页的表格文件，按列取子集

##https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570
if(F){
  #注：表格读取参数、文件列名不统一，活学活用，有的表格里没有symbol列，也有的GPL平台没有提供注释表格
  b = read.delim("GPL570-55999.txt", #read.table()的小分支
                 check.names = F, #不要把列名里的特殊字符转化为.
                 comment.char = "#") #开头的行作为注释行
  colnames(b)
  ids2 = b[,c("ID","Gene Symbol")]
  colnames(ids2) = c("probe_id","symbol") #为了后续的两个表格连接改列名
  k1 = ids2$symbol!="";table(k1) #去除没有基因对应的探针
  k2 = !str_detect(ids2$symbol,"///");table(k2) #去除非特异性探针（一个探针对应多个基因）
  ids2 = ids2[ k1 & k2,]
  # ids = ids2
  }

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

补充：AnnoProbe包

AnnoProbe包里的三个函数：

idmap：下载注释信息（探针与基因symbol的对应关系）获得ids

geoChina：下载GSE数据

annoGene：给基因提供注释

library(AnnoProbe)
?idmap
ids=idmap('GPL570')
?geoChina
geoChina(gse = "GSE2546", mirror = "tercent")
?annoGene

tinyarray包：geo_download函数

get_deg_all直接进行差异分析及可视化 需要找到group和ids

3.4 画PCA图+Top1000基因热图

3.4.1 PCA图

输入数据：表达矩阵+分组信息

PCA示例来源：http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials

dat=as.data.frame(t(exp)) #转置 转置以后都变成了矩阵 再从矩阵变为数据框
library(FactoMineR)
library(factoextra) 
dat.pca <- PCA(dat, graph = FALSE) 
pca_plot <- 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_plot
save(pca_plot,file = "pca_plot.Rdata") #方便后续与其他图进行拼图

画图思路：

3.4.2 热图

cg=names(tail(sort(apply(exp,1,sd)),1000)) #隐式循环取方差最大的1000个基因
n=exp[cg,] 

# 直接画热图，对比不鲜明 直接拿表达量数据画图
library(pheatmap)
annotation_col=data.frame(group=Group)
rownames(annotation_col)=colnames(n) 
pheatmap(n,
         show_colnames =F,
         show_rownames = F,
         annotation_col=annotation_col
)

# 按行标准化 让热图对比更加鲜明 突出显示行内部的变化
pheatmap(n,
         show_colnames =F,
         show_rownames = F,
         annotation_col=annotation_col, #在上方添加注释条
         scale = "row",
         breaks = seq(-3,3,length.out = 100) #设置颜色分配范围→避免离群值对整张图产生的影响
         ) 
dev.off()

3.5 差异分析

差异分析得到的结果：

rm(list = ls()) 
load(file = "step2output.Rdata")
#差异分析，用limma包来做
#需要表达矩阵和Group，不需要改
library(limma) #输入数据Group和exp，输出数据deg 把以下几行代码看作一个整体
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 = ids[!duplicated(ids$symbol),]
#其他去重方式在zz.去重方式.R
deg <- inner_join(deg,ids,by="probe_id")
nrow(deg)

#3.加change列,标记上下调基因
logFC_t=1
P.Value_t = 0.05
k1 = (deg$P.Value < P.Value_t)&(deg$logFC < -logFC_t)
k2 = (deg$P.Value < P.Value_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) #转换依据 Hs为人类
s2e <- bitr(deg$symbol, 
            fromType = "SYMBOL", #从……基因类型转换为……基因类型
            toType = "ENTREZID",
            OrgDb = org.Hs.eg.db)#人类
#其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
save(Group,deg,logFC_t,P.Value_t,gse_number,file = "step4output.Rdata")

探针注释：多个探针对应一个基因

随机去重
保留行和/行平均值最大的探针
取多个探针的平均值

如何实现随机去重？

3.6 画图

rm(list = ls()) 
load(file = "step1output.Rdata")
load(file = "step4output.Rdata")
#1.火山图----
library(dplyr)
library(ggplot2)
dat  = deg[!duplicated(deg$symbol),]

p <- ggplot(data = dat, 
            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.Value_t),lty=4,col="black",linewidth=0.8) +
  theme_bw()
p
#在火山图上添加感兴趣的基因 叠加一个新的图层
#ggplot叠加的不同图层可以使用不同的数据
for_label <- dat%>% 
  filter(symbol %in% c("HADHA","LRRFIP1"))

volcano_plot <- p +
  geom_point(size = 3, shape = 1, data = for_label) +
  ggrepel::geom_label_repel(
    aes(label = symbol),
    data = for_label,
    color="black"
  )
volcano_plot

#2.差异基因热图----

load(file = 'step2output.Rdata')
# 表达矩阵行名替换
exp = exp[dat$probe_id,]
rownames(exp) = dat$symbol
if(F){
  #全部差异基因
  cg = dat$symbol[dat$change !="stable"]
  length(cg)
}else{
  #取前10上调和前10下调
  library(dplyr)
  dat2 = dat %>%
    filter(change!="stable") %>% 
    arrange(logFC) 
  cg = c(head(dat2$symbol,10),
         tail(dat2$symbol,10))
}
n=exp[cg,]
dim(n)

#差异基因热图
library(pheatmap)
annotation_col=data.frame(group=Group)
rownames(annotation_col)=colnames(n) 
heatmap_plot <- pheatmap(n,show_colnames =F,
                         scale = "row",
                         #cluster_cols = F, 
                         annotation_col=annotation_col,
                         breaks = seq(-3,3,length.out = 100)
) 
heatmap_plot
#拼图
library(patchwork)
volcano_plot +heatmap_plot$gtable

# 3.感兴趣基因的相关性----
library(corrplot)
g = deg$symbol[1:10] # 换成自己感兴趣的基因
g

M = cor(t(exp[g,]))
pheatmap(M)

library(paletteer)
my_color = rev(paletteer_d("RColorBrewer::RdYlBu"))
my_color = colorRampPalette(my_color)(10)
corrplot(M, type="upper",
         method="pie",
         order="hclust", 
         col=my_color,
         tl.col="black", 
         tl.srt=45)
library(cowplot)
cor_plot <- recordPlot() 

# 拼图
load("pca_plot.Rdata")
plot_grid(pca_plot,cor_plot,
          volcano_plot,heatmap_plot$gtable)
dev.off()
#保存
pdf("deg.pdf")
plot_grid(pca_plot,cor_plot,
          volcano_plot,heatmap_plot$gtable)
dev.off()

以上内容为自己学习生信技能树课程的学习笔记，GEO数据挖掘这部分内容实在是太多啦！笔记字数破万啦！兜兜转转写了一个星期才写完，写完回过头看才觉得应该把理论知识和代码分开来记录，后面的复杂数据分析部分就放在下一篇笔记吧~

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

学习笔记

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

学习笔记

登录后参与评论

0 条评论

热度