GEO数据分析流程之芯片1补充

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) }

#数据下载rm(list = ls())#清空数据library(GEOquery)#加载包#先去网页确定是否是表达芯片数据，不是的话不能用本流程。gse_number = "GSE56649"eSet <- getGEO(gse_number, destdir = '.', getGPL = F)#下载并读取class(eSet)#列表length(eSet)#列表长度eSet = eSet[[1]]#提取列表#(1)提取表达矩阵expexp <- exprs(eSet) #提取dim(exp)exp[1:4,1:4] #正常范围在0-20之间#检查矩阵是否正常，如果是空的就会报错，空的和有负值的、有异常值的矩阵需要处理原始数据。#如果表达矩阵为空，大多数是转录组数据，不能用这个流程(后面另讲)。#自行判断是否需要logexp = log2(exp+1) #这一行下载的数据如果取过log，就不要运行。如果在几百，需要取log。如果在0-20，不用取log。为了避免0值，exp+1.甲基化+0.0001，为了保证取log后不会有负无穷。boxplot(exp)#(2)提取临床信息pd <- pData(eSet) #pd文件是样本信息#(3)让exp列名与pd的行名顺序完全一致p = identical(rownames(pd),colnames(exp));p #exp是表达矩阵，identical表示判断是否一致if(!p) exp = exp[,match(rownames(pd),colnames(exp))] #match#(4)提取芯片平台编号gpl_number <- eSet@annotation;gpl_number#@和$都表示提取子集，可以通过补齐判断，或者用str().save(gse_number,pd,exp,gpl_number,file = "step1output.Rdata")

rm(list = ls())library(GEOquery)library(AnnoProbe)#先去网页确定是否是表达芯片数据，不是的话不能用本流程。gse_number = "GSE56649"eSet = geoChina(gse_number)