R语言GEO数据挖掘01-数据下载及提取表达矩阵
欢迎来到医科研，这里是白介素2的读书笔记，跟我一起聊临床与科研的故事, 生物医学数据挖掘，R语言，TCGA、GEO数据挖掘。
这一节的内容包括应用
GEOquery包下载芯片数据，提取表达矩阵，提取metadata信息。解决一个探针对应多个基因的问题
GEO数据下载-GEOquery
安装GEOquery包
1options(stringsAsFactors = F)##避免将character转换为因子
2if (!requireNamespace("BiocManager", quietly = TRUE))
3 install.packages("BiocManager")
4if(!require("GEOquery")) BiocManager::install("GEOquery")
5library(GEOquery)
6library(dplyr)
7browseVignettes("GEOquery")##获取帮助
选择一个数据集GSE7765演示分析
数据集基本情况
芯片平台GPL96,GPL97; sample数12个
下载表达矩阵
1gse <- getGEO("GSE7765", GSEMatrix = TRUE)
2show(gse)
GSE7765中包括两个平台，两个数据集
提取表达矩阵及metadata
1class(gse)
2str(gse)
3a<-gse[[1]]
4b<-gse[[2]]
5class(gse[[1]])##ExpressionSet
6
7##提取第一个数据集的phenodata
8dim(pData(gse[[1]]))
9metdata<-pData(gse[[1]])
10metdata[1:5,1:5]
11colnames(metdata)##phenodata信息很多，但用得上的很少
12
13##提取第一个表达矩阵
14expma<-exprs(a)
15dim(expma)
16expma[1:5,1:5]
17save(metdata,expma,file = "expma.Rdata")
18 GSM188013 GSM188014 GSM188016 GSM188018 GSM188020
191007_s_at 15630.200 17048.800 13667.500 15138.800 10766.600
201053_at 3614.400 3563.220 2604.650 1945.710 3371.290
21117_at 1032.670 1164.150 510.692 5061.200 452.166
22121_at 5917.800 6826.670 4562.440 5870.130 3869.480
231255_g_at 224.525 395.025 207.087 164.835 111.609
平台注释信息处理
芯片数据分析中很重要的内容即平台信息处理，获取相应的平台
这里我们选择GPL96
1if(F){load(file="expma.Rdata")}
2GPL="GPL96"##下载平台注释
3gpl <- getGEO(GPL,destdir = getwd()) %>%
4 Table() %>% ##转换为data.frame格式
5 save(file = "GPL96_annot.Rdata")
6if(F){load(file = "GPL96_annot.Rdata")}
7head(gpl)
8colnames(gpl)
9
10##取出注释信息
11probe<-gpl %>%
12 select("ID","Gene Symbol","ENTREZ_GENE_ID")
13head(probe)
14dim(probe)##22283个
15 ID Gene Symbol ENTREZ_GENE_ID
161 1007_s_at DDR1 /// MIR4640 780 /// 100616237
172 1053_at RFC2 5982
183 117_at HSPA6 3310
194 121_at PAX8 7849
205 1255_g_at GUCA1A 2978
216 1294_at MIR5193 /// UBA7 7318 /// 100847079
22[1] 22283 3
23Modify 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
24Show in New WindowClear OutputExpand/Collapse Output
25[1] 20878 3
到这里我们发现一个情况：1个探针存在对应多个基因的情况，///
我们的解决思路是有几种，第一种是直接将存在///的信息去掉
第二种是将///的数据拆开来，然后再把有重复的删去
先按第一种解决，这种解决比较简单，我的需要是只要找出///
1probe<-probe[!grepl(" /// ",probe$`Gene Symbol`),]
2dim(probe)##数量减少到20878
第二种方法是，先拆解再删除
拆解时需要将///分割开，再与ID相连
值得注意的是这种方法其实与第一种是有区别的，这种方法仍然保留了探针对应多个基因中的
一种情况，所有得出的probe注释要多，这里不纠结这个内容
1if(F){
2library(tidyverse)
3probe2<-apply(probe,1,function(x){
4 paste(x[1],
5 str_split(x[2]," /// ",simplify = T),##分割
6 sep = "|")##连接符号
7 }
8) %>%
9 unlist()##得到的是个list
10
11###
12head(probe2)
13length(probe2)##展开后得到24807个探针及对应关系
14probe2<-as_tibble(probe2)
15
16##注意这里的\\是用于转义 匹配分割 "|"，达到分割目的
17probe2<-probe2 %>% separate(value,c("ID","GeneName"),sep = "\\|")
18dim(probe2)##增加到24807行
19
20## 找出重复ID，两个table的妙用
21table(table(probe2$ID))##探针找出对应一个基因的有20878个，与grepl法得出的结果相同
22
23## 下一步的目的即筛选出对应一个基因的探针
24test2<-probe2 %>% count(ID) %>% filter(n==1) %>% ## count计数有点类似于table
25 inner_join(probe2,by="ID") %>% ## 内连接只保留x,y中观测相同的变量
26 select(-n)##remove "n" column
27dim(test2)
28head(test2)
29probe2<-test2##将最终得到的结果赋值给probe2
30}