Affymetrix 芯片数据处理流程及原理简介

版权所有，转载必究！

本文作者:任红雷, 联系邮箱: renhongleiz@126.com

本文整理自个人汇报的PPT

首先，本文将结合代码讲解常用的Affymetrix 芯片数据处理流程

以GSE11787为例,数据为关于副猪嗜血杆菌对猪炎症发生影响的数据，请下载数据，并将其解压到任意的文件夹，以供之后使用：

http://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE11787&format=file

首先介绍一下关于基因表达芯片的背景。（最后附上流程代码）

1.提纲

提纲

2.中心法则(可自行百度)

中心法则

3.DNA 芯片是什么，作用如何？

芯片的含义与作用

5.Affymetrix 芯片

Affymetrix 芯片

6.芯片上的信息

芯片上的信息

7.PM和MM

PM和MM

8.芯片扫描结果

芯片扫描结果及其工作原理

9..CEL文件格式

.CEL文件格式

10.intensity(信号强度)的实际含义

intensity(信号强度)的实际含义

11..CEL文件以及.CDF文件

.CEL文件以及.CDF文件

12.数据处理流程

数据处理流程

13.拓展学习书籍列表

（本教程内容整理自Bioconductor Case Studies）

1.A (very) short introduction to R

https://cran.r-project.org/doc/contrib/Torfs+Brauer-Short-R-Intro.pdf

2.An Introduction to R

https://cran.r-project.org/doc/manuals/R-intro.pdf

3.Bioconductor Case Studies

lots of bioconductor workflows with source codes

链接: http://pan.baidu.com/s/1qYBtvIk

密码: eq85

Affymetrix 芯片数据处理流程的源代码

1.Read data

##1.ReadData##
#add bioconductor source

source("https://bioconductor.org/biocLite.R")

#install affy package from bioconductor

biocLite("affy")

#load affy package

library("affy")

#set the directory that you are working with,this can be replaced by your own CEL files path

setwd("/Users/Ren/Documents/Rcode/GSE11787_RAW")

#set the .CEL files path,this can be replaced by your own CEL files path

celpath="/Users/Ren/Documents/Rcode/GSE11787_RAW/"

#read .CEL files from directory of celfile.path

data = ReadAffy(celfile.path =celpath)

#replace the sample names of the data by trim off the ".CEL" suffix

sampleNames(data) = sub("\\.CEL$", "", sampleNames(data))

#read the phenoData of the CEL file.The "type.csv" file is consitituded with .csv format after you comprehensed the samples grouping on the GEO: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE11787

samplestype <- read.csv(paste(celpath,"type.csv",sep = ""),header = TRUE,row.names = 1,na.strings = "NA",sep = ',',stringsAsFactors = F,as.is = !default.stringsAsFactors())

#replace the rownames of the samplestype by its sampleID

rownames(samplestype) = samplestype$sampleID

#match the rownames of your phenoData to the samplenames of experiments data,which will return the matched index

mt = match(rownames(samplestype), sampleNames(data))

#give a completed column description of the phenoData

vmd = data.frame(labelDescription = c("Sample ID", "Samples type: Haemophilus Parasuis infected or control"))

#make the matched sampletype rows and varMetadata to your data's phenoData

phenoData(data) = new("AnnotatedDataFrame", data = samplestype[mt, ], varMetadata = vmd)

#erase samples which has no type

data = data[,!is.na(data$type)]

其中type.csv的文件形式如下:

type.csv

2.Quality Control

##2.Quality Control##

#install the essential packages from Bioconductor

source("https://bioconductor.org/biocLite.R")

biocLite(c("affyQCReport","simpleaffy"))

#load the pakages of affyQCReport and simpleaffy

library("affyQCReport")

library("simpleaffy")

#Execute the Quality Control

saqc = qc(data)

#plot the quality control result

plot(saqc)

QC 之后的结果,如下图：

Quality Control

3.Data Normalization

##3.Data Normalization##
#packaged affyPLM provides another set of diagnostics that can be used to help assess array quality

source("https://bioconductor.org/biocLite.R")

biocLite(c("affyPLM"))

library("affyPLM")

#fit the basic probe-level model

dataPLM = fitPLM(data)

#plot normalized unscaled standard error (NUSE)

boxplot(dataPLM, main="NUSE", outline = FALSE, col="lightblue", las=3, whisklty=0, staplelty=0)

#plot relative log expression (RLE)

Mbox(dataPLM, main="RLE", outline = FALSE, col="mistyrose", las=3, whisklty=0, staplelty=0)

#if some sample was significantly elevated or more spread out,you can remove it as bad arrays

#here we don not have bad array ,the code next line is just a example to remove the sample "GSM298263",to execute this ,you can remove the sharp at the head of next line

#badArray = data[,-match(c("GSM298263"), sampleNames(data)) ]

normalized unscaled standard error (NUSE) 的各样本结果。

NUSE

relative log expression (RLE) 的各样本结果。

RLE

4.Data preprocessing

##4.Data preprocessing##
#load affy package

library("affy")

#rma preprocessing:which concluded 1.background correction,2.between array normalization 3.reporter summarization

datarma = rma(data)

source("https://bioconductor.org/biocLite.R")

#install the porcines' db for get its specices annotation of the probe set

biocLite(c("porcine.db"))

#load the package of porcine.db

library("porcine.db")

#to filter out probe sets with no Entrez Gene identifiers and Affymetrix control probes

datafiltered = nsFilter(datarma, remove.dupEntrez=FALSE, var.cutof =0.5)$eset


#t-tests for rows of a matrix, intended to be speed efficient

#The function computes the t-statistic, the difference of the group means between the two disease groups, and the corresponding p-value by using the t-distribution.

datatt = rowttests(datafiltered, "type")

#draw the volcano plot using x-axis:dm(difference of the group means),y-axis:group log(p-value) 

lod = -log10(datatt$p.value)

plot(datatt$dm, lod, pch=".", xlab="log-ratio",ylab=expression(-log[10]~p))

#line indicates an untransformed p-value of 0.01, so points above it will be significant

abline(h=2)

source("https://bioconductor.org/biocLite.R")

#install the porcine's db for get its specices annotation of the probe set

biocLite(c("limma"))

library("limma")

#make the design of the experiment as a matrix with its type

design = model.matrix(~datafiltered$type)

datalim = lmFit(datafiltered, design)

#When there are few replicates, the variance is not well estimated and the t-statistic can perform poorly. 

#It can be solved by a improved method called empirical Bayes to give a sensible estimation of variance and t-statistic

#When sample sizes are moderate or large, say ten or more in each group, there is generally no advantage (but also no disadvantage) to using the Bayesian approach.

dataeb = eBayes(datalim)

火山图的结果(Volcano plot),黑线之上的基因代表p值显著的基因

Volcano plot

5.GO analysis

##GO analysis##

source("https://bioconductor.org/biocLite.R")

#install the porcine's db for get its specices annotation of the probe set

biocLite(c("topGO"))

library(topGO)

#multiple testing problem corrected by multtest package

#Alternatively,the function topTable in the limma package provides multiple testing adjustment methods, including Benjamini and Hochberg’s false discovery rate (FDR), simple Bonferroni correction, and several others.

#for more information you can access:http://www.jianshu.com/p/9e97e9b351fd,and http://www.jianshu.com/p/a262cf3d18b9

#get the top 10 differential expressed gene by multiple testing correction of Benjamini-Hochberg (FDR) method

tabofallgene = topTable(dataeb, coef=2, adjust.method="BH", n=length(datatt[,1]))

#get the p-value of all genes named with the probe set name

geneList=setNames(tabofallgene[,5],rownames(tabofallgene))

annotation_db_name=paste(annotation(datafiltered),".db",sep="")

#load the annotation db package

library(package = annotation_db_name, character.only = TRUE)

#method to extract top differential expressed gene

topDiffGenes <- function(allScore) {
  return(allScore < 0.01)
}

#get the count of the differential expressed gene set

sum(topDiffGenes(geneList))

#new a topGOdata type object and use biology process(a gene ontology category),all gene is geneList, differential expressed gene is  topDiffGenes

sampleGOdata <- new("topGOdata", description = "Simple session", ontology = "BP",allGenes = geneList, geneSel = topDiffGenes,nodeSize = 10,annot = annFUN.db, affyLib = annotation_db_name)

#Performing the enrichment tests

#Fisher’s exact test which is based on gene counts

resultFisher <- runTest(sampleGOdata, algorithm = "classic", statistic = "fisher")

#Kolmogorov-Smirnov like test which computes enrichment based on gene scores.

resultKS <- runTest(sampleGOdata, algorithm = "classic", statistic = "ks")

resultKS.elim <- runTest(sampleGOdata, algorithm = "elim", statistic = "ks")

#make all GO enrichment result to a table called allRes

allRes <- GenTable(sampleGOdata, classicFisher = resultFisher,classicKS = resultKS, elimKS = resultKS.elim,orderBy = "elimKS", ranksOf = "classicFisher", topNodes = 10)

#draw the GO tree graph,with a depth of 10

tree_depth=10

showSigOfNodes(sampleGOdata, score(resultKS.elim), firstSigNodes = tree_depth, useInfo = 'all')

GO terms的DAG图：

GO

6.可选步骤.绘制样本集的热图(heatmap)和level plot

##prestep:1.read data from .CEL files and 2.Quality control

##generate sample's level plot and heatmap##

#calculate the distance of a n*n matrix,and set the diagonal to zero

dd = dist2(log2(exprs(data)),diagonal=0)

#Hierarchical clustering of dd

dd.row <- as.dendrogram(hclust(as.dist(dd)))

#x or y axis sample order index in dd matrix

row.ord <- order.dendrogram(dd.row)

source("https://bioconductor.org/biocLite.R")

biocLite(c("latticeExtra"))

library("latticeExtra")

#add lengend

legend = list(top=list(fun=dendrogramGrob, args=list(x=dd.row, side="top")))

#give a level plot

lp = levelplot(dd[row.ord, row.ord],scales=list(x=list(rot=90)), xlab="", ylab="", legend=legend)

lp

install.packages(c("gplots","RColorBrewer"))

library(gplots)

library(RColorBrewer)

# creates a own color palette from red to green

my_palette <- colorRampPalette(c( "green","yellow", "red"))(n = 299)

gplots:::heatmap.2(dd, Rowv =FALSE,Colv=dd.row, col = my_palette, tracecol=NA,density.info="none",cexRow=0.4,cexCol=0.4)

热图如下:

heatmap

level plot如下:

levelplot