原筆記地址：http://x2yline.gitee.io/microarray_data_analysis_note/

1. 芯片數(shù)據(jù)結(jié)構(gòu)

芯片數(shù)據(jù)的結(jié)構(gòu)如下

芯片數(shù)據(jù)的結(jié)構(gòu)

1.1 Sample Information

對于Sample Information酝掩，我們至少需要知道樣本的組別妄呕，在一些更加復(fù)雜的使用設(shè)計中，我們可能還需要知道更多關(guān)于樣本的信息如：

• 對于研究癌癥和正常人的實驗設(shè)計中虱朵，我們可能需要記錄每個個體的性別布持，年齡等信息
• 對于動物實驗很多相關(guān)的處理因素可能需要進行記錄

1.2 R中的數(shù)據(jù)框

在R的數(shù)據(jù)框類型數(shù)據(jù)中豌拙，每一行代表一個實驗樣本，每一列代表每一個我們感興趣的因素（包括處理因素）鳖链，dataframe可以同時處理文本因子類型和數(shù)值類型的值因而非常實用

通常我們在excel中創(chuàng)建一個儲存樣本信息的表格【每一行代表一個樣本姆蘸，每一列代表一個因素】墩莫，并存儲為csv或tab分隔的文本文件，再用R讀取為dataframe【read.csv逞敷，read.delim狂秦，read.table函數(shù)】

1.3 AnnotatedDataFrame對象（樣本信息）

Bioconductor的Biobase包通過一個叫做AnnotatedDataFrame對象來存儲樣本信息，這個對象包括：

• 每個樣本的“phenotype”信息推捐，即不同的處理或因素的信息
• 其他meta-information裂问，如對各個處理或不同因素的描述

如我們有8個樣本的芯片數(shù)據(jù)，每個芯片都有200個基因探針：

fake.data <- matrix(rnorm(8*200), ncol=8)

接下來我們需要創(chuàng)建AnnotatedDataFrame對象來描述這8個樣本的表型：

• 包括一個數(shù)據(jù)框用來區(qū)分8個樣本（即8列）的表型信息
• 一個meta數(shù)據(jù)框用了描述表型數(shù)據(jù)框中每一個因素的詳細信息

## 樣本的表型信息
## 注意data.frame函數(shù)會自動把字符串向量轉(zhuǎn)化為因子（factor）
sample.info <- data.frame(
  spl=paste('A', 1:8, sep=''),
  stat=rep(c('cancer', 'healthy'), each=4))

## 對上一個dataframe列名的詳細描述
meta.info <- data.frame (labelDescription =
                           c('Sample Name', 'Cancer Status'))

## 創(chuàng)建AnnotatedDataFrame對象
library(Biobase)
pheno <- new("AnnotatedDataFrame",
             data = sample.info,
             varMetadata = meta.info)

可以使用pData查看AnnotatedDataFrame對象的樣本表型信息

pheno
## An object of class 'AnnotatedDataFrame'
##   rowNames: 1 2 ... 8 (8 total)
##   varLabels: spl stat
##   varMetadata: labelDescription

pData(pheno)
##   spl    stat
## 1  A1  cancer
## 2  A2  cancer
## 3  A3  cancer
## 4  A4  cancer
## 5  A5 healthy
## 6  A6 healthy
## 7  A7 healthy
## 8  A8 healthy

1.4 ExpressionSet對象（表達矩陣與樣本信息的整合）

在bioconductor中牛柒，ExpressionSet對象可以把表達矩陣和樣本信息整合在一起

my.experiments <- new("ExpressionSet",
                      exprs=fake.data,
                      phenoData=pheno)

可以直接查看ExpressionSet的一些基本信息(基因數(shù)堪簿，樣本數(shù)等)

my.experiments
## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 200 features, 8 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: 1 2 ... 8 (8 total)
##   varLabels: spl stat
##   varMetadata: labelDescription
## featureData: none
## experimentData: use 'experimentData(object)'
## Annotation:  

注意：如果使用上述方法創(chuàng)建一個ExpressionSet對象，則一定要保證表達矩陣的列與樣本信息的行相對應(yīng)皮壁，軟件沒有辦法對這些信息進行檢查

1.5 基因信息

我們創(chuàng)建好了一個ExpressionSet對象椭更，但是這個對象卻沒有對基因的注釋及相關(guān)信息，bioconductor的一個瑕疵就是它允許用戶把表達數(shù)據(jù)與基因信息完全分離開（但大多數(shù)人對此都欣然接受）

1.6 ExpressionSet對象的其他可選信息

處理表達數(shù)據(jù)（exprs）和樣本信息（phenoData）外ExpressoinSet對象還允許添加其他的一些可選信息

• annotation：基因注釋信息的環(huán)境（這里有點不太理解蛾魄？虑瀑？？）
• se.exprs：表達矩陣中每個元素對應(yīng)的標準差滴须，因此該數(shù)據(jù)也是一個矩陣
• notes：描述該實驗相關(guān)信息的字符串
• description：MIAMEclass的一個對象舌狗，即Minimun Information About a Microarray Experiment

1.7 內(nèi)置示例數(shù)據(jù)及簡單分析

對于Affymetrix公司的芯片數(shù)據(jù)，bioconductor提供了一個特殊的ExpressionSet數(shù)據(jù)類型扔水，即AffyBatch痛侍，affydata包提供了一個示例數(shù)據(jù)Dilution，Dilution即為一個AffyBatch對象

# BiocInstaller::biocLite("affydata")
library(affydata)
data(Dilution)
Dilution
## AffyBatch object
## size of arrays=640x640 features (30414 kb)
## cdf=HG_U95Av2 (12625 affyids)
## number of samples=4
## number of genes=12625
## annotation=hgu95av2
## notes=

phenoData(Dilution)
## An object of class 'AnnotatedDataFrame'
##   sampleNames: 20A 20B 10A 10B
##   varLabels: liver sn19 scanner
##   varMetadata: labelDescription

pData(Dilution)
##     liver sn19 scanner
## 20A    20    0       1
## 20B    20    0       2
## 10A    10    0       1
## 10B    10    0       2

注：如果是第一次查看Dilution這個數(shù)據(jù)魔市，bioconductor會自動聯(lián)網(wǎng)并更新基因的注釋信息主届，這可能會花一段時間。

首選我們可以看一下芯片的掃描圖像嘹狞，如果出現(xiàn)了異常亮區(qū)或異常暗區(qū)岂膳，則說明質(zhì)量可能存在問題（被刮了一下）如下圖

par(mfrow=c(2, 2))
image(Dilution)
par(mfrow=c(1, 1))

也可以用箱式圖及密度分布圖查看每個芯片的表達量分布（注意表達量自動進行l(wèi)og2處理）誓竿，如果出現(xiàn)了某個芯片與其他芯片的差異很大磅网，說明該芯片質(zhì)量可能存在問題

boxplot(Dilution, col=1:4)
    
hist(Dilution, col=1:4, lty=1)

affy包也提供提取特定探針數(shù)據(jù)的函數(shù)

geneNames(Dilution)[1:3]
## [1] "100_g_at" "1000_at"  "1001_at" 

probeNames(Dilution)[1:48]
##  [1] "100_g_at" "100_g_at" "100_g_at" "100_g_at" "100_g_at" "100_g_at" "100_g_at" "100_g_at"
##  [9] "100_g_at" "100_g_at" "100_g_at" "100_g_at" "100_g_at" "100_g_at" "100_g_at" "100_g_at"
## [17] "1000_at"  "1000_at"  "1000_at"  "1000_at"  "1000_at"  "1000_at"  "1000_at"  "1000_at" 
## [25] "1000_at"  "1000_at"  "1000_at"  "1000_at"  "1000_at"  "1000_at"  "1000_at"  "1000_at" 
## [33] "1001_at"  "1001_at"  "1001_at"  "1001_at"  "1001_at"  "1001_at"  "1001_at"  "1001_at" 
## [41] "1001_at"  "1001_at"  "1001_at"  "1001_at"  "1001_at"  "1001_at"  "1001_at"  "1001_at"

可看出一個序列（一個基因或RNA對應(yīng)一個probset）對應(yīng)16個探針（通常一個probset對應(yīng)11~20種探針，每種探針會重復(fù)數(shù)百次）

# probset函數(shù)截取特定探針的數(shù)據(jù) 
# 產(chǎn)生一個Probset對象組成的list 
# probset對象包含了表達數(shù)據(jù)
probeset(Dilution, c("100_g_at", "34803_at"))
## $`100_g_at`
## ProbeSet object:
##   id=100_g_at
##   pm= 16 probes x  4  chips
## 
## $`34803_at`
## ProbeSet object:
##   id=34803_at
##   pm= 16 probes x  4  chips
## 

ps <- probeset(Dilution, c("100_g_at", "34803_at"))[[2]]

# perfect match
pm(ps)
##          20A   20B   10A   10B
## 235143 435.3 267.0 315.0 192.3
## 121393 206.0 143.8 155.0 118.0
## 99468  349.8 302.8 277.0 171.8
## 196621 567.0 399.0 480.0 322.3
## 256831 160.5 140.5 137.3  97.8
## 351719 474.0 312.8 401.0 240.0
## 361833 561.5 365.5 468.8 319.0
## 111309 742.8 503.0 622.5 346.0
## 407736 713.5 440.0 503.0 303.0
## 407737 736.0 468.0 564.5 327.0
## 238835 808.3 527.5 516.0 321.0
## 173838 463.3 322.3 372.0 228.5
## 320847 240.0 177.0 172.0 136.8
## 195435 129.0  80.0  92.0  75.3
## 327084 193.0 128.3 150.0  96.3
## 12704  120.0  98.0 109.0  80.5

# mis-match
mm(ps)
##          20A   20B   10A   10B
## 235783 230.0 140.3 180.0 111.0
## 122033 109.5  78.0  94.5  67.8
## 100108 162.0 133.8 119.0  85.0
## 197261 223.0 132.3 174.8 102.5
## 257471 164.0 104.0 130.5  76.3
## 352359 371.0 267.0 301.0 187.5
## 362473 487.0 275.5 352.3 200.0
## 111949 439.0 286.0 311.5 183.0
## 408376 393.0 240.3 282.3 154.0
## 408377 655.8 415.8 453.8 282.0
## 239475 885.0 548.3 511.0 338.0
## 174478 215.8 137.0 145.0 100.0
## 321487 134.3  96.0 102.0  88.0
## 196075 113.0  67.3 107.0  62.3
## 327724 132.3  69.3  91.5  67.0
## 13344  107.0 137.0  95.3  64.0

對該基因的探針數(shù)據(jù)（未經(jīng)過log2處理）進行可視化：

plot(c(1,16), c(50, 900), type='n',
     xlab='Probe', ylab='Intensity')
for (i in 1:4) lines(pm(ps)[,i], col=i)
for (i in 1:4) lines(mm(ps)[,i], col=i, lty=2)

## pm-mm的值進行可視化
plot(c(1,16), c(-80, 350), type='n',
     xlab='Probe Pair', ylab='PM - MM')
temp <- pm(ps) - mm(ps)
for (i in 1:4) lines(temp[,i], col=i)

同時affy包也能查看RNA的降解情況筷屡，因為RNA的降解都是從5`向3`開始的涧偷，這16（通常每個probset即每個RNA對應(yīng)這11~20個不同位置的探針）個探針的靶序列也是從RNA的5`向3`方向逐漸靠近的

以雜交的信號為縱坐標，探針的位置作為橫坐標作圖（其中縱坐標進行了縮放和移動以便對不同樣本的斜率進行比較）毙死，可以推測5`端的探針序列對應(yīng)的RNA相對與3`來說密度較低燎潮，信號也越低，因此曲線整體呈上升趨勢

比斜率更重要的是芯片之間的比較扼倘，理想狀況下各芯片的曲線是平行的确封，而降解得很多的樣品其5' 到 3' 線段的斜率會比其他樣品大除呵，該數(shù)據(jù)顯示每個位置的探針降解速度都不太一致（斜率呈分段變化），但樣品之間的降解一致^[1][2]爪喘，通常從5‘到3’端的第5個探針認為是比較穩(wěn)定的

degrade <- AffyRNAdeg(Dilution)
plotAffyRNAdeg(degrade, col=1:4)

下圖為存在異常降解芯片脫穎而出的數(shù)據(jù)：

2. 使用bioconductor讀取Affymetrix芯片數(shù)據(jù)

affy包提供了ReadAffy函數(shù)使創(chuàng)建AffyBatch對像變得簡單颜曾，該函數(shù)還提供了GUI（graphical user interface）界面方便用戶的使用，然而我們需要對分析過程進行文檔化秉剑，為了分析的可重復(fù)性我們應(yīng)該避免使用GUI界面

由CEL文件創(chuàng)建AffyBatch對象

library(affy)
# There is a function list.celfiles that 
# makes it easy to get a list of the CEL files
fns <- list.celfiles("data/CEL_file_dir", 
                     full.names = TRUE)
fns[1:2]
## [1] "CEL_file_dir/CL2001011101AA.CEL" "CEL_file_dir/CL2001011102AA.CEL"

# read CEL data
my.data <- ReadAffy(filenames=fns) 

構(gòu)建樣本信息的數(shù)據(jù)AnnotatedDataFrame對象

其中樣本的信息文件（data_file_list.txt）格式如下：

scan.name   sample_name type
CL2001011101AA  ALL_1   A
CL2001011104AA  ALL_2   A
CL2001011105AA  ALL_3   A   

使用函數(shù)Biobase包的read.AnnotatedDataFrame讀取注釋文件泛豪，設(shè)置第二列為行名：

adf <- read.AnnotatedDataFrame(file.path("data", "/data_file_list.txt"),
                               header = TRUE, sep = "\t", row.names = 2)

adf
## An object of class 'AnnotatedDataFrame'
##   rowNames: ALL_1 ALL_2 ... AML_28 (72 total)
##   varLabels: scan.name type
##   varMetadata: labelDescription

# 對列名的注釋信息（可選）
varMetadata(adf) <- data.frame(labelDescription=
                                 c("file name", "disease classification"))

phenoData(my.data) <- adf

也可以再構(gòu)建一個description，即MIAME對象

miame <- new("MIAME", name = "someone",
             lab = "your lab name",
             title = "Gene expression data ...",
             pubMedIds = c("your pubmed ID"))

description(my.data) <- miame

description(my.data)
## Experiment data
##   Experimenter name: someone 
##   Laboratory: your lab name 
##   Contact information:  
##   Title: Gene expression data ... 
##   URL:  
##   PMIDs: your pubmed ID 
##   No abstract available.

另一種方法為之間使用read.affybatch函數(shù)之間構(gòu)建AffyBatch對象

abatch <- read.affybatch(filenames = fns, 
                         phenoData = adf, 
                         description = miame)

3. 數(shù)據(jù)預(yù)處理

bioconductor把low-level processing（即把探針級別【一個基因?qū)?yīng)多個探針】的數(shù)據(jù)轉(zhuǎn)為為基因級別的數(shù)據(jù)侦鹏，從而使一個probset對應(yīng)一個表達量值）分為了4個步驟诡曙，每個步驟都可以選擇不同的算法，這些方法的不同極有可能會導(dǎo)致后面high-level analyses的差異

這四個步驟為：

• Background correction
• Normalization (on the level of features = probes)
• PM-correction（包括只使用PM略水，丟棄掉MM价卤，或者把PM-MM作為校正的PM等）
• Summarization（把一個probset中的多個probe對應(yīng)的多值總結(jié)為1個值）

3.1 背景信號的校正

背景校正的目的是為了消除人為的或芯片本身存在的噪聲所帶來的誤差

bgcorrect.methods()
## [1] "bg.correct" "mas"        "none"       "rma"   

• mas：使用 MAS 5.0 算法
• none：不做校正
• rma：使用 RMA 算法

MAS 5.0 算法

MAS 5.0 的背景校正方法為把每個芯片（每個CEL文件）分為16個區(qū)，在每個區(qū)中最最暗的2%作為每個區(qū)的背景水平渊涝，在對該芯片的每個探針用這16個背景值的加權(quán)平均進行校正荠雕，每個區(qū)的權(quán)重取決于該探針與芯片中這16個區(qū)域中心的距離。

RMA算法

RMA算法只對PM進行背景校正驶赏，MM的值不做校正（放棄MM炸卑，因為幾乎有1/3的MM值是大于對應(yīng)的PM的），該模型的假設(shè)為PM的值是由兩個部分組成的即S=X+Y煤傍，其中S為實際檢查到的信號盖文，X代表有探針雜交產(chǎn)生的信號服從指數(shù)分布Exp(α)，Y代表背景信號服從正態(tài)分布N(μ, δ^2)蚯姆，Y只取大于0的部分五续，并且X與Y是獨立的

最后校正過后的X值=E(X|S = s)，其中s為該探針的信號強度龄恋，E(X|S = s)表示在檢測到的信號強度為s的情況下X的期望值疙驾，具體推導(dǎo)可以看大牛的博士論文：http://bmbolstad.com/Dissertation/Bolstad_2004_Dissertation.pdf，暫時沒有時間來研究郭毕，先擱在這里它碎。

還有一個gcrma的算法，考慮到不同序列g(shù)c含量與序列親和力的關(guān)系進而估計不同背景信號大小显押，需要安裝gcrma包bgc.gcrma = gcrma(Data)扳肛，該命令包括了背景校正，標準化和summarization這幾個步驟乘碑。

不同算法的結(jié)果

library(affydata)
data(Dilution)
d.mas <- bg.correct(Dilution[, 1], "mas")
d.rma <- bg.correct(Dilution[, 1], "rma")
bg.with.mas <- pm(Dilution[, 1]) - pm(d.mas)
bg.with.rma <- pm(Dilution[, 1]) - pm(d.rma)

par(mfrow=c(1,2))
hist(bg.with.mas)
hist(bg.with.rma)
par(mfrow=c(1,1))

summary(data.frame(bg.with.mas, bg.with.rma))
##       X20A           X20A.1     
##  Min.   :74.53   Min.   : 72.4  
##  1st Qu.:93.14   1st Qu.:113.7  
##  Median :94.35   Median :114.9  
##  Mean   :94.27   Mean   :112.1  
##  3rd Qu.:95.80   3rd Qu.:114.9  
##  Max.   :97.67   Max.   :114.9 

在該示例數(shù)據(jù)中挖息，RMA算法估計的背景值比MAS 5.0估計的背景值大，而且RMA算法的背景值比較恒定

# 兩種算法的差值分布
hist(bg.with.rma - bg.with.mas)

差值大部分為20兽肤，可以由芯片的信號值分布看一看20究竟有多大

tmp <- data.frame(pm(Dilution[, 1]), mm(Dilution[, 1]))
colnames(tmp) <- c("PM", "MM")
summary(tmp)
##        PM                MM         
##  Min.   :   76.0   Min.   :   77.3  
##  1st Qu.:  137.0   1st Qu.:  120.3  
##  Median :  225.0   Median :  164.5  
##  Mean   :  507.3   Mean   :  323.5  
##  3rd Qu.:  489.0   3rd Qu.:  313.0  
##  Max.   :23356.3   Max.   :17565.3 
boxplot(tmp, ylim=c(0, 1000))

可以看出20越占大部分信號值的1/10套腹，不同算法的差值還是比較大的

3.2 summarization=Quantication

暫時先不管normalization和PM correction這一步绪抛，summarization的目的就是把所有同一probset（即同一基因mRNA）的PM和MM探針對轉(zhuǎn)化為一個數(shù)值，這個數(shù)值是對目標基因表達量的最佳估計电禀，其方法有：

express.summary.stat.methods()
## [1] "avgdiff"      "liwong"       "mas"          "medianpolish" "playerout"   

使用不同的summarization方法對probeset進行定量睦疫，expresso函數(shù)把所有預(yù)處理步驟都整合在了一起

tempfun <- function(method) {
  expresso(Dilution, bg.correct = FALSE,
           normalize = FALSE, pmcorrect.method = "pmonly",
           summary.method = method)
}

ad <- tempfun("avgdiff")
al <- tempfun("liwong")
am <- tempfun("mas")
ar <- tempfun("medianpolish")

MA plot

MA plot即M-versus-A plot，在芯片數(shù)據(jù)處理出現(xiàn)之前也稱為Bland-Altman plot鞭呕，是由發(fā)明者名字命名的蛤育，而MA plot是對M與A作圖而得名，M是minus的縮寫葫松，代表兩個值之差瓦糕，A是add的縮寫，代表兩個值之和咕娄。有研究者也把MA plot稱為Ratio-Intensity (RI) plots。

MA plot的作用是為了展示兩個值幾乎處處相等的變量（x和y）之間的關(guān)系珊擂，為了展示兩個變量之間的變化關(guān)系圣勒，大多數(shù)人的思維都是把x與y分別作為橫軸和縱軸進行繪圖，如果y=x摧扇，則該圖呈45度角的直線（如下圖中左邊圖的藍色直線）圣贸，可以通過查看點形成的直線偏離預(yù)期直線的多少來衡量系統(tǒng)偏差，然而該圖存在以下幾個缺點：

1. 人的視覺對水平線比更敏感
2. 不同坐標軸的刻度可能會使預(yù)期參考直線偏離45度
3. 很難從直觀上衡量偏離一條線性的大小

MA plot的處理方法是把該直線順時針旋轉(zhuǎn)45度扛稽，把參考對角線變?yōu)橹本€吁峻，具體做法是把(x+y)/2作為橫軸，(y-x)作為縱軸在张，則參考直線變?yōu)橐粭l水平線用含，如下方右圖，這樣可以很清楚的在視覺上展示兩個相等的變量之間偏離參考值的大小帮匾，即存在的系統(tǒng)誤差的大小

bioconductor中的MA plot

affy包的函數(shù)mva.pairs可以快速地做出MA plot（也可以使用MAplot函數(shù)）啄骇，接下來我們使用MA plot來比較不同summarize方法之間的差異，mva.pairs函數(shù)的輸入?yún)?shù)就是各變量形成的一個dataframe

## ar進行了log2處理
temp <- data.frame(exprs(ad)[, 1], 
                   exprs(al)[,1], 
                   exprs(am)[, 1],
                   2^exprs(ar)[, 1])
dimnames(temp)[[2]] <- c("Mas4", "dChip",
                         "Mas5", "RMA")
affy::mva.pairs(temp)

上圖中的差異也有可能是由于我們沒有在summarization之前進行校正瘟斜，接下來我們使用RMA方法進行背景校正缸夹，使用百分位數(shù)標準化方法，并且只用PM值進行數(shù)據(jù)定量匯總

tempfun <- function(method) {
  expresso(Dilution,
         bgcorrect.method = "rma",
         normalize.method = "quantiles",
         pmcorrect.method = "pmonly",
         summary.method = method)
}

ad <- tempfun("avgdiff")
al <- tempfun("liwong")
am <- tempfun("mas")
ar <- tempfun("medianpolish")

## ar進行了log2處理
temp <- data.frame(exprs(ad)[, 1], 
                   exprs(al)[,1], 
                   exprs(am)[, 1],
                   2^exprs(ar)[, 1])
dimnames(temp)[[2]] <- c("Mas4", "dChip",
                         "Mas5", "RMA")
mva.pairs(temp)

summarization的詳細過程

每個算法都包括了四個步驟的不同處理方式（即一個算法包含4個小算法）哼转，如MAS 4.0 包括MAS明未，constant，subtractmm壹蔓，avgdiff

不同算法的組合如下表（摘抄自：NAR文獻Comparison of algorithms for the analysis of Affymetrix microarray data as evaluated by co-expression of genes in known operons）：

MAS 4.0 算法（MicroArray Suite 4.0）

MAS4.0代表一個完整的數(shù)據(jù)處理算法，對背景校正猫态，標準化佣蓉，PM校正披摄，probeSet定量4個完整步驟都有其特定的計算過程：

• 背景校正算法也叫做MAS，即使用2%的信號作為背景如前所述
• 標準化算法為constant勇凭，即把所有芯片的中位數(shù)縮放致一個相等的常數(shù)
• PM校正為PM-MM疚膊，即subtractmm
• summarization使用叫做avgdiff的算法，即把PM-MM的離群值（偏離均值大于3倍標準差）刪除后計算平均值為基因的表達量值

對avgdiff的評價：

• 對所有的探針都一視同仁（沒有加權(quán)重）
• 對PM-MM的值只是簡單的相加
• 可能會把一些有趣的探針給刪除掉（離群值）
• 可能會產(chǎn)生負值（大約有1/3的MM大于PM）
• 不能有效利用其他芯片的數(shù)據(jù)（沒有建立一個模型）

MAS 5.0 算法（MicroArray Suite 5.0）

mas5.0的4個步驟中虾标，除了標準化為constant外寓盗，其余均名為mas

PM校正方法名為mas：即把PM-CT作為校正的PM，CT為change threshold是一個調(diào)整后的MM璧函，以保證CT值小于PM值傀蚌，CT值算法為

• 如果MM小于PM則CT=MM
• 如果MM大于PM，則CT校正為一個比PM較小的數(shù)蘸吓，如下圖所示

summarization方法也稱為mas：即Signal = anti-log of Tukeybiweight(log (PM –CT))

點評：mas 仍然沒有用到不同芯片的數(shù)據(jù)即建立一個模型

開發(fā)了mas 5.0后善炫，affymetrix公司不再投入大量精力開發(fā)算法，但提供了一套可以檢驗算法的數(shù)據(jù)The Affy Latin Square Experiment.用這個數(shù)據(jù)可以證明mas 5.0比mas 4.0（avgdiff）表現(xiàn)更優(yōu)越

dChip（一個芯片數(shù)據(jù)分析軟件库继，哈佛）算法也較Li-Wong算法

Li-Wong算法【由現(xiàn)北大研究員李程與現(xiàn)斯坦福教授王永雄開發(fā)】包括三個小部分箩艺，不包括背景校正過程，標準化為invariantset（以后再了解）宪萄，PM校正方法為只使用PM值或PP和MM值都使用（可自定義）艺谆，summarization為一個統(tǒng)計學(xué)模型稱為Li-Wong算法竭宰，所有芯片的數(shù)據(jù)都得到了很好的利用祟昭，但是如果芯片數(shù)目比較少時可信度較差

[圖片上傳失敗...(image-b9fc36-1525698429039)]

[圖片上傳失敗...(image-b1b7e2-1525698429039)]

其模型最終PM-MM值為探針親和力$\varphi_{j}$與表達量$\theta_{i}$的乘積，i代表樣本已旧，j代表探針聊记，$v_{i}$代表非特異性雜交的基線水平撒妈，文章地址為：http://www.pnas.org/content/98/1/31

RMA（Robust Multichip Analysis）算法

RMA算法是2003年由Rafael A. Irizarry（Bekerley）等人提出的，前面已經(jīng)大致介紹了背景校正的方法排监，標準化使用的是一種叫做median polish的方法（以后再介紹）狰右，該方法完全丟棄了MM值（他們給出的理由是因為太多MM信號值大于PM值的情況，因此考慮MM的值會引入很大的變異）【他們的想法或許是正確的】

像dChip一樣舆床，該算法的summarization方法是建立在一個統(tǒng)計學(xué)模型之上的：

[圖片上傳失敗...(image-a12c9f-1525698429039)]

其中芯片編號為i棋蚌，j代表探針

模型的參數(shù)由多個芯片的數(shù)據(jù)進行擬合得到，與dChip不同的是$\epsilon$的值是經(jīng)過對數(shù)處理的挨队，這充分考慮到了值越大的探針其變異也越大的特點

PDNN模型（考慮雜交的熱力學(xué)因素）

以上算法都是從數(shù)學(xué)層面上進行處理谷暮，而沒有試圖解釋不同探針之間的信號差異來源，以及是什么決定了非特異性結(jié)合盛垦。通常這些都與探針的序列湿弦，及雜交的熱力學(xué)性質(zhì)相關(guān)

Li Zhang提出了Position-Dependent Nearest Neighbor(PDNN) 模型(Nat Biotech, 2003; 21:818).不像dChip和RMA，PDNN模型參數(shù)的估計只需要單個芯片的數(shù)據(jù)（很大程度上是由于參數(shù)數(shù)目比較少）腾夯。該模型考慮序列對雜交的影響颊埃，把信號強度值與序列聯(lián)系起來蔬充。

該模型的參數(shù)有：

1. 堿基對在序列中的位置k
2. 與臨近堿基對的作用，除了知道k班利，還需要知道k-1和k+1

該模型的實現(xiàn)需要單獨安裝一個包affypdnn饥漫，在載入這個包之后，會自動載入相關(guān)函數(shù)罗标，并更新可選的模型

# registering new summary method 'pdnn'.
# registering new pmcorrect method 'pdnn' and 'pdnnpredict'.

express.summary.stat.methods()
## [1] "avgdiff"      "liwong"       "mas"          "medianpolish" "playerout"    "pdnn" 

需要注意的是PDNN模型沒有按照標準的expresso4個步驟處理數(shù)據(jù)庸队，該模型直接從百分位數(shù)校正開始，可以只使用PM或PM和MM同時進行擬合闯割，背景的估計整合到了能量和表達參數(shù)里面彻消。

需要使用expressopdnn這個expresso的變體實現(xiàn)PDNN模型的運算。

4. 質(zhì)量控制（QC）

4.1 簡單質(zhì)控

即前面的簡單分析部分

第一步是看圖形纽谒，有沒有異常的區(qū)域

接著可以看看表達值的分布情況证膨，看有沒有異常的芯片

image(Dilution)
boxplot(Dilution, col = 2:5)
hist(Dilution, col = 2:5, lty = 1)

4.2 simpleaffy包進行質(zhì)控

simpleaffy質(zhì)控得到的結(jié)果都是按照MAS5.0算法產(chǎn)生的，如校正因子鼓黔，背景信號強度等都是MAS 5.0算法的結(jié)果

require(simpleaffy)
Dil.qc <- qc(Dilution)

首先我們可以看看每個芯片的背景值是否有可比性央勒，差別非常大的異常值

# 各芯片背景均值
sort(avbg(Dil.qc))
##      10B      20B      10A      20A 
## 54.25830 63.63855 80.09436 94.25323 

# 各芯片背景最大值
sort(maxbg(Dil.qc))
##      10B      20B      10A      20A 
## 57.62283 68.18998 83.24646 97.66280 

# 各芯片背景最小值（MAS 5.0算法）
sort(minbg(Dil.qc))
##      10B      20B      10A      20A 
## 49.22574 60.01397 77.32196 89.52555 

芯片標準化時，標準的Affymetrix流程為乘上一個校正因子（scaling factors）使各芯片的中位數(shù)相同澳化，Affymetrix要求可比較的芯片之間的校正因子不能超過3

# 標準校正流程中的校正因子
sfs(Dil.qc)
## [1] 0.8934013 1.2653627 1.1448430 1.8454067

max(sfs(Dil.qc))/min(sfs(Dil.qc))
## [1] 2.065597

接下來是percent of present值（即有多少基因按照Detection calls算法是確實存在表達且被檢測到的）是否有可比性崔步，如果有離群值或低于30，高于60的值說明可能存在質(zhì)量問題

percent.present(Dil.qc)
## 20A.present 20B.present 10A.present 10B.present 
##    48.79208    49.82178    49.37822    49.75842 

Affymetrix設(shè)計了一些參考管家基因（通常是GAPDH或beta-actin）的3‘缎谷，middle和5’的探針井濒，其中actin的3‘/5'的值應(yīng)該是小于3，GAPDH的3'/5'值小于1（5’端可能存在部分降解列林，可與RNA降解圖結(jié)合起來看）瑞你，否則可能存在問題

ratios(Dil.qc)
##     actin3/actin5 actin3/actinM gapdh3/gapdh5 gapdh3/gapdhM
## 20A     0.6961423     0.1273385     0.4429746   -0.06024147
## 20B     0.7208418     0.1796231     0.3529890   -0.01366293
## 10A     0.8712069     0.2112914     0.4326566    0.42375270
## 10B     0.9313709     0.2725534     0.5726650    0.11258237

也可以之間對qc結(jié)果進行作圖查看

中間的黑色實線代表fold change為0的參考線，黑色虛線代表上調(diào)或下調(diào)基因的參考線即fold change絕對值超過3

藍色區(qū)域代表各樣本的scale factor理論上的范圍跨度（其確定方式為以各芯片的scale factor均值為中心先兩側(cè)各延伸1.5個fold change單位希痴，即跨度剛好為3）

橫線可是藍色或紅色（藍色代表芯片的scale factor是可接受的者甲，紅色代表有芯片的scale factor超過了理論范圍不具有可比性，下圖為藍色說明各芯片的scale factor比較相似具有可比性）砌创；

actin的3‘/5'的值應(yīng)該是小于3虏缸，GAPDH的3'/5'值小于1，任何超過該范圍的標記為紅色否則為藍色嫩实；

左邊的數(shù)值上面的是percent of present（差異超過10%則標記為紅色）刽辙，下面的是average background（差異超過20個單位則標記為紅色，但是這個值可能需要調(diào)整）

plot(Dil.qc)

4.3 RNA降解圖使用affy包查看甲献，可參考之前的部分

4.4 依賴統(tǒng)計學(xué)模型的質(zhì)控包affyPLM

affyPLM是一個與RMA模型相似的probe-level model (PLM)從而進行質(zhì)控

library("affyPLM")
pset <- fitPLM(Dilution)

模型擬合的殘差可以用image函數(shù)來可視化顯示宰缤，可以很清楚地判斷是否有異常殘差值的出現(xiàn)，網(wǎng)址http://www.plmimagegallery.bmbolstad.com/收集了很多特征性的圖片，而Dilution這個數(shù)據(jù)集基本正常撵溃，通常該圖片沒有瑕疵的情況很少疚鲤，小的瑕疵一般都可以接受锥累，通過RMA這個穩(wěn)健的算法基本可以把誤差消除掉缘挑，而當污點比較大時，我們需要再看其他的質(zhì)控數(shù)據(jù)決定是否應(yīng)該舍棄該芯片桶略。

par(mfrow = c(2, 2))

# 第三張芯片的原始圖片
image(Dilution[, 3])

# 第三張芯片的校正后值圖片（默認為type=“weight"语淘，值越大顏色越淡）
image(pset, type = "weights", which = 3)

# 第三張芯片的殘差的圖片（紅色漸變色表示正值越大，白色表示接近0际歼，藍色漸變色表示負值越大）
image(pset, type = "resids", which = 3)

# 第三張芯片的殘差的圖片（紅色表示所有正值惶翻，白色表示接近0，藍色表示所有負值）
image(pset, type = "sign.resids", which = 3)

下面這個就是有問題的圖片鹅心，除原始圖片不太清楚外吕粗，其他圖片都明顯顯示有人為誤差，即一個“環(huán)”旭愧，就需要看下一步的QC數(shù)據(jù)是否有較大異常決定去留

還可以用箱式圖查看M值（相對于表達量為中等的那個芯片）分布也叫做Relative Log Expression的分布颅筋。 如果有芯片的中位數(shù)不是0或跨度明顯大于其他芯片，則可能存在誤差输枯，這個兩個圖都是一模一樣的议泵，不過名稱不同而已

require(RColorBrewer)
cols <- brewer.pal(8, "Set1")
par(mfrow=c(1,2))
Mbox(pset, main="M value for Dilution dataset", col=cols)
RLE(pset,
    main="RLE for Dilution dataset", 
    ylim=c(-0.7, 1.1),
    col=cols)
par(mfrow=c(1,1))

RLE(pset, type = "stats")
##                20A          20B         10A         10B
## median 0.002359773 -0.007542337 0.003051789 0.001902185
## IQR    0.092560253  0.089580148 0.076335304 0.093406833

RLE(pset, type = "values")[1:4, ]
##                  20A          20B          10A         10B
## 100_g_at  0.08993939 -0.074287948  0.026305114 -0.02630511
## 1000_at   0.11551276 -0.007420133  0.007420133 -0.04293690
## 1001_at   0.08565396 -0.071191808 -0.043903380  0.04390338
## 1002_f_at 0.11275165 -0.045783998 -0.013287804  0.01328780

接下來就是看Normalized Unscaled Standard Error (NUSE)值的分布是否異常，與M值（RLE值）的思想差不多桃熄，NUSE是每個基因的標準差與該基因中所有芯片中的標準差中位數(shù)的比值先口，因此芯片具有可比性的條件是其所有基因的標準差幾乎相同，即NUSE的中位數(shù)為1瞳收，且分布大致相當碉京。

par(mfrow=c(1,2))
boxplot(pset, main="NUSE for Dilution dataset", col=cols)
NUSE(pset, main="NUSE for Dilution dataset", 
     col=cols, 
     ylim=c(0.88, 1.23))
par(mfrow=c(1,1))

我們再來看一個異常的RLE和NUSE圖，其中第二個芯片與image一樣都偏離了正常值

本文為MD Anderson芯片分析課程的lecture 7筆記螟深，最后的異常值圖片參考子bioconductor書籍Bioinformatics and Computational Biology Solutions Using R and Bioconductor

其他參考：
http://repository.countway.harvard.edu/xmlui/bitstream/handle/10473/4713/hummel_affy_normalization.pdf
http://www.math.usu.edu/~jrstevens/stat5570/Bowles.pdf

1. http://arrayanalysis.org/main.html#RD ?

2. http://www.reibang.com/p/859a7a7d040f ?