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)來(lái)分別代指生物學(xué)變異和技術(shù)變異。
我們?cè)谶M(jìn)行統(tǒng)計(jì)推斷來(lái)解釋有差異時(shí)，需要注意不要混淆生物學(xué)變異和技術(shù)變異，因?yàn)樗鼈兯淼囊饬x是不同的。例如，當(dāng)我們分析技術(shù)重復(fù)時(shí)，樣本僅僅只有一個(gè)，而不是代表某個(gè)總體，例如人類或小鼠。這一節(jié)中，我們通過(guò)設(shè)計(jì)一個(gè)實(shí)驗(yàn)來(lái)說(shuō)明一下技術(shù)重復(fù)和生物重復(fù)。
合并實(shí)驗(yàn)數(shù)據(jù)
我們將要研究的數(shù)據(jù)集來(lái)源于基因表達(dá)微陣列。在這個(gè)實(shí)驗(yàn)中，我們從兩個(gè)品系的小鼠中分別隨機(jī)選擇了12只小鼠，并提取了它們的RNA。我們將得到的24個(gè)樣本使用微陣列的手段進(jìn)行檢測(cè)，引外，我們還構(gòu)成了兩個(gè)庫(kù)，這兩個(gè)庫(kù)中的樣本則是分別來(lái)源于這兩個(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ù)來(lái)源于j只小鼠的樣本i。通過(guò)行名我們可以知道品系信息（不推薦這種方法來(lá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ù)，用于說(shuō)明小鼠的樣本信息，如下所示：
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)或來(lái)源于12只小鼠的數(shù)據(jù)進(jìn)行檢驗(yàn)。我們可以提取出這些合并的樣本，因?yàn)閬?lái)自于每個(gè)品系的所有小鼠都都可以用0或1來(lái)表示，因此那些來(lái)源于某品系（例如用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ù)這些值來(lái)進(jìn)行t檢驗(yàn)，那么我們就會(huì)得到非常顯著的結(jié)果，也就是說(shuō)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ì)算的過(guò)程就是在重復(fù)原來(lái)的實(shí)驗(yàn)流程。我們?yōu)槊總€(gè)合并的樣本創(chuàng)建了4個(gè)技術(shù)重復(fù)?；蛟SGene 1在小鼠品系內(nèi)存在著很高的變異，而Gene 2的變異則比較低，但是我們不能看到這一點(diǎn)，因?yàn)樾∈笾g的變異被合并后的數(shù)據(jù)所掩蓋了（注：這里作者只是假設(shè)了Gene 1可能存在著比較高的組內(nèi)變異）。
我們還有每只小鼠的微陣列數(shù)據(jù)。對(duì)于每個(gè)品系的小鼠來(lái)說(shuō)，我們有12個(gè)生物學(xué)重復(fù)。我們可以通過(guò)查詢1來(lái)找到這個(gè)品系，如下所示：
individuals=which(rowSums(pd)==1)
事實(shí)證明，某些單獨(dú)的小鼠本身就含有技術(shù)重復(fù)，因此現(xiàn)在我們將這種情況剔除掉，從而用于說(shuō)明僅用生物學(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è)品系中的方差，并與通過(guò)技術(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ī)樣本來(lái)做實(shí)驗(yàn)，你認(rèn)識(shí)他會(huì)發(fā)現(xiàn)哪個(gè)基因有差異（Gene 1還是Gene 2）？如果我們希望我們的結(jié)論是關(guān)于小鼠品系的通用結(jié)論（通用結(jié)論這里我的理解就是，無(wú)論放哪個(gè)實(shí)驗(yàn)室做，都能做出這個(gè)結(jié)果，得到這個(gè)結(jié)論），而不僅僅是我們這次實(shí)驗(yàn)的結(jié)論（也就是說(shuō)，有可能就是只有自己實(shí)驗(yàn)室能做出這個(gè)結(jié)果，那么這個(gè)結(jié)果就不太可信，有可能是假陽(yáng)性），那么檢測(cè)生物學(xué)的變異則非常重要。
在這一節(jié)中，我們使用了兩種分析方法。第一種分析方法是將這兩個(gè)品系的生物學(xué)重復(fù)做為總體總體來(lái)進(jìn)行分析析（也就是說(shuō)把某個(gè)品系的12只小鼠數(shù)據(jù)匯總起來(lái)分析）。第二種則是將合并后的數(shù)據(jù)來(lái)作為總體進(jìn)行分析（我的理解就是把12只小鼠的RNA混合在一起，作為一個(gè)樣本，然后使用4個(gè)技術(shù)重復(fù)來(lái)進(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.
