首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >GEO表达芯片数据分析

GEO表达芯片数据分析

原创
作者头像
小叮当aka
修改2023-03-23 15:58:04
2.7K1
修改2023-03-23 15:58:04
举报

title: "GEO表达芯片数据分析"

output: html_document

date: "2023-03-20"


关于该流程代码的说明:

(1)本流程仅适用于GEO芯片表达数据,以"GSE56649"为例

(2)先在GEO数据库中确定是否为"Expression profiling by array",不是的话不能使用本流程!

(3)注意需要自行修改或判断的代码一般放在了两个空行之间

(4)代码的注释有一丢丢多,目的是为了更好地帮助大家理解

1.下载数据,提取表达矩阵、临床信息和GPL编号

rm(list = ls())
library(GEOquery)

gse_number = "GSE56649" #自行修改编号!

eSet <- getGEO(gse_number, destdir = '.', getGPL = F)
#下载方法备选:geoChina函数,优点下载速度快,缺点较新数据未收入
# 2018年及以前的芯片数据可以用geoChina快速下载
#library(AnnoProbe)
#eSet <- geoChina(gse_number, destdir = '.') 
class(eSet)
#> [1] "list"
length(eSet)
#> [1] 1
eSet = eSet[[1]] #当一个数据集有多个平台时,分别提取子集,当成多个数据分析

#(1)提取表达矩阵信息exp
exp <- exprs(eSet)
dim(exp)
#> [1] 54675    22
exp[1:4,1:4]
#>           GSM1366348 GSM1366349 GSM1366350 GSM1366351
#> 1007_s_at    279.156    202.866    406.818    232.668
#> 1053_at      487.700    409.018    393.810    397.127
#> 117_at       666.869    388.589    704.633    953.481
#> 121_at       240.646    361.198    305.229    374.757
# 目的:
#1)检查矩阵是否正常:
#如果是空的就会报错,空的和有异常值的矩阵需要处理原始数据。
#如果表达矩阵为空,大多数是转录组数据,不能用这个流程(后面另讲)
#2)判断是否需要取log:
#几百、几千的话就需要取log,一般0-20的说明取过log了,不用再取
#没取过log,有负值(光信号值没负数)→错误数据(弃用或处理原始数据)
#取过log,有负值(数值偏小时)→正常

exp = log2(exp+1)  #自行判断,不需要取log时注释掉

boxplot(exp)
# 目的:
#1)再次判断是否需要取log:
#箱线图如果表现为压得很扁,说明没取log
#2)检查并处理异常样本:
#如果有样本的表达量整体低于其他的一大截,运行如下代码去除异常
#exp = limma::normalizeBetweenArrays(exp)
#boxplot(exp)
#3)表达矩阵的数字一半正一半负:做了标准化,弃用或处理原始数据

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

#(3)让exp列名与pd的行名顺序完全一致
p = identical(rownames(pd),colnames(exp));p
#> [1] TRUE
if(!p) exp = exp[,match(rownames(pd),colnames(exp))]
# 分组信息来自临床信息,分组信息需要与表达矩阵一一对应
# 因此临床信息的行需要与表达矩阵的列一一对应

#(4)提取芯片平台编号
gpl_number <- eSet@annotation;gpl_number
#> [1] "GPL570"
save(gse_number,pd,exp,gpl_number,file = "step1output.Rdata")

2.获取实验分组(Group)和探针注释(ids)

rm(list = ls())  
load(file = "step1output.Rdata")
library(stringr)
#(1)获取实验分组(Group)-----------------
# 标准流程代码是二分组,多分组数据的分析不适用,后面另讲
# 生成Group向量有三种常规方法,根据数据三选一即可
# 选谁就把第几个逻辑值写成T,另外两个为F
# 如果三种办法都不适用,可以继续往后写else if

if(F){
  # 第一种方法:有现成的可以用来分组的列
  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")
}
# 只运行了第三种方法,因为只有第三种方法的if后面是T

Group = factor(Group,levels = c("control","RA")) #修改level!

Group
#>  [1] RA      RA      RA      RA      RA      RA      RA     
#>  [8] RA      RA      RA      RA      RA      RA      control
#> [15] control control control control control control control
#> [22] control
#> Levels: control RA
# Group是一个分类型向量,分类型数据可以用因子factor组织
# 需要把Group转换成因子,并设置参考水平,指定levels
# 注意levels顺序要设置成对照组在前,处理组在后
# 因子正文与levels不对应时,会产生NA

#(2)获取探针注释(ids)-----------------
# 探针注释:探针与基因的对应关系,probe_id和symbol
# 捷径注释的方法(首选,能找到70%-80%,失败后选方法1,以此类推)
library(tinyarray)
find_anno(gpl_number) #打出找注释的代码
#> [1] "`library(hgu133plus2.db);ids <- toTable(hgu133plus2SYMBOL)` and `ids <- AnnoProbe::idmap('GPL570')` are both avaliable"

ids <- AnnoProbe::idmap('GPL570')  #修改芯片平台编号!

head(ids) #调试
#>         probe_id symbol
#> 193731   1053_at   RFC2
#> 193732    117_at  HSPA6
#> 193733    121_at   PAX8
#> 193734 1255_g_at GUCA1A
#> 193735   1316_at   THRA
#> 193736   1320_at PTPN21
#四种探针注释的方法,先看捷径,方法1里找不到就从方法2找,以此类推-----------------
#方法1 Bioconductor R包(最常用,但不是万能的!)
if(F){ 
  #要用方法1时才需要改成T
gpl_number 
#http://www.bio-info-trainee.com/1399.html
if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db")
library(hgu133plus2.db)
ls("package:hgu133plus2.db")  #看R包有哪些数据,看SYMBOL的出处
ids <- toTable(hgu133plus2SYMBOL) #提取R包的注释结果
head(ids)
}
# 方法2 读取GPL网页的表格文件,按列取子集
##https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570
if(F){ 
  #要用方法2时才需要改成T
  #注:表格读取参数、文件列名不统一,活学活用,有的表格里没有symbol列,也有的GPL平台没有提供注释表格
  b = read.delim("GPL570-55999.txt",
                 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
}
# 方法3 官网下载注释文件并读取(看运气,不一定能找到)
##http://www.affymetrix.com/support/technical/byproduct.affx?product=hg-u133-plus
# 方法4 自主注释 
#https://mp.weixin.qq.com/s/mrtjpN8yDKUdCSvSUuUwcA

# 保存数据-----------------
save(exp,Group,ids,gse_number,file = "step2output.Rdata")

3.绘制PCA图和热图

rm(list = ls())  
#输入数据:exp和Group
load(file = "step1output.Rdata")
load(file = "step2output.Rdata")

#(1)PCA 图(样本聚类图)----
# Principal Component Analysis
# 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")

#(2)top 1000 sd 热图---- 
cg=names(tail(sort(apply(exp,1,sd)),1000))
n=exp[cg,] #筛选标准差最大的1000个基因(最活跃、离散的基因)

# 直接画热图,对比不鲜明
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
)
# 按行标准化(scale = "row")
#比较的是一行即一个基因在不同分组间的差别
#放大行内部的差别,抹平行与行之间的差别,图例范围缩小
#使基因在不同分组的表达差异更明显
pheatmap(n,
         show_colnames =F,
         show_rownames = F,
         annotation_col=annotation_col,
         scale = "row",
         breaks = seq(-3,3,length.out = 100) #设置色带范围,避免极端值影响整体配色
         ) 
dev.off()
#> null device 
#>           1

4.差异分析(用limma包来做)

rm(list = ls()) 
load(file = "step2output.Rdata")
#需要表达矩阵exp和分组信息Group,不需要改
library(limma)
design=model.matrix(~Group)           #构建一个模型矩阵
fit=lmFit(exp,design)                 #线性拟合
fit=eBayes(fit)                       #贝叶斯检验
deg=topTable(fit,coef=2,number = Inf) #提取差异分析结果
#可以将上面几句代码运用简化思维组织成一个函数function(exp,Group)

#为deg数据框添加几列
#(1)加probe_id列,把行名变成一列----
library(dplyr)
deg <- mutate(deg,probe_id=rownames(deg))

#(2)加上探针注释symbol列,随机去重更客观----
# 随机去重,因为存在多个探针对应一个symbol
ids = ids[!duplicated(ids$symbol),]      #方法1
ids = distinct(ids,symbol,.keep_all = T) #方法2
# 随机去重更客观,其他去重方式----
if(F){
  #保留最大值去重
  exp2 = exp[ids$probe_id,]
  identical(ids$probe_id,rownames(exp2))
  ids = ids[order(rowSums(exp2),decreasing = T),]
  ids = ids[!duplicated(ids$symbol),];nrow(ids)
  # 拿这个ids去inner_join
}
if(F){
  exp3 = exp[ids$probe_id,]
  rownames(exp3) = ids$symbol
  exp3[1:4,1:4]
  exp4 = limma::avereps(exp3)
  # 此时拿到的exp4已经是一个基因为行名的表达矩阵
  # 直接差异分析,不再需要inner_join 
}
# 用inner_join加上symbol列----
deg <- inner_join(deg,ids,by="probe_id")
nrow(deg)
#> [1] 20188

#(3)加change列,标记上下调基因----
#logFC_t和p_t可自行调整
logFC_t=1 #差异基因过多,可调大(2.2,2,1.5,1.2,0.585=log2(1.5))
p_t = 0.05#P值越小,差异越显著,可调小(0.01)
#P.Value可调成adj.P.Val的列名
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 
#>    796  18753    639

#(4)加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)#转换依据的数据库(人类)
# symbol和entrezid并非一一对应,损失或增加一些很正常
# 其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
save(Group,deg,logFC_t,p_t,gse_number,file = "step4output.Rdata")

5.可视化——火山图、差异基因热图、相关性热图

rm(list = ls()) 
load(file = "step1output.Rdata")
load(file = "step4output.Rdata")
#(1)火山图----
library(dplyr)
library(ggplot2)
dat  = distinct(deg,symbol,.keep_all = T) #按symbol列去重复
#去掉添加ENTREZID后多的行
#防止需要添加到图中的symbol重复出现

#1)作火山图
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_t),lty=4,col="black",linewidth=0.8) +
  theme_bw() #去掉灰色背景
p
#2)添加火山图的标签:
#① 自己添加感兴趣的基因
#② P.Value的top10
#③ 根据logFC挑也行

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(T){
  #取前10上调和前10下调,可以根据logFC(更合适),也可以根据Pvalue
  library(dplyr)
  dat2 = dat %>%
    filter(change!="stable") %>% 
    arrange(logFC)              #将logFC从小到大排序
  cg = c(head(dat2$symbol,10),  #logFC最小的10个基因
         tail(dat2$symbol,10))  #logFC最大的10个基因
}else{ #没有运行
  #全部差异基因
  cg = dat$symbol[dat$change !="stable"]
  length(cg)
}
n=exp[cg,]
dim(n)
#> [1] 20 22
#绘制差异基因热图
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 
#$gtable提取热图的主体部分,才能完成拼图的操作

#(3)感兴趣基因的相关性热图----
# 注意:这里的流程只适用于多个基因的相关性
# 两个基因的相关性需要用相关性点图表示
library(corrplot)
g = sample(deg$symbol[1:500],10) 
# 这里是随机取样,注意换成自己感兴趣的基因
g
#>  [1] "YWHAB"   "P3H2"    "COIL"    "ARL6IP5" "MYDGF"   "R3HDM4" 
#>  [7] "ALYREF"  "RIF1"    "PBX2"    "TDG"
M = cor(t(exp[g,]))  #相关系数矩阵
#计算的是列与列的相关性,需要转置才能计算基因相关性
pheatmap(M)        
library(paletteer) #配色R包
my_color = rev(paletteer_d("RColorBrewer::RdYlBu"))
#选出一组颜色,通过逆转配色,使红色代表正相关,蓝色代表负相关
my_color = colorRampPalette(my_color)(10) #切10种颜色
corrplot(M, type="upper",
         method="pie",
         order="hclust", 
         col=my_color,
         tl.col="black", 
         tl.srt=45)
library(cowplot)         #cowplot可以拼不同类型的图
cor_plot <- recordPlot() #把画板上的图记录下来

# 拼图
load("pca_plot.Rdata")
plot_grid(pca_plot,cor_plot,
          volcano_plot,heatmap_plot$gtable)
          #右下角画板太小时会报错
dev.off()
#> null device 
#>           1
#保存
pdf("deg.pdf", width = 10, height = 10)
#width,height可以调整导出图片的比例
plot_grid(pca_plot,cor_plot,
          volcano_plot,heatmap_plot$gtable)
dev.off()
#> null device 
#>           1

6.检验差异分析的分组是否颠倒

# 通过判断某个基因的箱线图和logFC是否匹配
check_gene = deg$symbol[1]
boxplot(exp[check_gene,]~Group) #处理组<对照组→logFC<0
deg$logFC[1]                    #logFC<0成立
#> [1] -1.633518

7.富集分析

# 富集分析:衡量每个通路里的基因在差异基因里是否足够多
rm(list = ls())  
load(file = 'step4output.Rdata')
library(clusterProfiler)
library(ggthemes)
library(org.Hs.eg.db)
library(dplyr)
library(ggplot2)
library(stringr)
library(enrichplot)

# 1.GO 富集分析----
#(1)输入数据
gene_up = deg$ENTREZID[deg$change == 'up'] 
gene_down = deg$ENTREZID[deg$change == 'down'] 
gene_diff = c(gene_up,gene_down) #差异基因(上下调)

#(2)富集
#以下限速步骤耗时很长,设置了存在即跳过
#判断一个文件在工作目录下是否存在,存在就不运行
f = paste0(gse_number,"_GO.Rdata")
if(!file.exists(f)){
  ego <- enrichGO(gene = gene_diff,
                  OrgDb= org.Hs.eg.db,
                  ont = "ALL",
                  readable = TRUE) #直接转换为symbol
  ego_BP <- enrichGO(gene = gene_diff,
                  OrgDb= org.Hs.eg.db,
                  ont = "BP",
                  readable = TRUE)
  #ont参数:One of "BP", "MF", and "CC" subontologies, or "ALL" for all three.
  save(ego,ego_BP,file = f)
}
load(f)

#(3)可视化
#1)条带图
barplot(ego)
barplot(ego, split = "ONTOLOGY", font.size = 10, 
        showCategory = 5) + 
  facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y") 
#分面,根据自己的选择决定是否运行

#2)气泡图
dotplot(ego)
dotplot(ego, split = "ONTOLOGY", font.size = 10, 
        showCategory = 5) + 
  facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y") 
#3)展示top通路的共同基因,要放大看。
#gl 用于设置下图的颜色
gl = deg$logFC
names(gl)=deg$ENTREZID
#Gene-Concept Network,要放大看
cnetplot(ego,
         # layout = "star",
         color.params = list(foldChange = gl),
         showCategory = 3)
#4)展示GO term 间的关系
goplot(ego_BP) #只能画一个子类
# https://zhuanlan.zhihu.com/p/99789859
# 有颜色的点是富集到的通路
# 灰色的点是根据有颜色的点挖掘到的通路

# 2.KEGG pathway analysis----
#上调、下调、差异、所有基因
#(1)输入数据(与GO完全一致)
gene_up = deg[deg$change == 'up','ENTREZID'] 
gene_down = deg[deg$change == 'down','ENTREZID'] 
gene_diff = c(gene_up,gene_down)

#(2)对上调/下调/所有差异基因进行富集分析
f2 = paste0(gse_number,"_KEGG.Rdata")
if(!file.exists(f2)){
  kk.up <- enrichKEGG(gene           = gene_up,
                      organism       = 'hsa')
  kk.down <- enrichKEGG(gene         =  gene_down,
                        organism     = 'hsa')
  kk.diff <- enrichKEGG(gene         = gene_diff,
                        organism     = 'hsa')
  save(kk.diff,kk.down,kk.up,file = f2)
}
load(f2)

#(3)看看富集到了吗?https://mp.weixin.qq.com/s/NglawJgVgrMJ0QfD-YRBQg
table(kk.diff@result$p.adjust<0.05)
#> 
#> FALSE 
#>   317
table(kk.up@result$p.adjust<0.05)
#> 
#> FALSE 
#>   292
table(kk.down@result$p.adjust<0.05)
#> 
#> FALSE 
#>   285
# 富集不到的话,很正常,但有补救方法:
# ①p.adjust改为P.Value
# ②换一种富集的方式GSEA
# ③调整差异基因的阈值

#(4)双向富集图
# 富集分析所有图表默认都是用p.adjust,富集不到可以退而求其次用p值,在文中说明即可
# 调用kegg_plot函数
library(ggthemes)
library(ggplot2)
kegg_plot <- function(kk.up,kk.down,p_t = 0.05,p.adjust = F){
  down_kegg <- kk.down@result %>%
    filter(pvalue<p_t) %>% #筛选行
    mutate(group=-1) #新增列
  
  up_kegg <- kk.up@result %>%
    filter(pvalue<p_t) %>%
    mutate(group=1)
  dat=rbind(up_kegg,down_kegg)
  if(p.adjust) dat$pvalue = dat$p.adjust
  colnames(dat)
  dat$pvalue = -log10(dat$pvalue)
  dat$pvalue=dat$pvalue*dat$group 
  
  dat=dat[order(dat$pvalue,decreasing = F),]
  
  gk_plot <- ggplot(dat,aes(reorder(Description, pvalue), y=pvalue)) +
    geom_bar(aes(fill=factor((pvalue>0)+1)),stat="identity", width=0.7, position=position_dodge(0.7)) +
    coord_flip() +
    scale_fill_manual(values=c("#0072B2", "#B20072"), guide="none") +
    labs(x="", y="" ) +
    theme_pander()  +
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          #axis.ticks.x = element_blank(),
          axis.line.x = element_line(linewidth = 0.3, colour = "black"),#x轴连线
          axis.ticks.length.x = unit(-0.20, "cm"),#修改x轴刻度的高度,负号表示向上
          axis.text.x = element_text(margin = margin(t = 0.3, unit = "cm")),##线与数字不要重叠
          axis.ticks.x = element_line(colour = "black",linewidth = 0.3) ,#修改x轴刻度的线                         
          axis.ticks.y = element_blank(),
          axis.text.y  = element_text(hjust=0),
          panel.background = element_rect(fill=NULL, colour = 'white')
    )
}
# 作图:
g_kegg <- kegg_plot(kk.up,kk.down,p_t = 0.01,p.adjust = F)
g_kegg
# 注意:
#①p.adjust可改为T,按照p.adjust作图
#②可以在括号里自己用tab加参数
#③根据不同的横坐标改变,因为理论上横坐标不可能为负数
#g_kegg +scale_y_continuous(labels = c(2,0,2,4,6)) 

引用自生信技能树

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

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

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

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

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 关于该流程代码的说明:
    • (1)本流程仅适用于GEO芯片表达数据,以"GSE56649"为例
      • (2)先在GEO数据库中确定是否为"Expression profiling by array",不是的话不能使用本流程!
        • (3)注意需要自行修改或判断的代码一般放在了两个空行之间
          • (4)代码的注释有一丢丢多,目的是为了更好地帮助大家理解
          • 1.下载数据,提取表达矩阵、临床信息和GPL编号
          • 2.获取实验分组(Group)和探针注释(ids)
          • 3.绘制PCA图和热图
          • 4.差异分析(用limma包来做)
          • 5.可视化——火山图、差异基因热图、相关性热图
          • 6.检验差异分析的分组是否颠倒
          • 7.富集分析
          领券
          问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档