作者：白介素2
- 這一節(jié)的內容包括應用
GEOquery包下載芯片數(shù)據(jù)挪哄，提取表達矩陣，提取metadata信息蟆盐。- 解決一個探針對應多個基因的問題
GEO數(shù)據(jù)下載-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")##獲取幫助
選擇一個數(shù)據(jù)集GSE7765演示分析
數(shù)據(jù)集基本情況
芯片平臺GPL96,GPL97; sample數(shù)12個
下載表達矩陣
gse <- getGEO("GSE7765", GSEMatrix = TRUE)
show(gse)
GSE7765中包括兩個平臺捐川，兩個數(shù)據(jù)集
提取表達矩陣及metadata
class(gse)
str(gse)
a<-gse[[1]]
b<-gse[[2]]
class(gse[[1]])##ExpressionSet
##提取第一個數(shù)據(jù)集的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
平臺注釋信息處理
芯片數(shù)據(jù)分析中很重要的內容即平臺信息處理，獲取相應的平臺
這里我們選擇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
到這里我們發(fā)現(xiàn)一個情況：1個探針存在對應多個基因的情況闺骚，///
我們的解決思路是有幾種彩扔，第一種是直接將存在///的信息去掉
第二種是將///的數(shù)據(jù)拆開來，然后再把有重復的刪去
先按第一種解決僻爽，這種解決比較簡單虫碉，我的需要是只要找出///
probe<-probe[!grepl(" /// ",probe$`Gene Symbol`),]
dim(probe)##數(shù)量減少到20878
第二種方法是，先拆解再刪除
拆解時需要將///分割開胸梆，再與ID相連
值得注意的是這種方法其實與第一種是有區(qū)別的敦捧，這種方法仍然保留了探針對應多個基因中的
一種情況，所有得出的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計數(shù)有點類似于table
inner_join(probe2,by="ID") %>% ## 內連接只保留x,y中觀測相同的變量
select(-n)##remove "n" column
dim(test2)
head(test2)
probe2<-test2##將最終得到的結果賦值給probe2
}
本節(jié)的內容就到這里秽荤，我是白介素2，下期再見
轉載請注明出處
