首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >geo数据挖掘-2

geo数据挖掘-2

作者头像
火星娃统计
发布2020-09-15 15:30:20
1.2K0
发布2020-09-15 15:30:20
举报
文章被收录于专栏:火星娃统计火星娃统计

geo数据挖掘-2

sunqi
2020/7/11

1.概述

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

2.数据下载

2.1 获得表达数据‘

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)
## [1] "list"
#长度
length(gset)
## [1] 1
# 因为只有一个平台,所以只有1个列表元素
class(gset[[1]])
## [1] "ExpressionSet"
## attr(,"package")
## [1] "Biobase"
# 列表的第一个元素为表达矩阵和biobase
gset
## $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
sprintf("显示下载的文件有6个样本,22397个位点 /n GPL6244")
## [1] "显示下载的文件有6个样本,22397个位点 /n GPL6244"
#获取列表元素,
a=gset[[1]]

#exprs函数获取表达矩阵
dat=exprs(a)
dim(dat)#查看维度
## [1] 33297     6
# GPL6244
# 查看前5行
head(dat[,1:5])
##         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
# 可以看到列名为样本好,行名为探针名
# 绘制箱式图
boxplot(dat,las=2)
#查看临床信息,包含6个患者的34个信息
pd=pData(a)

## 选择需要的临床信息
library(stringr)
# 通过空格分隔title,获得分组信息
group_list=str_split(pd$title,' ',simplify = T)[,4]
table(group_list)
## group_list
##     Control Vemurafenib
##           3           3
# 3个病例和3个对照

2.2 获得平台信息

# 查看平台信息探针信息
# 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)
## [1] 33297     6
dat=dat[ids$probe_id,]
dim(dat)
## [1] 19814     6
#对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)
## [1] 18821     6
#把ids的symbol这一列中的每一行给dat作为dat的行名
rownames(dat)=ids$symbol
# 保存数据集合分组信息
save(dat,group_list,file = 'step1-output.Rdata')

结束语

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

love&peace

本文参与 腾讯云自媒体分享计划,分享自微信公众号。
原始发表:2020-07-12,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 火星娃统计 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体分享计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • geo数据挖掘-2
    • 1.概述
      • 2.数据下载
        • 2.1 获得表达数据‘
        • 2.2 获得平台信息
      • 结束语
      领券
      问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档