Posted on

转载自:R语言GEO数据挖掘01-数据下载及提取表达矩阵 – 简书

  • 这一节的内容包括应用 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 % 
  Table() %>% ##转换为data.frame格式
  save(file = "GPL96_annot.Rdata")
if(F){load(file = "GPL96_annot.Rdata")}
head(gpl)
colnames(gpl)

##取出注释信息
probe% 
  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% 
  unlist()##得到的是个list

###
head(probe2)
length(probe2)##展开后得到24807个探针及对应关系
probe2% separate(value,c("ID","GeneName"),sep = "\|")
dim(probe2)##增加到24807行

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

## 下一步的目的即筛选出对应一个基因的探针
test2% 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
}

7 Replies to “GEO数据下载及提取表达矩阵”

  1. 你好,麻烦问一下如果表达矩阵里,范围是-3~20这是已经log过的吗?如果是这样的话2的20次方=1048576这样也太大了

    1. 您好,有负值的话是log2转化过的,芯片的原始数据是没有负值的,0-1的数据log转化之后就是负值。
      具体的标准化信息在你下载的数据里也可以看到

      1. 好的,谢谢进哥哥。麻烦再问一下怎么查看下载数据的处理方式呢?非常感谢

  2. 博主我用GEO的时候发现他们的箱式图不太整齐。我不确定是不是还要进行normalizeBetweenArrays

      1. 博主我看你新发的贴子, 如果我不是两个GEO数据集合并, 就一个数据集使用,箱式图不是很齐 我需要拉平吗? 我试了拉平和不拉平的差异很大

        1. 当然需要的,可能原因是没有进行log转化,导致数据中位数差得很大,如果发现数值有上千的 先log2转化 然后再看箱式图齐不齐 不齐在进行normalizeBetweenArrays。

发表评论

邮箱地址不会被公开。 必填项已用*标注