zhangmj123的个人博客分享 http://blog.sciencenet.cn/u/zhangmj123

博文

GEO分析第一弹

已有 1240 次阅读 2020-4-18 19:36 |个人分类:GEO|系统分类:科研笔记

#step 1--------下载GSE4107的芯片数据--------------- library(Biobase);library(BiocGenerics);library(parallel)#GEOquery的依赖包 library(GEOquery) gse="GSE4107"#将GSE4107赋值给gse eSet <- getGEO(gse,destdir = ".", getGPL = F)#使用getGEO函数下载表达数据(此时并不是表达矩阵) #(1)使用exprs()函数提取表达矩阵 class(eSet) length(eSet) eSet[[1]] length(eSet[[1]]) exp <- exprs(eSet[[1]]) class(exp) colnames(exp) rownames(exp) exp[1:4,1:4] #使用log2将芯片表达数据归一化exp=log2(exp+1) exp=log2(exp+1) exp[1:4,1:4] dim(exp) #(2)使用pData()函数从第一步中下载的eSet提取临床信息 pd <- pData(eSet[[1]]) class(pd) dim(pd) pd[1:4,1:4] colnames(pd) rownames(pd) colnames(exp) #(3)调整pd的行名和exp的列名完全一致 identical(rownames(pd),colnames(exp)) p <- identical(rownames(pd),colnames(exp))#检验pd行名是否和exp列名一致 match(rownames(pd),colnames(exp)) if(!p) exp = exp[,match(rownames(pd),colnames(exp))]#如不一致,需运行此行命令 #(4)提取芯片平台编号,平台编号列在annotation列 gpl <- eSet[[1]]@annotation #保存提取的表达矩阵、临床数据(主要是分组信息)、芯片平台编号供下一步使用 save(gse,exp,pd,gpl,file = "step1out.Rdata") #。。。。。。至此,第一步获取表达数据和分组信息已完成。。。。。。。。。。。。。。。 #。。。。。。。。。。。。。。。。。。。。。。。。。 #step 2------分组------------------- #(1)group_list(实验分组)和ids(芯片注释),主要是字符串操作,每次都需要改 rm(list = ls()) load("step1out.Rdata") library(stringr) pd$title#使用此命令查看pd(临床信息)的分组信息 group_list=ifelse(str_detect(pd$title,"control"),"control","treat") dim(pd) pd[1:4,1:4] #使用factor()函数设置参考水平,对照在前,处理在后,谁在前谁是参考 group_list <- factor(group_list,levels = c("control","treat")) #(2)--------ids--------- #方法1 BioconductorR包 gpl #http://www.bio-info-trainee.com/1399.html #此时需要网页搜索GPL对应的包,找到后就可以找到芯片探针ID和基因ID if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db") library(hgu133plus2.db) ls("package:hgu133plus2.db") #用toTable函数将带有symbol的gene生成数据框 ids <- toTable(hgu133plus2SYMBOL) class(ids) head(ids) #至此获得芯片探针ID和基因symbol #保存exp/ids/group_list/至step2.Rdata,下一步就是PCA/heatmap差异分析了 save(exp,ids,group_list,file = "step2out.Rdata") #--------------美妙的分割线--------------- #载入第一步和第二步的数据 rm(list = ls()) load("step1out.Rdata") load("step2out.Rdata") #PCA差异分析只需要exp(表达矩阵)和group_list(分组信息) #Principal Component Analysis #输入数据:exp和group_list #http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials #画图需要数据库,而exp是矩阵,需将其转化成数据库 #先看一下exp数据,确定横纵坐标,PCA横坐标是GSM(样本),纵坐标是group_list class(exp) colnames(exp) rownames(exp) dat=as.data.frame(t(exp))#先将exp表达矩阵转置,后将其转化成数据框 library(FactoMineR)#画主成分分析图需要加载这两个包 library(factoextra) # pca的统一操作走起 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_list, # color by groups #palette = c("#00AFBB", "#E7B800"), addEllipses = TRUE, # Concentration ellipses legend.title = "Groups" ) pca_plot #使用ggplot 中的ggsave()函数保存PCA ggsave(plot = pca_plot,filename = paste0(gse,"PCA.png")) #将图保存为R.data,以供后续拼图用 save(pca_plot,file = "pca_plot.Rdata") #热图 ------ cg=names(tail(sort(apply(exp,1,sd)),1000)) n=exp[cg,] #绘制热图 annotation_col=data.frame(group=group_list) rownames(annotation_col)=colnames(n) library(pheatmap) heatmap=pheatmap(n, show_colnames =F, show_rownames = F, annotation_col=annotation_col, scale = "row") ggsave(plot = heatmap,filename = paste0(gse,"heatmap.png")) save(heatmap,file = "heatmap.Rdata") save(heatmap,cg,n,file = "heatmap.Rdata") #至此初步分析芯片分组情况结束,下一步是挑选表达差异基因 #-----------------华丽丽的分割线--------------

https://blog.sciencenet.cn/blog-3414136-1229004.html

上一篇:学习R语言的第一步:安装软件+安装包
下一篇:我的第一个转录组分析

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2022-5-21 03:54

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部