#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://wap.sciencenet.cn/blog-3414136-1229004.html

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