二、CEO数据 ID转换方法

一、首先获取探针id 和 基因名和对应关系

1、方法1：通过   ”library(hugene10sttranscriptcluster.db)"；这个包中含有基因的探针id和基因的对应关系

library(hugene10sttranscriptcluster.db)
#toTable这个函数：通过看hgu133plus2.db这个包的说明书知道提取probe_id（探针名）和的对应关系的表达矩阵的函数为toTable
ids=toTable(hugene10sttranscriptclusterSYMBOL) 
head(ids)#head为查看前六行
# probe_id     symbol
# 1  7896746   MTND1P23
# 2  7896754 SEPTIN7P13
# 3  7896759  LINC01128
# 4  7896761     SAMD11
# 5  7896779     KLHL17
# 6  7896798    PLEKHN1

2、方法2：通过   ”library(AnnoProbe)"；

library(AnnoProbe)
GEO_file = 'GSE19188_eSet.Rdata'
if(!file.exists(GEO_file)){
  
  # gset = getGEO("GSE19188",destdir = '.',AnnotGPL = F,getGPL = F)
  gset = geoChina("GSE19188")
  save(gset,file =GEO_file )
}
load("./GSE19188_eSet.Rdata") 

# 提取数据
eset <- gset[[1]]
probes_expr <- exprs(eset)
dim(probes_expr)
# [1] 54675   156

#查看表达矩阵最大值是否超过50,limma 包只能分析数值小于50的样本
#如果数值大于50 ，那么就对数值进行log转化
if(max(probes_expr) >50){
  probes_expr = log2(probes_expr+1)
}

#提取表型数据信息
phenoDat <- pData(eset)

# 对表达芯片的探针进行基因注释
(gpl=eset@annotation)
#[1] "GPL570"
# 检查下载的文件中是否包含注释文件信息
checkGPL(gpl)
# [1] TRUE
# 获取下载的文件的相关说明信息
printGPLInfo(gpl)
# 132                                                           
# V1 "GPL570"                                                      
# V2 "[HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array"
# V3 "Homo sapiens"
# 提取注释信息，获取probe id 与symbel id的对应关系的数据框
probe2gene= AnnoProbe::idmap(gpl)

二、id转换

1、方法1：使用AnnoProbe 包自带的函数 filterEM 函数

缺点：filterEM 函数中，如果遇到多个探针对应1个基因时或者说一个基因有多个探针，会随机选择一个基因对应的探针，能会导致在数据集查找不到该探针对应的基因，导致数据不可靠

# 源码如下
filterEM <- function(probes_expr,probe2gene){
  colnames(probe2gene) <- c("probeid","symbol")
  probe2gene$probeid=as.character(probe2gene$probeid)
  probe2gene$symbol=trimws(probe2gene$symbol)
  # 去除NA
  probe2gene=na.omit(probe2gene)
  # if one probe mapped to many genes, we will only keep one randomly.
  probe2gene=probe2gene[!duplicated(probe2gene$probeid),]
  # 
  probe2gene = probe2gene[probe2gene$probeid %in% rownames(probes_expr),]
  
  probes_expr=probes_expr[as.character(probe2gene$probeid),]
  # probes_expr[1:4,1:4]
  
  probe2gene$median=apply(probes_expr,1,median)
  probe2gene=probe2gene[order(probe2gene$symbol,probe2gene$median,decreasing = T),]
  

   #主要是这里，根据基因名进行去重，随机选取一对，
   #比如说有 1-gene1;2-gene1;我选取了1-gene1，但在数据集中探针是2，不是1，
   #那么就匹配不到这个基因
  probe2gene=probe2gene[!duplicated(probe2gene$symbol),]
  
  # 判断筛选结果中是否有重复值
  # any(duplicated(b$symbol))
  # [1] TRUE
  genes_expr=probes_expr[as.character(probe2gene$probeid),]
  rownames(genes_expr)=probe2gene$symbol
 
  return(genes_expr)
}

2、方法2 ：自定义一个函数，将多个相同基因的样本，根据需要，自定义取平均值或者中位数，这里默认是平均值

# 新方法，自定义处理这些多个基因的情况
# 首先要明确的2个概念：
#1、针对一个基因有多个探针
#2、数据框的行名和列名是唯一的
filterEM2 <- function(probes_expr,probe2gene,method = 'mean'){
  colnames(probe2gene) <- c("probeid","symbol")
  probe2gene$probeid=as.character(probe2gene$probeid)
  probe2gene$symbol=trimws(probe2gene$symbol)
  # 去除NA
  probe2gene=na.omit(probe2gene)
  
  #去除重复的探针名，此时得到的探针是唯一的
  #我使用[!duplicated(),] 结果还是有果还有重复值，不能完全筛选出唯一值，不知道为什莫
  #所以我用unique()代替
  rownames(probe2gene)=probe2gene$probeid
  probe2gene=probe2gene[unique(probe2gene$probeid),]
  #筛选数据集中含有的探针对应的probe2gene
  probe2gene = probe2gene[probe2gene$probeid %in% rownames(probes_expr),]
  
  probes_expr <- as.data.frame(probes_expr)
  
  # 将探针和基因的对应关系的数据框和表达样本的数据，根据探针名进行合并，探针对应的基因多为新数据框的一列，此时会出现
  #symbol 这一列的基因名有多个是重复的
  # 不能直接把symbol这一列作为行名，因为行名是唯一的，不能重复
  genes_exprs <- tibble::rownames_to_column(probes_expr,var = 'probeid') |> merge(probe2gene,by='probeid')
  
  #多同一个基因的样本信息求平均值，或者中位数
  
  #1、删除第一列的探针名信息
  genes_exprs <- genes_exprs[-1]
  #2、用aggregate()分组概括数据框
  genes_exprs <- aggregate(
    genes_exprs[,1:ncol(genes_exprs)-1],
    by=list(genes_exprs$symbol),
    FUN = method,
    na.rm=T)
  #将第一列symbol 转化为行名
  genes_exprs_switch <- tibble::column_to_rownames(genes_exprs,var = 'Group.1')
  
  return(genes_exprs_switch)
}