【结合文献】——Affymatrix芯片数据预处理

【实战】

数据来源：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")