Affy包 estrogen包

下载安装

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("affy")
BiocManager::install("estrogen")

affy 包为 Bioconductor 之中一个用于数据预处理的包。

affy 包的功能只有一个：读取 affymetrix 的基因表达芯片数据，即CEL格式文件，然后处理成表达矩阵。

先解释以下为什么这个包叫这个名字，这是因为affy其实是一个生产芯片的公司，这个公司做的芯片所产生的的数据肯定不是拿过来就能用的，不同的实验组，不同的PM和MM都是需要处理的。所以他们官方出了一个工具包叫做affy，专门用来处理原始的实验数据。

Affymetrix 的探针（proble）一般是长度为25bp的寡聚核苷酸；探针总是以 perfect match 和 mismatch 成对出现，其信号值称为 PM 和 MM，成对的 perfect match 和 mismatch 有一个共同的 affyID。

CEL文件：信号值和定位信息

CDF文件：探针对在芯片上的定位信息

CEL文件的下载地址（以GSE1438为例），一般在网页的最下方，下载后解压到指定目录即可。

 

 下载到本地并解压之后可以用下面的代码读取 cel 文件

读取过程比较费事，这里使用了 mas5() 对表达数据进行 normalization，当然也可以使用 rma() 函数

> library(affy)
> dir_cels <- "cel_files"  # cel文件存放的路径
> affy_data <- ReadAffy(celfile.path = dir_cels)
> eset.mas5 <- mas5(affy_data)

理论上，处理得到的数据应该和直接从官网下载的表达量数据是一样的 即 XXX_series_matrix.txt.gz 文件

严格讲，这个芯片得到的表达矩阵，是需要过滤的。

setwd('../')
library(affy)
dir_cels='GSE34824_RAW'
data <- ReadAffy(celfile.path=dir_cels)
eset <- rma(data)
calls <- mas5calls(data) # get PMA calls
calls <- exprs(calls)
absent <- rowSums(calls == 'A') # how may samples are each gene 'absent' in all samples
absent <- which (absent == ncol(calls)) # which genes are 'absent' in all samples
rmaFiltered <- eset[-absent,] # filters out the genes 'absent' in all samples

54675 features 经过过滤后，剩下 42482 features

 

基因芯片的原始数据是一个 .cel 文件，当然在我们的 estrogen 包中自带了一小部分 .cel 数据，可以用来学习。

 

参考

http://treeh.cn/?id=22

https://www.jianshu.com/p/ca908850bd84