Affymetrix的表达量芯片的cel文件如何处理

生信技能树

发布于 2024-05-27 17:16:16

7590

如果是常规的geo表达量芯片数据集代码，比如illumina的芯片，我们汇总了系列代码：

https://www.jianguoyun.com/p/DdqkaeUQ1pC6BhixiLAFIAA
表达量芯片是非常适合锻炼大家的r编程基础的，新的一年从这3个gse数据集开始吧：
2015-GSE67936-AML-illumina
2016-GSE65409-AML-illumina
2019-GSE114868-AML-hta2.0

而且绝大部分表达量芯片并不需要从原始数据开始，比如affymetrix的芯片，一般来说就是读取作者给出来的表达量矩阵文件即可，比如 GSE30122_series_matrix.txt.gz 文件是 7.0M ，可以看到它在线链接是有规律的：https://ftp.ncbi.nlm.nih.gov/geo/series/GSE30nnn/GSE30122/matrix/GSE30122_series_matrix.txt.gz

读取作者给出来的表达量矩阵文件的标准代码如下所示：

library(AnnoProbe)
library(GEOquery)  
getOption('timeout')
options(timeout=10000) 
list.files() 
if(T){gset <- geoChina(gse_number)}
if(F){ gset = getGEO("GSE30122", destdir = '.', getGPL = F) } 
a=gset[[1]] 
dat=exprs(a) #a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数
dim(dat)#看一下dat这个矩阵的维度
dat[1:4,1:4] #查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列
pd = pData(a)
head(pd)

boxplot(dat[,1:4],las=2)

可以很明显的看到作者给出来的表达量矩阵里面的表达量是被zscore的：

表达量是被zscore的

这个时候就需要下载这个项目的raw文件了，因为是affymetrix芯片，所以绝大部分是cel格式的文件，在线链接仍然是有规律的：https://ftp.ncbi.nlm.nih.gov/geo/series/GSE30nnn/GSE30122/suppl/GSE30122_RAW.tar

同样的，很简单的读取它：

library(oligo) 
celFiles <- list.celfiles('../cel_files/',listGzipped = T,
                          full.name=TRUE)
celFiles 
options(BioC_mirror="https://mirrors.westlake.edu.cn/bioconductor")
exon_data <- oligo::read.celfiles( celFiles ) 
dim(exprs(exon_data)) 
exprs(exon_data)[1:4,1:4]
### 标准化(一步完成背景校正、分位数校正标准化和RMA (robust multichip average) 表达整合)
exon_data_rma <- oligo::rma(exon_data  ) 
exp_probe <- Biobase::exprs(exon_data_rma)
exp_probe[1:4,1:4]

dat=exp_probe
dim(dat)#看一下dat这个矩阵的维度
dat[1:4,1:4] #查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列

如下所示看起来就很正常的表达量矩阵分布范围啦：