回顾基因表达芯片分析(数据处理)

去年的这个时候,我参加了生信技能树Jimmy办的线上R语言处理表达芯片教程,这也是我进入生信领域,第一个学习的完整流程。那时对生信的一些基础概念还是处于模糊的状态,只是略懂了点R语言,勉强的听完了整个教程。

后来似乎也再也没接触过芯片数据了,对于芯片数据的概念还停留在表达矩阵那块,以致有次在找公共数据时没找到表达矩阵,只找到了芯片的raw data,却因无法分析而无能为力(现在看来是由于在学习教程时并没有完全吃透

有幸不久之前看到一篇朋友推荐的芯片教程Differential gene expression analysis I (microarray data),看了后感觉受益不少,然后再从头看了一遍Jimmy的教程,粗略了解Affymetrix平台芯片的如何进行前期预处理以及获取表达矩阵

基础知识

一组探针(或探针组probe set)来自于一个基因,通常是由11对探针组成,每对探针是由匹配探针(PM)和错配探针(MM)组成,叫做探针对(Probe Pair),样本RNA则叫做靶标。芯片中MM探针作用是检测非特异杂交信号,理论上,MM只有非特异杂交信号,而不会有特异性杂交。MM的信号值永远小于其对应的PM信号值,因此可以做一个PM-MM或PM/MM即可去除背景噪声的影响

数据下载

以芯片教程中的GSE27114作为实例进行后续操作。在GSE27114里有两个study,一个是比较旧的Affymetrix平台数据(GPL1261):Affymetrix Mouse Genome 430 2.0 Array(Affymetrix 3’ end arrays),另一个是相对较新的平台数据(GPL6246):Affymetrix Mouse Gene 1.0 ST Array(Affymetrix arrays Whole Transcript Gene Arrays),下面先以GPL1261为例:

在数据下载之前还需要先了解自己需要哪一个芯片数据,也就是说是需要芯片原始数据(CEL文件)还是表达矩阵数据(matrix)

如果是只想获得芯片数据已经做完标准化等预处理过的表达矩阵,那么可以使用GEOquery包的getGEO函数,表达矩阵文件的压缩包也会下载到R的工作目录下

library(GEOquery)
gset <- getGEO("GSE27114", destdir = ".")

由于GSE27114里有2个study,所以我先看看有哪两个study

names(gset)
##[1] "GSE27114-GPL1261_series_matrix.txt.gz" "GSE27114-GPL6246_series_matrix.txt.gz"

由于下面会以GPL1261为例,所以我先取出第一个study,可以看到gset_one的数据是储存在ExpressionSet对象中,这个对象里包含了N多的信息,可以通过slotNames列出可以查看的信息,然后根据列出来的函数一一查看,比如assayData用于保存表达数据信息,experimentData用于保存实验设计信息

gset_one <- gset[[1]]
class(gset_one)
##[1] "ExpressionSet"
##attr(,"package")
##[1] "Biobase"

接下来就可以通过exprs函数获取表达矩阵,有45101个探针,6个样本

exprSet <- exprs(gset_one)
dim(exprSet)
##[1] 45101     6

还可以通过pData获取这个study的更多详细信息,比如样本信息、分组信息、rawdata信息,比如这个例子就有32列信息。。。

pdata <- pData(gset_one)
dim(pdata)
##[1]  6 32

当然如果只想获取表达矩阵,并且其他信息并不想知晓(或者已经在GEO网站上查到),那么通过GEO网站手动下载表达矩阵也是可以的:点击GSE27114网页下方的Series Matrix File(s)即可进入ftp下载界面ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE27nnn/GSE27114/matrix/,而且这个ftp的URL可以根据你的GSE号自行进行修改,再配合脚本就可以批量下载N个GSE样本的表达矩阵,很好用。


如果是想获取GSE27114的RAW DATA,那么可以选择使用oligo包(Affymetrix 平台),主要用于读取Affymetrix平台的基因表达数据,处理为表达矩阵;如果是Illumina平台,教程则推荐使用beadarray或lumi

首先通过download.file函数下载raw data,至于下载URL如何获取有两个途径:

  1. 通过上述的pdata里查看,有一列是各个样本的rawdata下载地址,然后写个循环即可批量下载数据

    pdata$supplementary_file
    
  2. 通过GSE27114网站的最下方GSE27114_RAW.tar的,点击http即可下载所有文件的压缩包(网速慢或者不稳定可以首选这个,至少下载一个算一个,慢慢来。。。)

数据读入及标准化

将所有样本的raw data文件放置在R工作目录下(Microarray文件夹),然后用list.celfiles函数查看所有文件,再用read.celfiles函数读入R中,并将上述的pdata信息也一并写入

library(oligo)
file_CELs <- list.celfiles("Microarray", listGzipped = TRUE, full.name = TRUE)
rawAffy <- read.celfiles(filenames = file_CELs, phenoData = phenoData(gset_one), sampleNames = rownames(pData(gset_one)))

注:读入文件时会检查R中是否安装了该数据对应平台的注释包,比如这个例子需要这个pd.mouse430.2包

粗略看下,rawAffy数据中包含了芯片上所有探针的强度值信息

head(exprs(rawAffy), n = 2)

我们可以拿这个数据来进行质控,但是一般质控会在芯片公司那边就做了,所以我们拿公共数据的用的话,一般也不用进行质控了。所以我们接下来就对rawAffy数据进行标准化处理,因为每个样本芯片的绝对光密度是不一样的,不做处理的话无法将各个样本间进行数据分析,同时也是为了消除测量间的非实验差异

芯片常用的标准化方法也有好几种(如:RMA、MAS5),这里就不一一尝试了,就拿最常用的RMA(Robust Multichip Average algorithm),引用芯片教程,RMA标准化过程主要分为3步:

  1. Background correction (removes array auto-fluorescence ) 背景矫正
  2. Quantile normalization (makes all intensity distributions identical) 标准化
  3. Probeset summarization (calculates one representative value per probeset) 汇总

由于rawAffy是ExpressionFeatureSet对象,则rma函数默认是probeset level class(rawAffy) ##1 “ExpressionFeatureSet” ##attr(,”package”) ##1 “oligoClasses” eset <- rma(rawAffy)

我们可以看看标准化后的结果

rma_boxplot.png rma_hist.png

最后我们需要得到表达矩阵,也是用exprs函数,从结果上来看,跟之前getGEO得到的表达矩阵是一样的(说明用的是同一种标准化方法?以及相同的QC处理?)

exprSet <- exprs(eset)
dim(exprSet)
##[1] 45101     6

参考资料:
Differential gene expression analysis I (microarray data)

本文出自于http://www.bioinfo-scrounger.com转载请注明出处