geo数据挖掘-2

sunqi

2020/7/11

1.概述

对下载的数据进行处理，提取表达矩阵，并匹配探针信息，基因名教程来自：https://github.com/jmzeng1314/GEO/

2.数据下载

2.1 获得表达数据‘

代码语言：javascript复制

rm(list=ls())
# 设置默认转换因子为否
options(stringsAsFactors = F)
# 目标文件
f='GSE42872_eSet.Rdata'
# 上章的geo包
library(GEOquery)
# 下载文件，如果存在则不进行下载
if(!file.exists(f)){
  gset <- getGEO('GSE42872', destdir=".",
                 AnnotGPL = F,     ## 注释文件
                 getGPL = F)       ## 平台文件
  save(gset,file=f)   ## 保存到本地
}
# 为了避免自动化错误，这里再次导入
load('GSE42872_eSet.Rdata')  ## 载入数据
# 查看数据类型为list
class(gset)

代码语言：javascript复制

## [1] "list"

代码语言：javascript复制

#长度
length(gset)

代码语言：javascript复制

## [1] 1

代码语言：javascript复制

# 因为只有一个平台，所以只有1个列表元素
class(gset[[1]])

代码语言：javascript复制

## [1] "ExpressionSet"
## attr(,"package")
## [1] "Biobase"

代码语言：javascript复制

# 列表的第一个元素为表达矩阵和biobase
gset

代码语言：javascript复制

## $GSE42872_series_matrix.txt.gz
## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 33297 features, 6 samples
##   element names: exprs
## protocolData: none
## phenoData
##   sampleNames: GSM1052615 GSM1052616 ... GSM1052620 (6 total)
##   varLabels: title geo_accession ... cell type:ch1 (34 total)
##   varMetadata: labelDescription
## featureData: none
## experimentData: use 'experimentData(object)'
##   pubMedIds: 24469106
## Annotation: GPL6244

代码语言：javascript复制

sprintf("显示下载的文件有6个样本，22397个位点 /n GPL6244")

代码语言：javascript复制

## [1] "显示下载的文件有6个样本，22397个位点 /n GPL6244"

代码语言：javascript复制

#获取列表元素，
a=gset[[1]]

#exprs函数获取表达矩阵
dat=exprs(a)
dim(dat)#查看维度

代码语言：javascript复制

## [1] 33297     6

代码语言：javascript复制

# GPL6244
# 查看前5行
head(dat[,1:5])

代码语言：javascript复制

##         GSM1052615 GSM1052616 GSM1052617 GSM1052618 GSM1052619
## 7892501    7.24559    6.80686    7.73301    6.18961    7.05335
## 7892502    6.82711    6.70157    7.02471    6.20493    6.76554
## 7892503    4.39977    4.50781    4.88250    4.36295    4.18137
## 7892504    9.48025    9.67952    9.63074    9.69200    9.91324
## 7892505    4.54734    4.45247    5.11753    4.87307    5.15505
## 7892506    6.80701    6.90597    6.72472    6.77028    6.77058

代码语言：javascript复制

# 可以看到列名为样本好，行名为探针名
# 绘制箱式图
boxplot(dat,las=2)

代码语言：javascript复制

#查看临床信息,包含6个患者的34个信息
pd=pData(a)

## 选择需要的临床信息
library(stringr)
# 通过空格分隔title，获得分组信息
group_list=str_split(pd$title,' ',simplify = T)[,4]
table(group_list)

代码语言：javascript复制

## group_list
##     Control Vemurafenib
##           3           3

代码语言：javascript复制

# 3个病例和3个对照

2.2 获得平台信息

代码语言：javascript复制

# 查看平台信息探针信息
# GPL6244
# 需要下载时，改为T
if(F){
  library(GEOquery)
  gpl <- getGEO('GPL6244', destdir=".")
  colnames(Table(gpl))
  head(Table(gpl)[,c(1,12)])
  probe2gene=Table(gpl)[,c(1,12)]
  head(probe2gene)
  save(probe2gene,file='probe2gene.Rdata')
}
# 获得平台的所有探针
load(file='probe2gene.Rdata')



# 需要的时候通过biocmanager进行安装
# library(BiocManager)
# BiocManager::install("hugene10sttranscriptcluster.db")
# 这个包用来处理探针信息，每个平台用到的包不一样
library(hugene10sttranscriptcluster.db)
# 获得平台探针信息对应的基因名
ids=toTable(hugene10sttranscriptclusterSYMBOL) #toTable提取probe_id（探针名）和symbol（基因名）
# 重命名
colnames(ids)=c('probe_id','symbol')
# 去掉基因缺失的探针
ids=ids[ids$symbol != '',]
# 匹配平台探针和样本探针
ids=ids[ids$probe_id %in%  rownames(dat),]

# 按照探针名，选取行,去掉样本中平台上没有的探针
dim(dat)

代码语言：javascript复制

## [1] 33297     6

代码语言：javascript复制

dat=dat[ids$probe_id,]
dim(dat)

代码语言：javascript复制

## [1] 19814     6

代码语言：javascript复制

#对dat按行操作，取每一行的中位数将结果返回到median
ids$median=apply(dat,1,median)
#对ids$symbol按照ids$median中位数从大到小排列的顺序排序
ids=ids[order(ids$symbol,ids$median,decreasing = T),]
#对symbol进行去重
ids=ids[!duplicated(ids$symbol),]
#按照探针id取dat的子集
dat=dat[ids$probe_id,]
dim(dat)

代码语言：javascript复制

## [1] 18821     6

代码语言：javascript复制

#把ids的symbol这一列中的每一行给dat作为dat的行名
rownames(dat)=ids$symbol
# 保存数据集合分组信息
save(dat,group_list,file = 'step1-output.Rdata')

结束语

到这里需要分析的数据已经下载并预处理完成，后面的文章将会基于现在保存的结果进行下一步的主成分分析、差异基因表达分析、GO、KEGG的富集分析等。

love&peace

数据挖掘 go geo 教程

0 人点赞