R语言GEO数据挖掘01-数据下载及提取表达矩阵

作者：白介素2

• 这一节的内容包括应用 GEOquery包下载芯片数据，提取表达矩阵，提取metadata信息。
• 解决一个探针对应多个基因的问题

GEO数据下载-GEOquery

安装GEOquery包

options(stringsAsFactors = F)##避免将character转换为因子
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
if(!require("GEOquery")) BiocManager::install("GEOquery")
library(GEOquery)
library(dplyr)
browseVignettes("GEOquery")##获取帮助

选择一个数据集GSE7765演示分析

数据集基本情况
芯片平台GPL96,GPL97; sample数12个

下载表达矩阵

gse <- getGEO("GSE7765", GSEMatrix = TRUE) 
show(gse)

GSE7765中包括两个平台，两个数据集

提取表达矩阵及metadata

class(gse)
str(gse)
a<-gse[[1]]
b<-gse[[2]]
class(gse[[1]])##ExpressionSet

##提取第一个数据集的phenodata
dim(pData(gse[[1]]))
metdata<-pData(gse[[1]])
metdata[1:5,1:5]
colnames(metdata)##phenodata信息很多，但用得上的很少

##提取第一个表达矩阵
expma<-exprs(a)
dim(expma)
expma[1:5,1:5]
save(metdata,expma,file = "expma.Rdata")
          GSM188013 GSM188014 GSM188016 GSM188018 GSM188020
1007_s_at 15630.200 17048.800 13667.500 15138.800 10766.600
1053_at    3614.400  3563.220  2604.650  1945.710  3371.290
117_at     1032.670  1164.150   510.692  5061.200   452.166
121_at     5917.800  6826.670  4562.440  5870.130  3869.480
1255_g_at   224.525   395.025   207.087   164.835   111.609

平台注释信息处理

芯片数据分析中很重要的内容即平台信息处理，获取相应的平台

这里我们选择GPL96

if(F){load(file="expma.Rdata")}
GPL="GPL96"##下载平台注释
gpl <- getGEO(GPL,destdir = getwd()) %>% 
  Table() %>% ##转换为data.frame格式
  save(file = "GPL96_annot.Rdata")
if(F){load(file = "GPL96_annot.Rdata")}
head(gpl)
colnames(gpl)

##取出注释信息
probe<-gpl %>% 
  select("ID","Gene Symbol","ENTREZ_GENE_ID")
head(probe)  
dim(probe)##22283个
         ID      Gene Symbol     ENTREZ_GENE_ID
1 1007_s_at DDR1 /// MIR4640  780 /// 100616237
2   1053_at             RFC2               5982
3    117_at            HSPA6               3310
4    121_at             PAX8               7849
5 1255_g_at           GUCA1A               2978
6   1294_at MIR5193 /// UBA7 7318 /// 100847079
[1] 22283     3
Modify Chunk OptionsRun All Chunks AboveRun Current ChunkModify Chunk OptionsRun All Chunks AboveRun Current ChunkModify Chunk OptionsRun All Chunks AboveRun Current ChunkModify Chunk OptionsRun All Chunks AboveRun Current Chunk
Show in New WindowClear OutputExpand/Collapse Output
[1] 20878     3

到这里我们发现一个情况：1个探针存在对应多个基因的情况，///

我们的解决思路是有几种，第一种是直接将存在///的信息去掉
第二种是将///的数据拆开来，然后再把有重复的删去

先按第一种解决，这种解决比较简单，我的需要是只要找出///

probe<-probe[!grepl(" /// ",probe$`Gene Symbol`),]
dim(probe)##数量减少到20878

第二种方法是，先拆解再删除

拆解时需要将///分割开，再与ID相连
值得注意的是这种方法其实与第一种是有区别的，这种方法仍然保留了探针对应多个基因中的
一种情况，所有得出的probe注释要多，这里不纠结这个内容

if(F){
library(tidyverse)
probe2<-apply(probe,1,function(x){
  paste(x[1],
  str_split(x[2]," /// ",simplify = T),##分割
  sep = "|")##连接符号
 }  
) %>% 
  unlist()##得到的是个list

###
head(probe2)
length(probe2)##展开后得到24807个探针及对应关系
probe2<-as_tibble(probe2)

##注意这里的\\是用于转义 匹配分割 "|"，达到分割目的
probe2<-probe2 %>% separate(value,c("ID","GeneName"),sep = "\\|")
dim(probe2)##增加到24807行

## 找出重复ID，两个table的妙用
table(table(probe2$ID))##探针找出对应一个基因的有20878个，与grepl法得出的结果相同

## 下一步的目的即筛选出对应一个基因的探针
test2<-probe2 %>% count(ID) %>% filter(n==1) %>% ## count计数有点类似于table
  inner_join(probe2,by="ID") %>% ## 内连接只保留x,y中观测相同的变量
  select(-n)##remove "n" column
dim(test2)
head(test2)
probe2<-test2##将最终得到的结果赋值给probe2
}

本节的内容就到这里，我是白介素2，下期再见
转载请注明出处