title: GDAS010-生物學(xué)變異與技術(shù)變異
date: 2019-09-10 12:0:00
type: "tags"
tags:
- Bioconductor
- RNA-seq
categories: - Genomics Data Analysis Series
前言
在后面的章節(jié)中啥箭，我們將介紹基因組學(xué)實(shí)驗(yàn)中的統(tǒng)計(jì)推斷岳遥。我們會(huì)使用前面幾節(jié)中介紹的一些概念刻诊，其中包括t檢驗(yàn)和多重比較似忧。在這一節(jié)中，我們先介紹一個(gè)在基因組學(xué)數(shù)據(jù)分析中一個(gè)極為重要的概念：生物學(xué)變異與技術(shù)變異(biological and technical variability)以及它們之間的區(qū)別骑脱。
通常而言煤辨，我們?cè)谝粋€(gè)種群里所觀察到的生物學(xué)單位沦偎，例如某個(gè)個(gè)體的變異稱為生物學(xué)變異(biological variability)。我們將那些對(duì)同一個(gè)生物學(xué)單位進(jìn)行測(cè)量時(shí)所觀察到的變異稱為技術(shù)性變異(technical variability吹缔，例如同一個(gè)樣本測(cè)3次商佑，這3次就是技術(shù)重復(fù)，3次數(shù)據(jù)的變異稱為技術(shù)變異)厢塘。由于在基因組學(xué)研究中經(jīng)常會(huì)用到新技術(shù)茶没，因此技術(shù)重復(fù)(technical replicates)在很多時(shí)候都會(huì)用于評(píng)估實(shí)驗(yàn)數(shù)據(jù)肌幽。理想情況下，同一樣本的檢測(cè)結(jié)果應(yīng)該是相同的礁叔，但在實(shí)際檢測(cè)中同一樣本的多次檢測(cè)會(huì)出現(xiàn)一定的偏差牍颈，因此使用技術(shù)重復(fù)我們就能我們就能夠評(píng)估新技術(shù)的技術(shù)變異。我還經(jīng)常使用生物學(xué)重復(fù)(biological replicates)與技術(shù)重復(fù)(technical replicates)來分別代指生物學(xué)變異和技術(shù)變異琅关。
我們?cè)谶M(jìn)行統(tǒng)計(jì)推斷來解釋有差異時(shí)煮岁，需要注意不要混淆生物學(xué)變異和技術(shù)變異，因?yàn)樗鼈兯淼囊饬x是不同的涣易。例如画机，當(dāng)我們分析技術(shù)重復(fù)時(shí)，樣本僅僅只有一個(gè)新症，而不是代表某個(gè)總體步氏，例如人類或小鼠。這一節(jié)中徒爹，我們通過設(shè)計(jì)一個(gè)實(shí)驗(yàn)來說明一下技術(shù)重復(fù)和生物重復(fù)荚醒。
合并實(shí)驗(yàn)數(shù)據(jù)
我們將要研究的數(shù)據(jù)集來源于基因表達(dá)微陣列。在這個(gè)實(shí)驗(yàn)中隆嗅，我們從兩個(gè)品系的小鼠中分別隨機(jī)選擇了12只小鼠界阁，并提取了它們的RNA。我們將得到的24個(gè)樣本使用微陣列的手段進(jìn)行檢測(cè)胖喳，引外泡躯，我們還構(gòu)成了兩個(gè)庫，這兩個(gè)庫中的樣本則是分別來源于這兩個(gè)品系的12只小鼠的混合樣本丽焊。
現(xiàn)在我們安裝以下包较剃，這個(gè)包中含有這些數(shù)據(jù)，如下所示：
library(devtools)
install_github("genomicsclass/maPooling")
使用pData函數(shù)我們就能查看實(shí)驗(yàn)設(shè)計(jì)技健。每行表示一個(gè)樣本写穴，每列表示1只小鼠。例如i,j單元格雌贱，中如果是1确垫，則表示RNA數(shù)據(jù)來源于j只小鼠的樣本i。通過行名我們可以知道品系信息（不推薦這種方法來進(jìn)行分組帽芽，你可以向phenoData中添加一個(gè)變量删掀，用于指定小鼠品系信息），數(shù)據(jù)如下所示：
library(Biobase)
library(maPooling)
data(maPooling)
head(pData(maPooling))
## a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12 a14 b2 b3 b5 b6 b8 b9 b10 b11
## a10 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
## a10a11 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0
## a10a11a4 0 0 1 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0
## a11 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
## a12 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
## a12a14 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0
## b12 b13 b14 b15
## a10 0 0 0 0
## a10a11 0 0 0 0
## a10a11a4 0 0 0 0
## a11 0 0 0 0
## a12 0 0 0 0
## a12a14 0 0 0 0
現(xiàn)在我們創(chuàng)建一個(gè)函數(shù)导街，用于說明小鼠的樣本信息披泪，如下所示：
library(rafalib)
mypar()
flipt <- function(m) t(m[nrow(m):1,])
myimage <- function(m,...) {
image(flipt(m),xaxt="n",yaxt="n",...)
}
myimage(as.matrix(pData(maPooling)),col=c("white","black"),
xlab="experiments",
ylab="individuals",
main="phenoData")
需要注意的是，最終我們只對(duì)兩個(gè)品系小鼠的差異基因感興趣搬瑰，也就是數(shù)據(jù)中標(biāo)為0的品系和1的品系款票。我們可以對(duì)合并樣本的技術(shù)重復(fù)進(jìn)行檢驗(yàn)或來源于12只小鼠的數(shù)據(jù)進(jìn)行檢驗(yàn)控硼。我們可以提取出這些合并的樣本，因?yàn)閬碜杂诿總€(gè)品系的所有小鼠都都可以用0或1來表示艾少，因此那些來源于某品系（例如用1表示的這個(gè)品系）代表了這些樣本卡乾，因此我們提取出矩陣行之和為12的基因，這段話不太理解缚够，原文如下：
Note that ultimately we are interested in detecting genes that are differentially expressed between the two strains of mice which we will refer to as strain 0 and 1. We can apply tests to the technical replicates of pooled samples or the data from 12 individual mice. We can identify these pooled samples because all mice from each strain were represented in these samples and thus the sum of the rows of experimental design matrix add up to 12:
如下所示：
data(maPooling)
pd=pData(maPooling)
pooled=which(rowSums(pd)==12)
我們可以從列名確定品系幔妨，如下所示：
factor(as.numeric(grepl("b",names(pooled))))
## [1] 0 0 0 0 1 1 1 1
## Levels: 0 1
如果我們要比較每個(gè)基因在不同品系之間的差異，我們就會(huì)發(fā)現(xiàn)有幾個(gè)基因出現(xiàn)了比較一致的差異谍椅，如下所示：
###look at 2 pre-selected genes for illustration
i=11425;j=11878
pooled_y=exprs(maPooling[,pooled])
pooled_g=factor(as.numeric(grepl("b",names(pooled))))
mypar(1,2)
stripchart(split(pooled_y[i,],pooled_g),vertical=TRUE,method="jitter",col=c(1,2),
main="Gene 1",xlab="Group",pch=15)
stripchart(split(pooled_y[j,],pooled_g),vertical=TRUE,method="jitter",col=c(1,2),
main="Gene 2",xlab="Group",pch=15)
這里要注意一下误堡，如果我們根據(jù)這些值來進(jìn)行t檢驗(yàn)，那么我們就會(huì)得到非常顯著的結(jié)果雏吭，也就是說p值極低锁施，如下所示：
library(genefilter)
pooled_tt=rowttests(pooled_y,pooled_g)
pooled_tt$p.value[i]
## [1] 2.075617e-07
pooled_tt$p.value[j]
## [1] 3.400476e-07
但是，如果我們?cè)俅螐膬蓚€(gè)品系的小鼠中各選擇12個(gè)小鼠的話杖们，這樣的結(jié)論會(huì)成立嗎悉抵？請(qǐng)注意，t檢驗(yàn)的定義包括要比較的總體的標(biāo)準(zhǔn)差摘完。我們這里做的計(jì)算結(jié)果是否能夠反映出總體差異姥饰？（注：這里要回顧一下t檢驗(yàn)的定義，樣本與方差的差異等概念）描焰。
需要注意的是媳否，我們現(xiàn)在重復(fù)計(jì)算的過程就是在重復(fù)原來的實(shí)驗(yàn)流程栅螟。我們?yōu)槊總€(gè)合并的樣本創(chuàng)建了4個(gè)技術(shù)重復(fù)荆秦。或許Gene 1在小鼠品系內(nèi)存在著很高的變異力图，而Gene 2的變異則比較低步绸，但是我們不能看到這一點(diǎn)，因?yàn)樾∈笾g的變異被合并后的數(shù)據(jù)所掩蓋了（注：這里作者只是假設(shè)了Gene 1可能存在著比較高的組內(nèi)變異）吃媒。
我們還有每只小鼠的微陣列數(shù)據(jù)瓤介。對(duì)于每個(gè)品系的小鼠來說，我們有12個(gè)生物學(xué)重復(fù)赘那。我們可以通過查詢1來找到這個(gè)品系刑桑，如下所示：
individuals=which(rowSums(pd)==1)
事實(shí)證明，某些單獨(dú)的小鼠本身就含有技術(shù)重復(fù)募舟，因此現(xiàn)在我們將這種情況剔除掉祠斧，從而用于說明僅用生物學(xué)重復(fù)的分析結(jié)果，如下所示：
##remove replicates
individuals=individuals[-grep("tr",names(individuals))]
y=exprs(maPooling)[,individuals]
g=factor(as.numeric(grepl("b",names(individuals))))
我們可以計(jì)算每個(gè)基因的在某個(gè)品系中的方差拱礁，并與通過技術(shù)重復(fù)獲得的標(biāo)準(zhǔn)差進(jìn)行比較琢锋，如下所示：
technicalsd <- rowSds(pooled_y[,pooled_g==0])
biologicalsd <- rowSds(y[,g==0])
LIM=range(c(technicalsd,biologicalsd))
mypar(1,1)
boxplot(technicalsd,biologicalsd,names=c("technical","biological"),ylab="standard deviation")
從上圖上我們可以看到辕漂，生物學(xué)方差(variance)遠(yuǎn)大于技術(shù)方差。并且對(duì)于方差的變異而言吴超，生物學(xué)方差更大钉嘹。這只是我們前面提到的2個(gè)基因的情況，現(xiàn)在我們展示每只小鼠上測(cè)得的這2個(gè)基因的表達(dá)值鲸阻，如下所示：
mypar(1,2)
stripchart(split(y[i,],g),vertical=TRUE,method="jitter",col=c(1,2),xlab="Gene 1",pch=15)
points(c(1,2),tapply(y[i,],g,mean),pch=4,cex=1.5)
stripchart(split(y[j,],g),vertical=TRUE,method="jitter",col=c(1,2),xlab="Gene 2",pch=15)
points(c(1,2),tapply(y[j,],g,mean),pch=4,cex=1.5)
現(xiàn)在p值就發(fā)生了一點(diǎn)變化跋涣，如下所示：
library(genefilter)
tt=rowttests(y,g)
tt$p.value[i]
## [1] 0.0898726
tt$p.value[j]
## [1] 1.979172e-07
當(dāng)我們?cè)诓煌废档男∈笾g比較這2個(gè)基因的差異時(shí)，我們對(duì)哪個(gè)基因的差異更有信心呢（Gene 1還是Gene 2）赘娄？如果另一位研究人員另外一批隨機(jī)樣本來做實(shí)驗(yàn)仆潮，你認(rèn)識(shí)他會(huì)發(fā)現(xiàn)哪個(gè)基因有差異（Gene 1還是Gene 2）？如果我們希望我們的結(jié)論是關(guān)于小鼠品系的通用結(jié)論（通用結(jié)論這里我的理解就是遣臼，無論放哪個(gè)實(shí)驗(yàn)室做性置，都能做出這個(gè)結(jié)果，得到這個(gè)結(jié)論）揍堰，而不僅僅是我們這次實(shí)驗(yàn)的結(jié)論（也就是說鹏浅，有可能就是只有自己實(shí)驗(yàn)室能做出這個(gè)結(jié)果，那么這個(gè)結(jié)果就不太可信屏歹，有可能是假陽性）隐砸，那么檢測(cè)生物學(xué)的變異則非常重要。
在這一節(jié)中蝙眶，我們使用了兩種分析方法季希。第一種分析方法是將這兩個(gè)品系的生物學(xué)重復(fù)做為總體總體來進(jìn)行分析析（也就是說把某個(gè)品系的12只小鼠數(shù)據(jù)匯總起來分析）。第二種則是將合并后的數(shù)據(jù)來作為總體進(jìn)行分析（我的理解就是把12只小鼠的RNA混合在一起幽纷，作為一個(gè)樣本式塌，然后使用4個(gè)技術(shù)重復(fù)來進(jìn)行分析，這種分析方法其實(shí)檢測(cè)的就是技術(shù)變異）友浸。在科學(xué)研究中峰尝，我們通常關(guān)注的是總體（也就是第一種分析方法）。作為一個(gè)非常實(shí)際的例子收恢，這里需要注意的是武学，如果另外一個(gè)實(shí)驗(yàn)室重復(fù)該實(shí)驗(yàn)，他們會(huì)有另外一組12個(gè)小鼠的樣本伦意，因此關(guān)于我們實(shí)驗(yàn)中對(duì)總體的推斷也會(huì)在另外的實(shí)驗(yàn)室中得到重復(fù)火窒。
An analysis with biological replicates has as a population these two strains of mice. An analysis with technical replicates has as a population the twelve mice we selected and the variability is related to
the measurement technology. In science we typically are concerned with populations. As a very practical example, note that if another lab performs this experiment they will have another set of twelve mice and thus inferences about populations are more likely to be reproducible.
