从零开始的异世界生信学习 GEO数据库数据挖掘--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) 
}

#清空环境，加载R包
rm(list = ls())
library(GEOquery)
#下载GEO数据先去GEO数据库网页确定是否是表达芯片数据，不是的话不能用本流程。
gse_number = "GSE56649"
eSet <- getGEO(gse_number, destdir = '.', getGPL = F) 
##getGEO函数可以下载到工作目录下和读取GSE文件，
class(eSet)
length(eSet)
eSet = eSet[[1]]

#(1)提取表达矩阵exp
exp <- exprs(eSet)
##exprs函数是expressionset文件中提取表达矩阵的函数
dim(exp)
exp[1:4,1:4]
#检查矩阵是否正常，如果是空的就会报错，空的和有负值的、有异常值的矩阵需要处理原始数据。
#如果表达矩阵为空，大多数是转录组数据，不能用这个流程(后面另讲)。
#自行判断是否需要log
exp = log2(exp+1)
## 这个步骤在表达矩阵中表达量+1为了防止数据中有0，取log2后出现负无穷
boxplot(exp)
##exp = limma::normalizeBetweenArrays(exp) 可以通过这句代码进行对表达矩阵处理
#(2)提取临床信息
pd <- pData(eSet)
##表达矩阵的列名和临床信息的行名必须一致才能进行后续分析操作
#(3)让exp列名与pd的行名顺序完全一致
p = identical(rownames(pd),colnames(exp));p
if(!p) exp = exp[,match(rownames(pd),colnames(exp))]
#(4)提取芯片平台编号
gpl_number <- eSet@annotation;gpl_number ##狭义的数据类型取子集需要用@或者$符号，可以进行实验能否补齐
save(gse_number,pd,exp,gpl_number,file = "step1output.Rdata")

rm(list = ls())
library(AnnoProbe)
library(GEOquery)
gse_number = "GSE56649"
eSet = geoChina(gse_number)

# Group(实验分组)和ids(探针注释)
rm(list = ls())  
load(file = "step1output.Rdata")
library(stringr)
# 标准流程代码是二分组，多分组数据的分析后面另讲
# 生成Group向量的三种常规方法，三选一，选谁就把第几个逻辑值写成T，另外两个为F。如果三种办法都不适用，可以继续往后写else if
if(F){
  # 1.Group----一般实验分组为一个单词
  # 第一种方法，有现成的可以用来分组的列
  Group = pd$`disease state:ch1` 
   ## pd$后tab补齐，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$source_name_ch1,"control"),
               "control",
               "RA")
  ##Group=ifelse(str_detect(pd$source_name_ch1,"patient")，"RA", "control")
}
## str_detect用来搜索关键词
##第三种方法需要临床信息中的字符串中有分组信息的文字

# 需要把Group转换成因子，并设置参考水平，指定levels，对照组在前，处理组在后
Group = factor(Group,levels = c("control","RA"))
Group
## factor（Group）生成因子是默认按照首字母顺序排序
##Group = factor(Group,levels = c("control","RA")) 按照代码中的顺序进行排序，control组在第一个位置上

#2.探针注释的获取-----------------
#捷径，以下代码可以检索使用那种方法获得基因注释
library(tinyarray)
find_anno(gpl_number) #打出找注释的代码
ids <- AnnoProbe::idmap('GPL570')

#四种方法，方法1里找不到就从方法2找，以此类推。
#方法1 BioconductorR包(最常用)
gpl_number 
#http://www.bio-info-trainee.com/1399.html
#获取了GPL编号后，登陆网站，搜索使用那个R包进行注释，注意R包名称后面有.db后缀
if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db")
library(hgu133plus2.db)
ls("package:hgu133plus2.db") ##加载R包后，查看R包中哪部分是所需要的注释，R包无法自动补齐，注意
ids <- toTable(hgu133plus2SYMBOL)  ##使用toTable函数加载R包中的SYMBOL，并生成数据框
head(ids)

# 方法2 读取GPL网页的表格文件，按列取子集
##https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570
if(F){
  #注：表格读取参数、文件列名不统一，活学活用，有的表格里没有symbol列，也有的GPL平台没有提供注释表格
  #read.delim函数是read.table的替代函数，用来读取，sep默认\t,默认header=T
  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) ##symbol列部分的空格为空字符串，取不要空格的行
  k2 = !str_detect(ids2$symbol,"///");table(k2) ##
  ids2 = ids2[ k1 & k2,]
  # ids = ids2
}
##GPL网站下载的表格文件中可能存在多余的行，探针没有对应genesymbol

# 方法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")