作者：白介素2

• 這一節(jié)的內(nèi)容包括應(yīng)用 GEOquery包下載芯片數(shù)據(jù)，提取表達(dá)矩陣，提取metadata信息。
• 解決一個(gè)探針對(duì)應(yīng)多個(gè)基因的問題

GEO數(shù)據(jù)下載-GEOquery

安裝GEOquery包

options(stringsAsFactors = F)##避免將character轉(zhuǎn)換為因子
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
if(!require("GEOquery")) BiocManager::install("GEOquery")
library(GEOquery)
library(dplyr)
browseVignettes("GEOquery")##獲取幫助

選擇一個(gè)數(shù)據(jù)集GSE7765演示分析

數(shù)據(jù)集基本情況
芯片平臺(tái)GPL96,GPL97; sample數(shù)12個(gè)

下載表達(dá)矩陣

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

GSE7765中包括兩個(gè)平臺(tái)，兩個(gè)數(shù)據(jù)集

提取表達(dá)矩陣及metadata

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

##提取第一個(gè)數(shù)據(jù)集的phenodata
dim(pData(gse[[1]]))
metdata<-pData(gse[[1]])
metdata[1:5,1:5]
colnames(metdata)##phenodata信息很多，但用得上的很少

##提取第一個(gè)表達(dá)矩陣
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

平臺(tái)注釋信息處理

芯片數(shù)據(jù)分析中很重要的內(nèi)容即平臺(tái)信息處理，獲取相應(yīng)的平臺(tái)

這里我們選擇GPL96

if(F){load(file="expma.Rdata")}
GPL="GPL96"##下載平臺(tái)注釋
gpl <- getGEO(GPL,destdir = getwd()) %>% 
  Table() %>% ##轉(zhuǎn)換為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個(gè)
         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)一個(gè)情況：1個(gè)探針存在對(duì)應(yīng)多個(gè)基因的情況，///

我們的解決思路是有幾種，第一種是直接將存在///的信息去掉
第二種是將///的數(shù)據(jù)拆開來，然后再把有重復(fù)的刪去

先按第一種解決，這種解決比較簡單，我的需要是只要找出///

probe<-probe[!grepl(" /// ",probe$`Gene Symbol`),]
dim(probe)##數(shù)量減少到20878

第二種方法是，先拆解再刪除

拆解時(shí)需要將///分割開，再與ID相連
值得注意的是這種方法其實(shí)與第一種是有區(qū)別的，這種方法仍然保留了探針對(duì)應(yīng)多個(gè)基因中的
一種情況，所有得出的probe注釋要多，這里不糾結(jié)這個(gè)內(nèi)容

if(F){
library(tidyverse)
probe2<-apply(probe,1,function(x){
  paste(x[1],
  str_split(x[2]," /// ",simplify = T),##分割
  sep = "|")##連接符號(hào)
 }  
) %>% 
  unlist()##得到的是個(gè)list

###
head(probe2)
length(probe2)##展開后得到24807個(gè)探針及對(duì)應(yīng)關(guān)系
probe2<-as_tibble(probe2)

##注意這里的\\是用于轉(zhuǎn)義 匹配分割 "|"，達(dá)到分割目的
probe2<-probe2 %>% separate(value,c("ID","GeneName"),sep = "\\|")
dim(probe2)##增加到24807行

## 找出重復(fù)ID，兩個(gè)table的妙用
table(table(probe2$ID))##探針找出對(duì)應(yīng)一個(gè)基因的有20878個(gè)，與grepl法得出的結(jié)果相同

## 下一步的目的即篩選出對(duì)應(yīng)一個(gè)基因的探針
test2<-probe2 %>% count(ID) %>% filter(n==1) %>% ## count計(jì)數(shù)有點(diǎn)類似于table
  inner_join(probe2,by="ID") %>% ## 內(nèi)連接只保留x,y中觀測(cè)相同的變量
  select(-n)##remove "n" column
dim(test2)
head(test2)
probe2<-test2##將最終得到的結(jié)果賦值給probe2
}

本節(jié)的內(nèi)容就到這里，我是白介素2，下期再見
轉(zhuǎn)載請(qǐng)注明出處