前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布

GEO

原创
作者头像
浅念
发布2023-04-04 16:17:12
1.5K0
发布2023-04-04 16:17:12
举报
文章被收录于专栏:syj生信

一、数据下载

1.加载我们将要使用的R包

代码语言:javascript
复制
library(GEOquery)

2.先去网页确定是否是表达芯片数据,不是的话不能用本流程。

代码语言:txt
复制
gse_number = "GSE56649"
eSet <- getGEO(gse_number, destdir = '.', getGPL = F)#getGEO有从GEO中下载数据到工作目录下,并将数据读取到R中。

#若网络不好,下载数据慢,可以换一个R包

代码语言:txt
复制
library(AnnoProbe)
gse_number = "GSE56649"
eSet=geoChina(gse_number)

3.查看数据

代码语言:txt
复制
 class(eSet)#列表
 length(eSet)#列表的长度
 eSet = eSet[[1]]
 exp <- exprs(eSet)#(1)提取表达矩阵exp
 dim(exp)#矩阵几行几列
 exp[1:4,1:4]#看数据是否正常,看数据是否取过log,如果取过log,则数据在0~20中间差不多。还可以画箱线图看是否正常,正确情况下每个样本的整体都差不多。

#若数据没有取log

代码语言:txt
复制
exp = log2(exp+1)#之所以要+1是因为害怕exp有数据=0,这样log2(0)就是负无穷了。
boxplot(exp)

二、提取临床信息

代码语言:txt
复制
 pd <- pData(eSet)

三、让exp列名与pd的行名顺序完全一致

分组信息的每一列与表达矩阵的每一行是对应关系

代码语言:txt
复制
p = identical(rownames(pd),colnames(exp));p#判断两个数据的行名和列名是否一致
if(!p) exp = exp[,match(rownames(pd),colnames(exp))]#如果不一致,进行该行操作。

四、提取芯片平台编号

代码语言:txt
复制
gpl_number <- eSet@annotation;gpl_number#提取的符号有@和$,该用哪个,我们可以一一试一试。
save(gse_number,pd,exp,gpl_number,file = "step1output.Rdata")#保存编号,ges_number是保存总的样本集编号,pd是临床信息,exp是表达矩阵、gpl_number是芯片平台编号

五、Group(实验分组)和ids(探针注释)

代码语言:javascript
复制
load(file = "step1output.Rdata")
library(stringr)

标准流程代码是二分组,多分组数据的分析后面另讲

生成Group向量的三种常规方法,三选一,选谁就把第几个逻辑值写成T,另外两个为F。如果三种办法都不适用,可以继续往后写else if

代码语言:txt
复制
 if(F){
  # 1.Group----
  # 第一种方法,有现成的可以用来分组的列
  Group = pd$`disease state:ch1` #这个`是在我们提取的一列名称中有空格时,R语言会自动生成`
}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$title,"control"),#str_detece是检查是否有这个字符串的意思,如果有是T,如果没有是F
               "control",#ifelse是如果括号中的检测结果是T,则用control替换T,否则用RA替换F
               "RA")
}

 

需要把Group转换成因子,因子相比group里的字符串少了双引号,并设置参考水平,指定levels,对照组在前,处理组在后

代码语言:txt
复制
Group = factor(Group,levels = c("control","RA"))# 如果不加level,那么将按因子名字的首字母排序
Group

六、探针注释的获取

捷径

代码语言:txt
复制
library(tinyarray)
find_anno(gpl_number) # 打出找注释的代码,根据gpl_number来找探针序号所对应的真正的探针基因序列
ids <- AnnoProbe::idmap('GPL570')

四种方法:

方法1 BioconductorR包(最常用)

http://www.bio-info-trainee.com/1399.html 在里面找自己想找的gpl_number,然后找到对应的R包,下载该R包(记住R包后面是.bd)

代码语言:txt
复制
if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db")
library(hgu133plus2.db)
ls("package:hgu133plus2.db") # 查看该包里面有什么
ids <- toTable(hgu133plus2SYMBOL)# symbol代表的是探针的ID和基因symbol,toTable是提取
head(ids)

方法2 读取GPL网页的表格文件,按列取子集

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570

代码语言:txt
复制
if(F){
  b = read.delim("GPL570-55999.txt",
               check.names = F,
               comment.char = "#")# read.delim是读取该文件
  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,]
}

方法3 官网下载注释文件并读取

http://www.affymetrix.com/support/technical/byproduct.affx?product=hg-u133-plus

方法4 自主注释

https://mp.weixin.qq.com/s/mrtjpN8yDKUdCSvSUuUwcA

保存

代码语言:javascript
复制
save(exp,Group,ids,gse_number,file = "step2output.Rdata")

七、PCA图

代码语言:txt
复制
dat=as.data.frame(t(exp))#转置,将数据框的横纵左边转置变成矩阵,之后再as.data.frame转成数据框
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")

八、热图

比较的应该是每行固定基因在不同样本中的表达量差异。用标准差大的基因进行画图,聚类和分组的差别可能会大一些,但若选择表达差异大的基因,聚类和分组会更一致。

代码语言:txt
复制
cg=names(tail(sort(apply(exp,1,sd)),1000))# 把标准差最大的1000个基因挑选出来,标准差大不一定是差异基因,只能代表其是活跃基因
n=exp[cg,]
代码语言:javascript
复制
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
         cluster_cols=F# 意思是不进行聚类,热图的顺序就是分组的顺序
)#这样得到的热图是表达矩阵里的所有数据都进行作图

按行标准化

代码语言:javascript
复制
pheatmap(n,
         show_colnames =F,
         show_rownames = F,
         annotation_col=annotation_col,
         scale = "row",
         breaks = seq(-3,3,length.out = 100)#设置色带分布范围,从-3到3,等差数列增加,形成100个色带值
         )                                #一般设置色带分布范围不能太大,这样太离群的数据会被突出显示,而正常数据差异区分的不明显。
dev.off()

九、差异分析

差异分析,用limma包来做,DEG:差异表达基因

需要表达矩阵和Group

代码语言:javascript
复制
library(limma)
design=model.matrix(~Group)
fit=lmFit(exp,design)
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf)
#这几行代码是做差异分析的,不需要改的,直接用

为deg数据框添加几列

1.加probe_id列,把行名变成一列

代码语言:javascript
复制
library(dplyr)
deg <- mutate(deg,probe_id=rownames(deg))

2.加上探针注释

代码语言:javascript
复制
ids = ids[!duplicated(ids$symbol),]# 或者用这种随机去重方式ids = diatinct(ids,symbol,.keep_all=T)
#其他去重方式在zz.去重方式.R
deg <- inner_join(deg,ids,by="probe_id")
nrow(deg)

若一个基因对应多个探针,解决办法有:

1.随机去重(推荐使用)

2.保留行和/行平均值最大的探针

3.取多个探针的平均值

3.加change列,标记上下调基因,用ifelse

代码语言:javascript
复制
logFC_t=1
p_t = 0.05 #设置logFC和p-value的阈值,把阈值调大一些,差异基因的数量就会相对增加一些
k1 = (deg$adj.P.Val < p_t)&(deg$logFC < -logFC_t)#这里用agj.p.value做纵坐标
k2 = (deg$adj.P.Val < p_t)&(deg$logFC > logFC_t)#k1是下调基因的条件,k2是上调基因的条件
deg <- mutate(deg,change = ifelse(k1,"down",ifelse(k2,"up","stable")))
table(deg$change)

4.加ENTREZID列,用于富集分析

symbol转entrezid用的是bitr,然后inner_join),基因的命名有很多种,常见的有GeneSymbol、ENSEMBLID、EntrezID等,为了在不同的基因命名方式之间快速转换,使用OrgDb。

代码语言:javascript
复制
library(clusterProfiler)
library(org.Hs.eg.db)
s2e <- bitr(deg$symbol, 
            fromType = "SYMBOL",
            toType = "ENTREZID",
            OrgDb = org.Hs.eg.db)#人类
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))#将增加的那一列添加到表达数据框中,by=c(a,b)是Merge的两个数据框中列名的大小写不一样,这个操作是统一大小写
save(Group,deg,logFC_t,p_t,gse_number,file = "step4output.Rdata")

其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb

十、火山图

代码语言:javascript
复制
library(dplyr)
library(ggplot2)
p <- ggplot(data = dat, 
            aes(x = logFC, 
                y = -log10(adj.P.Val))) +
  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

dat  = deg[!duplicated(deg$symbol),]#去重复


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

十一、差异基因热图

代码语言:javascript
复制
load(file = 'step2output.Rdata')
# 表达矩阵行名替换
exp = exp[dat$probe_id,]
rownames(exp) = dat$symbol#把探针的id换成基因的名字
if(T){
  #取前10上调和前10下调
  library(dplyr)
  dat2 = dat %>%
    filter(change!="stable") %>% 
    arrange(logFC) 
  cg = c(head(dat2$symbol,10),#前10个
         tail(dat2$symbol,10))#后10个
}else{
  # 全部差异基因
  cg = dat$symbol[dat$change !="stable"]
  length(cg)
}
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",
                         #show_rownames=F,这行代码是运行全部差异基因时,去掉热图上基因的名字,要不然都叠到一起了
                         #cluster_cols = F, 这行代码是把聚类和分组对应
                         annotation_col=annotation_col,
                         breaks = seq(-3,3,length.out = 100)
) 
heatmap_plot

拼图

代码语言:javascript
复制
library(patchwork)
volcano_plot +heatmap_plot$gtable#热图和火山图拼到一起

十二、感兴趣的基因相关性

代码语言:javascript
复制
library(corrplot)
g = sample(deg$symbol[1:500],10) # 这里是随机取样,注意换成自己感兴趣的基因,从前500行中随机取10个
g

M = cor(t(exp[g,])) # t(exp[g,])是转置,行变成列,然后cor()计算列与列之间的相关性
pheatmap(M)

library(paletteer)#配色R包
my_color = rev(paletteer_d("RColorBrewer::RdYlBu")) # paletteer_d()是包里的颜色,rev()就逆转了包里颜色的顺序
my_color = colorRampPalette(my_color)(10) #把这一组颜色切分成10中颜色
corrplot(M, type="upper", #画的cowplot是不能赋值的
         method="pie",
         order="hclust", 
         col=my_color,
         tl.col="black", 
         tl.srt=45)
library(cowplot) 
cor_plot <- recordPlot() #只有当画图框里有图才能被记录

拼图

代码语言:javascript
复制
load("pca_plot.Rdata")
plot_grid(pca_plot,cor_plot,
          volcano_plot,heatmap_plot$gtable)
dev.off()

保存

代码语言:javascript
复制
pdf("deg.pdf")
plot_grid(pca_plot,cor_plot,
          volcano_plot,heatmap_plot$gtable)
dev.off()

十三、富集分析

提前下载的包

代码语言:javascript
复制
load(file = 'step4output.Rdata')
library(clusterProfiler)
library(ggthemes)
library(org.Hs.eg.db)
library(dplyr)
library(ggplot2)
library(stringr)
library(enrichplot)

1.输入数据

代码语言:javascript
复制
gene_up = deg$ENTREZID[deg$change == 'up'] 
gene_down = deg$ENTREZID[deg$change == 'down'] 
gene_diff = c(gene_up,gene_down)#合并上调基因和下调基因

2.富集

以下步骤耗时很长,设置了存在即跳过

代码语言:javascript
复制
f = paste0(gse_number,"_GO.Rdata");f
if(!file.exists(f)){
  ego <- enrichGO(gene = gene_diff,
                  OrgDb= org.Hs.eg.db,#物种来源
                  ont = "ALL",
                  readable = TRUE)#entrezid又变成sample,具有可读性
  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.可视化

条带图

代码语言:javascript
复制
barplot(ego)
barplot(ego, split = "ONTOLOGY", font.size = 10, 
        showCategory = 5) + 
  facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y") 

气泡图

代码语言:javascript
复制
dotplot(ego)

dotplot(ego, split = "ONTOLOGY", font.size = 10, 
        showCategory = 5) + 
  facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y") 

4.展示top通路的共同基因,要放大看

代码语言:javascript
复制
#gl 用于设置下图的颜色
gl = deg$logFC
names(gl)=deg$ENTREZID
#Gene-Concept Network,要放大看
cnetplot(ego,
         #layout = "star",
         color.params = list(foldChange = gl),
         showCategory = 3)

十四、KEGG pathway analysis

上调、下调、差异、所有基因

1.输入数据

代码语言:javascript
复制
gene_up = deg[deg$change == 'up','ENTREZID'] 
gene_down = deg[deg$change == 'down','ENTREZID'] 
gene_diff = c(gene_up,gene_down)

2.对上调/下调/所有差异基因进行富集分析

代码语言:javascript
复制
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

代码语言:javascript
复制
table(kk.diff@result$p.adjust<0.05)
table(kk.up@result$p.adjust<0.05)
table(kk.down@result$p.adjust<0.05)
#实验数据在KEGG中没有富集到很正常,因为KEGG的数据库本来就很小。

4.双向图

富集分析所有图表默认都是用p.adjust,富集不到可以退而求其次用p值,在文中说明即可

代码语言:javascript
复制
source("kegg_plot_function.R")#在不打开该文件的前提下全选运行
g_kegg <- kegg_plot(kk.up,kk.down)
g_kegg

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

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

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

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

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 一、数据下载
  • 二、提取临床信息
  • 三、让exp列名与pd的行名顺序完全一致
  • 四、提取芯片平台编号
  • 五、Group(实验分组)和ids(探针注释)
  • 六、探针注释的获取
    • 方法1 BioconductorR包(最常用)
      • 方法2 读取GPL网页的表格文件,按列取子集
        • 方法3 官网下载注释文件并读取
          • 方法4 自主注释
          • 七、PCA图
          • 八、热图
          • 九、差异分析
            • 1.加probe_id列,把行名变成一列
              • 2.加上探针注释
                • 3.加change列,标记上下调基因,用ifelse
                  • 4.加ENTREZID列,用于富集分析
                  • 十、火山图
                  • 十一、差异基因热图
                  • 十二、感兴趣的基因相关性
                  • 十三、富集分析
                    • 1.输入数据
                      • 2.富集
                        • 3.可视化
                          • 条带图
                          • 气泡图
                        • 4.展示top通路的共同基因,要放大看
                        • 十四、KEGG pathway analysis
                          • 1.输入数据
                            • 2.对上调/下调/所有差异基因进行富集分析
                              • 3.看看富集到了吗?
                                • 4.双向图
                                领券
                                问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档