转载自:R语言GEO数据挖掘01-数据下载及提取表达矩阵 – 简书
- 这一节的内容包括应用
GEOquery
包下载芯片数据,提取表达矩阵,提取metadata信息。- 解决一个探针对应多个基因的问题
GEO数据下载-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")##获取帮助
选择一个数据集GSE7765演示分析
数据集基本情况
芯片平台GPL96,GPL97; sample数12个
下载表达矩阵
gse <- getGEO("GSE7765", GSEMatrix = TRUE)
show(gse)
GSE7765中包括两个平台,两个数据集
提取表达矩阵及metadata
class(gse)
str(gse)
a<-gse[[1]]
b<-gse[[2]]
class(gse[[1]])##ExpressionSet
##提取第一个数据集的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
平台注释信息处理
芯片数据分析中很重要的内容即平台信息处理,获取相应的平台
这里我们选择GPL96
if(F){load(file="expma.Rdata")}
GPL="GPL96"##下载平台注释
gpl %
Table() %>% ##转换为data.frame格式
save(file = "GPL96_annot.Rdata")
if(F){load(file = "GPL96_annot.Rdata")}
head(gpl)
colnames(gpl)
##取出注释信息
probe%
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
到这里我们发现一个情况:1个探针存在对应多个基因的情况,///
我们的解决思路是有几种,第一种是直接将存在///的信息去掉
第二种是将///的数据拆开来,然后再把有重复的删去
先按第一种解决,这种解决比较简单,我的需要是只要找出///
probe<-probe[!grepl(" /// ",probe$`Gene Symbol`),]
dim(probe)##数量减少到20878
第二种方法是,先拆解再删除
拆解时需要将///分割开,再与ID相连
值得注意的是这种方法其实与第一种是有区别的,这种方法仍然保留了探针对应多个基因中的
一种情况,所有得出的probe注释要多,这里不纠结这个内容
if(F){
library(tidyverse)
probe2%
unlist()##得到的是个list
###
head(probe2)
length(probe2)##展开后得到24807个探针及对应关系
probe2% separate(value,c("ID","GeneName"),sep = "\|")
dim(probe2)##增加到24807行
## 找出重复ID,两个table的妙用
table(table(probe2$ID))##探针找出对应一个基因的有20878个,与grepl法得出的结果相同
## 下一步的目的即筛选出对应一个基因的探针
test2% count(ID) %>% filter(n==1) %>% ## count计数有点类似于table
inner_join(probe2,by="ID") %>% ## 内连接只保留x,y中观测相同的变量
select(-n)##remove "n" column
dim(test2)
head(test2)
probe2<-test2##将最终得到的结果赋值给probe2
}
博主,如果作者留的Supplementary file 里格式是tsv的,上面的series matrix files里面也是空数据集,这我应该怎么做呢,求助!(具体GEO accession number是GSE118089)
tsv也是制表符分割,打开看一下,如果只是数据文件,直接用,根据样本信息分组,进行差异分析,如果包含样本信息,把样本信息提取出来,只保留表达矩阵。然后一样的差异分析
有问题可以进交流群讨论
Jingle——科研互助群 – 王进的个人网站
https://www.jingege.wang/jingle_science/
博主你好,我问一下GEO上面的高通量测序数据这样下载只有样本信息没有表达矩阵怎么办
GEO accession number多少 我看看
GSE138963,GSE60873只要是高通量测序的数据都提不出来
是的 测序数据原始数据需要在SRA上下载 GEOquery包只是用芯片数据
对于测序数据 作者大部分已经上传了表达矩阵 直接下载即可,比如GSE138963 :https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE138963
底下Supplementary file 提供了GSE138963_all_genes_tpm.txt.gz处理好的表达矩阵
我也知道在那里,就是想快速批量下一些,我要下的有点多,看来还是不行只有去一个一个下了,感谢博主回复
我的数据集显示,有上千的,但是还有负数, 迷茫 不知道有没有标准化
根据你的描述是没有的 GEO一般下载的芯片数据 其实只用log2转换一下
你好,麻烦问一下如果表达矩阵里,范围是-3~20这是已经log过的吗?如果是这样的话2的20次方=1048576这样也太大了
您好,有负值的话是log2转化过的,芯片的原始数据是没有负值的,0-1的数据log转化之后就是负值。
具体的标准化信息在你下载的数据里也可以看到
好的,谢谢进哥哥。麻烦再问一下怎么查看下载数据的处理方式呢?非常感谢
博主我用GEO的时候发现他们的箱式图不太整齐。我不确定是不是还要进行normalizeBetweenArrays
您好 是需要的 刚刚转载了一个关于normalize的帖子,设计一个数据集内和数据集间GSE83521/GSE89143去除批次效应
博主我看你新发的贴子, 如果我不是两个GEO数据集合并, 就一个数据集使用,箱式图不是很齐 我需要拉平吗? 我试了拉平和不拉平的差异很大
当然需要的,可能原因是没有进行log转化,导致数据中位数差得很大,如果发现数值有上千的 先log2转化 然后再看箱式图齐不齐 不齐在进行normalizeBetweenArrays。