【实战】
数据来源:GSE98793
芯片平台:GPL570([HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array)
登录GEO数据库,在检索框中输入GSE98793,打开数据界面。
下载最下方的原始数据。
下载后大概有857MB大小。注意,使用浏览器下载时可能会出现数据下载不完整的问题(即下载的后的文件大小远小于原始文件大小)。不完整的数据无法使用。 如下图所示:
将该文件解压到Rworkplace的工作路径下,大体如下图所示,以“.CEL.gz”为文件后缀。
至此,前期准备工作完成。接着读取CEL数据,并对数据进行预处理。
首先要考虑用什么R包去读取CEL数据。
参考以下信息,我们的芯片平台前面也提到过为——GPL570 [HG-U133_Plus_2]。
即可以采用affy包提取数据。
Affymetrix表达谱芯片数据读取的方法分3种:
1、使用affy包读取。(HGU95/HGU133芯片)
2、使用oligo包读取。(Whole Transcriptome 芯片/ NimbleGen 芯片/ SNP芯片等)
3、使用simpleaffy包读取。(HGU95/HGU133芯片)
我也参考了文献中的芯片数据的处理方法,作者也是使用了affy进行CEL数据读取工作。
###############################读取文件夹下的CEL数据########################################
setwd("E://Rworkplace")
library(affy)
Data1 <- ReadAffy()
上述代码读取了Rworkplace下的所有CEL文件。
由于ReadAffy函数在读取数据的时候,将文件名默认为样本名,读取后的数据的样本名非常长。
于是做处理,简化,将最前方的GSM-作为新的样本名。
使用stringr包,将原始样本名的那一列,以“_”为分隔符,转化为10列。
然后去这10列中的第一列为新的样本名。代码如下所示,
library(stringr)
raw.names<-sampleNames(Data1)
new<-str_split_fixed(raw.names, "_", 10)
new.names<-new[,1]
sampleNames(Data1)<-new.names
pData(Data1)
新的样本名转化为如下,
接着使用RMA算法对原始数据进行预处理。
tum.rma <- rma(Data1)
rma算法主要执行如下步骤:
提取表达矩阵。
##############################提取表达矩阵################################################
tum_exp <- exprs(tum.rma)
发现,行号为探针名,而不是我们想要的基因名。
于是,下一步,将探针名转换为基因名。
对于GPL570平台而言,其注释的数据包为hgu133plus2.db。(对于平台与注释包的一一对应关系,可以参考链接)
下载平台注释包。
##############################对矩阵的probe进行注释########################################
##--加载平台注释包(version:3.2.3)
library(BiocInstaller)
biocLite("hgu133plus2.db")
library(hgu133plus2.db)
如下图所示,其中 hgu133plus2SYMBOL 是基因名与探针之间的意义对应关系。
使用toTable()函数来提取基因名与探针名之间的一一对应关系。
ids=toTable(hgu133plus2SYMBOL)
将数据集中的探针名与注释包中的基因名一一对应。
其中需要用到的一个中介是注释包中的探针名。
所以基本思路是,在注释包中找到与数据集中一样的探针名,然后选择与探针名对应的基因名。
其中数据集中的有些探针名在注释包中是没有注释的信息的,所以需要删去数据集中的这部分的探针。
##--将注释与当下数据集中所用的probe一一对应
length(unique(ids$symbol)) #注释包的探针数目:20188
table(rownames(tum_exp) %in% ids$probe_id)
#flase:12738;true:41937
data<-tum_exp[rownames(tum_exp) %in% ids$probe_id,]#去除了原始数据中未经注释的探针
dim(data)
ids<-ids[match(rownames(data),ids$probe_id),]
tmp = by(data,ids$symbol,function(x) rownames(x)[which.max(rowMeans(x))] )
probes<-as.character(tmp)
data<-data[rownames(data) %in% probes,]
ids<-ids[ids$probe_id %in% probes,]
row.names(data)<-ids$symbol
最后,将整理好的数据集,存入文件中,以备下次调用。
###########################将初始化的数据存入文档中###################################
write.table(data, file="GSE98793.txt")
完整代码如下:
###############################读取文件夹下的CEL数据########################################
setwd("E://Rworkplace")
library(affy)
Data1 <- ReadAffy()
#############################修改样本名####################################################
library(stringr)
raw.names<-sampleNames(Data1)
new<-str_split_fixed(raw.names, "_", 10)
new.names<-new[,1]
sampleNames(Data1)<-new.names
pData(Data1)
##############################用RMA算法进行预处理#########################################
tum.rma <- rma(Data1)
##############################提取表达矩阵################################################
tum_exp <- exprs(tum.rma)
##############################对矩阵的probe进行注释########################################
##--加载平台注释包(version:3.2.3)
library(BiocInstaller)
biocLite("hgu133plus2.db")
library(hgu133plus2.db)
ids=toTable(hgu133plus2SYMBOL)
##--将注释与当下数据集中所用的probe一一对应
length(unique(ids$symbol)) #注释包的探针数目:20188
table(rownames(tum_exp) %in% ids$probe_id)#flase:12738;true:41937
data<-tum_exp[rownames(tum_exp) %in% ids$probe_id,]#去除了原始数据中未经注释的探针
dim(data)
ids<-ids[match(rownames(data),ids$probe_id),]
tmp = by(data,ids$symbol,function(x) rownames(x)[which.max(rowMeans(x))] )
probes<-as.character(tmp)
data<-data[rownames(data) %in% probes,]
ids<-ids[ids$probe_id %in% probes,]
row.names(data)<-ids$symbol
###########################将初始化的数据存入文档中###################################
write.table(data, file="GSE98793.txt")