【生信技能树】GEO数据库挖掘

P5 4 ID转换技巧大全

【方法1】本章接导入，利用read.table函数导入gz文件后是一个data.frame

# 根据情况设置分隔符啥的
  a=read.table('GSE42872_series_matrix.txt.gz', sep = '\t',quote = "", fill = T,
               comment.char = "!", header = T)
  # 获取第一列的名字
  rownames(a)=a[,1]
  # 去掉第一列，避免重复
  a=a[,-1]
  class(a)
  str(a)

【方法2】本章利用getGEO方式导入，获取其表达矩阵。

如果使用getGEO获取的是一个对象，要取其中的[[1]]

这个GSE的文件包含了很多个平台的测序数据，需要str看一下他内部有几个对象，即有哪几个平台的数据在，然后获取gset[[1]]，这个相关信息也可以在这个GSE的样本界面去寻找。

找到这个对象，然后再对其进行expr转化为matrix

if(!file.exists(f)){
  gset <- getGEO('GSE42872', destdir=".",
                 AnnotGPL = F,     ## 注释文件
                 getGPL = F)       ## 平台文件
  save(gset,file=f)   ## 保存到本地，并且被直接命名为gset
}
b=gset[[1]]
a1=exprs(b) # 获得matrix

我们可以将代码打包成一个function：

downGSE <- function(studyID="GSE42872",destdir="."){
  library(GEOquery)

      gset <- getGEO(studyID, destdir=".",
                     AnnotGPL = F,     ## 注释文件
                     getGPL = F)       ## 平台文件
    # 对gset处理
    exprSet=exprs(gset[[1]]) # 获取表达矩阵
    pdata=pData(gset[[1]]) # 获取表达对应信息
    class(gset)
    length(gset)
    str(gset)
    write.csv(exprSet,paste0(studyID,"_exprSet.csv")) # 工作环境中导出表达矩阵
    write.csv(pdata,paste0(studyID,"_metadata.csv")) # 工作环境中导出相关信息
    reture(gset)
  }
# 该种function可以直接获取标本信息（pdata）以及修正过的表达矩阵
downGSE('GSE42872')

导入后需注意格式，必要时调整。

以下为获取包以及探针的内容，可在另一篇笔记中获取，等有空再补充学习完整其他方法。

# 然后看一下这个GSE数据所对应的平台的ID文件
library(hugene10sttranscriptcluster.db)
# 此时回到官网去搜这个包的GPL平台，进而去搜这个ID注释的包
ls("package:hugene10sttranscriptcluster.db") 

ids=toTable(hugene10sttranscriptclusterSYMBOL)
# 用exprSet来后续计算，方便一些。
exprSet <- as.matrix(GSE42872_exprSet)
# 注释文件命名为pdata
pdata <- read_csv("GSE42872_metadata.csv")
save(ids,exprSet,GSE42872_metadata,file = 'input.Rdata')
length(unique(ids$symbol))
tail(sort(table(ids$symbol)))
table(sort(table(ids$symbol)))

plot(table(sort(table(ids$symbol))))
# ids就是探针对应基因的list,ids中是基因的symble
# 由于可能每个symble也可能对应多个probe，即很多个symble后面都是相同的基因
# 因此要选取同一个symble（基因）对应多个probe中的处理过的量
# 后续在笔记中也可见。等时机成熟补充完全。

需要观察一下数据，有的symbol对应多个probe，因此需要洗一下。

Step-1：先选取这个包里有的probe的基因

# 清洗exprset中与id有交集的基因，%in%就是all in得意思
table(rownames(exprSet) %in% ids$probe_id)
dim(exprSet)
# 选取all in的基因
exprSet=exprSet[rownames(exprSet) %in% ids$probe_id,]
dim(exprSet)

Step-2：按exprset的顺序排列探针

# 更改id顺序，从ids中提取表达矩阵的每一个，这样就一致了
ids=ids[match(rownames(exprSet),ids$probe_id),]
# 注意ids和exprset的行数是否相同
head(ids)
# 查看一下
exprSet[1:5,1:5]

举个例子：用IGKC这个有很多probe的symbol来说明后续函数的运行逻辑

# 通过id的symbol对表达矩阵分类
# 用一个10个探针对应一个基因的情况来举例子
# IGKC基因：这里要在上述取相同基因之前学习
if(F){
  x=exprSet[ids[,2]=='IGKC',] # exprset的行==ids的列
  rownames(x)
  rowMeans(x) #  每行取平均值，因为这里是挑probe所以要看平均的表达量即可。
  which.max(rowMeans(x)) # 看看最大的在哪里
  rownames(x)
}

正式的代码：

# 获取新探针列表
  # by函数：把exprset数据集按照ids的symbol分类
  # 然后每个新数据集，做function：取这里的所有行平均值的最大值那个对应的rownames
  # 这里就可以根据需求自行调节
  tmp = by(exprSet,
           ids$symbol,
           function(x) rownames(x)[which.max(rowMeans(x))] )
  # 提取这个列表，教做probes
  probes = as.character(tmp)
  dim(exprSet)
  exprSet=exprSet[rownames(exprSet) %in% probes ,]
  dim(exprSet)
  
  rownames(exprSet)=ids[match(rownames(exprSet),ids$probe_id),2]
  exprSet[1:5,1:5]
  
  # 写一个函数包含上述的内容
  # 这里exprset和ids必须一一对应。
  if(F){
    cleaning <- function(exprSet,ids){
      tmp = by(exprSet,
               ids$symbol,
               function(x) rownames(x)[which.max(rowMeans(x))] )
    # 提取这个列表，教做probes
    probes = as.character(tmp)
    print(dim(exprSet))
    exprSet=exprSet[rownames(exprSet) %in% probes ,]
    print(dim(exprSet))
    return(exprSet)
    }
  }

在文中可以做一个function进行封存，但是需注意的是两个列表必须数量相等且名字一致。

这里by函数，以及后续再取个交集（以表达矩阵为核心）不展开。

这样就可以获得清洗后即只有一个最合适的表达量probe的symbol的矩阵，并且对应好了相应的symbol且只有1个。

【如果没有相应的包】

无法获取ids，需回到GPL平台获取，要看这个GSE的GPL信息有没有，可以直接下载，比较大。然后再通过getGEO，看看这个gpl内容，再选取其中有用的列，作为后续的ids进行分析。

# 如果无法找到对应的包
if(F){
  library(Biobase)
  library(GEOquery)
  # download GPL files,put in the directory and load
  # 这里的GPL应该是看GSE的矩阵对应的平台号码是谁,底下的表格里有
  gpl <- getGEO('GPL6480', destdir = ".")
  colnames(Table(gpl))
  # 选择所需的column，就是要探针和基因对应的关系，就可以作为ids进行后续分析
  head(Table(gpl)[,c()])
  write.csv(Table(gpl)[,c()],"name")
}