幾種常用的差異分析方法簡(jiǎn)介
如今在生物學(xué)研究中导盅，差異分析越來(lái)越普遍缠黍，也有許多做差異分析的方法可供選擇。但是在實(shí)際應(yīng)用中妻味，大多數(shù)人不知道該使用哪種方法來(lái)處理自己的數(shù)據(jù)正压，所以今天我就來(lái)介紹下目前幾種常用的差異分析方法及其適用場(chǎng)景。
1.方差分析责球、T檢驗(yàn)焦履、卡方檢驗(yàn)、秩和檢驗(yàn)
- 其實(shí)核心的區(qū)別在于：數(shù)據(jù)類型不一樣雏逾。如果是定類和定類嘉裤，此時(shí)應(yīng)該使用卡方分析；如果是定類和定量栖博，此時(shí)應(yīng)該使用方差或者T檢驗(yàn)屑宠。
- 方差和T檢驗(yàn)的區(qū)別在于，對(duì)于T檢驗(yàn)的X來(lái)講仇让，其只能為2個(gè)類別比如男和女典奉。如果X為3個(gè)類別比如本科以下，本科丧叽，本科以上卫玖；此時(shí)只能使用方差分析。
進(jìn)一步細(xì)分
T檢驗(yàn)
t檢驗(yàn)(student t檢驗(yàn))是應(yīng)用t分布的特征踊淳，將t作為檢驗(yàn)的統(tǒng)計(jì)量來(lái)進(jìn)行統(tǒng)計(jì)推斷方法假瞬。它對(duì)樣本要求較小(例如n＜30)。
主要用途：
- 樣本均數(shù)與總體均數(shù)的差異比較
- 兩樣本均數(shù)的差異比較
單樣本t檢驗(yàn)
單樣本t檢驗(yàn)主要用于判斷樣本均數(shù)與總體均數(shù)是否存在顯著差異。
適用條件
已知一個(gè)總體均數(shù)
已知一個(gè)樣本均數(shù)及該樣本標(biāo)準(zhǔn)差
樣本正態(tài)分布或近似正態(tài)總體
實(shí)際應(yīng)用中脱茉，當(dāng)數(shù)據(jù)量足夠大時(shí)剪芥，對(duì)樣本正態(tài)分布要求不再嚴(yán)格。只要數(shù)據(jù)分布不是嚴(yán)重偏態(tài)芦劣，一般來(lái)說單樣本t檢驗(yàn)都是適用的粗俱。
R語(yǔ)言中可以用t.test函數(shù)進(jìn)行t檢驗(yàn)
從某小學(xué)六年級(jí)抽取10名學(xué)生说榆，其身高(單位：cm),是否認(rèn)為該學(xué)校六年級(jí)平均身高130cm?
#生成數(shù)據(jù)
x <- c(130,120,130,110,130,135,129,124,130,134)
#t檢驗(yàn)
t.test(x,mu = 130)
One Sample t-test
data: x
t = -1.1884, df = 9, p-value =
0.2651
alternative hypothesis: true mean is not equal to 130
95 percent confidence interval:
121.8702 132.5298
sample estimates:
mean of x
127.2
#結(jié)果顯示虚吟，P=0.2651＞0.05。在統(tǒng)計(jì)學(xué)上說明樣本均數(shù)與總體均數(shù)沒有差別签财。
獨(dú)立樣本t檢驗(yàn)
獨(dú)立樣本t檢驗(yàn)主要檢驗(yàn)兩個(gè)樣本均數(shù)及其所代表的總體之間差異是否顯著串慰。
適用條件
- 獨(dú)立性，各觀察值之間相關(guān)獨(dú)立
- 正態(tài)性唱蒸，各樣本均來(lái)自正態(tài)分布的總體
- 方差齊性邦鲫，各樣本所在總體的方差相等
方差齊性可以用car包leveneTest函數(shù)檢驗(yàn)
leveneTest(y=,group =)
其中，y是兩組樣本組成的數(shù)據(jù)神汹，group是兩組樣本的分組情況庆捺。
方差齊性檢驗(yàn)之后，才可進(jìn)行獨(dú)立樣本t檢驗(yàn)屁魏。
用t.test(A,B,var.equal=TRUE滔以，paired=FALSE)
A、B為數(shù)據(jù)集氓拼，var.equal=TRUE為方差齊性你画。paired=FALSE非配對(duì)樣本。
示例：
(虛構(gòu))有兩組學(xué)生(每組10人)桃漾，一組采用傳統(tǒng)教育坏匪，一組采用素質(zhì)教育。一學(xué)期后撬统，兩組學(xué)生語(yǔ)文成績(jī)(滿分100)如下适滓。問兩組學(xué)生成績(jī)之間差別是否顯著。
A <- c(85,84,95,73,77,65,85,93,90,91)
B <- c(87,96,77,80,79,96,93,82,84,86)
#方差齊性檢驗(yàn)
#合并數(shù)據(jù)
y <- c(A,B)
#數(shù)據(jù)分組標(biāo)簽
group=as.factor(c(rep(1,10),rep(2,10)))
#載入car包
library(car)
#方差齊性檢驗(yàn)
leveneTest(y=y,group = group)
#結(jié)果顯示恋追，P=0.5505＞0.05凭迹。說明方差齊性。
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 0.3703 0.5505
18
#獨(dú)立樣本t檢驗(yàn)
t.test(A,B,paired = FALSE)
#結(jié)果顯示P=0.5639几于。說明兩者沒有區(qū)別蕊苗。
Welch Two Sample t-test
data: A and B
t = -0.589, df = 16.463, p-value = 0.5639
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-10.100024 5.700024
sample estimates:
mean of x mean of y
83.8 86.0
配對(duì)樣本t檢驗(yàn)
配對(duì)樣本t檢驗(yàn)同樣檢驗(yàn)兩個(gè)樣本均數(shù)及其所代表的總體之間差異是否顯著。
獨(dú)立樣本t檢驗(yàn)與配對(duì)樣本t檢驗(yàn)同屬于雙樣本t檢驗(yàn)沿彭，不同點(diǎn)在于配對(duì)樣本t檢驗(yàn)要求兩個(gè)樣本之間存在某些配對(duì)關(guān)系朽砰。
常見配對(duì)關(guān)系：
同一樣本兩種不同處理方法的檢驗(yàn)結(jié)果
同一樣本前后時(shí)間點(diǎn)的檢驗(yàn)結(jié)果
適用條件
正態(tài)性
示例
有20名女性分為10對(duì)，試吃兩種藥。經(jīng)過一段時(shí)間后瞧柔，藥效如下漆弄。問兩種藥是否有區(qū)別
#生成數(shù)據(jù)
drug1 <- c(4.4,5,5.8,4.6,4.9,4.8,6,5.9,4.3,5.1)
drug2 <- c(6.2,5.2,5.5,5,4.4,5.4,5,6.4,5.8,6.2)
#配對(duì)樣本t檢驗(yàn)
t.test(drug1,drug2,paired = TRUE)
#結(jié)果顯示，P=0.1575＞0.05,不能說兩者存在顯著差別造锅。
Paired t-test
data: drug1 and drug2
t = -1.5417, df = 9, p-value = 0.1575
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-1.0609306 0.2009306
sample estimates:
mean of the differences
-0.43
方差檢驗(yàn)
方差分析（analysis of variance ,ANOVA）就是通過檢驗(yàn)各總體的均值是否相等來(lái)判斷分類型自變量對(duì)數(shù)值型因變量是否有顯著影響撼唾。
示例
我們使用的是R里內(nèi)置的“npk”數(shù)據(jù)集，該數(shù)據(jù)集由24行和5列數(shù)據(jù)組成哥蔚，第一列代表區(qū)組（共6個(gè)）倒谷，N、P和K分別代表氮糙箍、磷和鉀元素的使用情況渤愁，yield代表豌豆產(chǎn)量，該數(shù)據(jù)集主要是用來(lái)研究不同肥料對(duì)豌豆產(chǎn)量的影響深夯。
fit <- aov(yield ~ N, data=npk)
fit <- aov(yield ~ N + block, data=npk)
卡方檢驗(yàn)
卡方檢驗(yàn)就是統(tǒng)計(jì)樣本的實(shí)際觀測(cè)值與理論推斷值之間的偏離程度抖格，實(shí)際觀測(cè)值與理論推斷值之間的偏離程度就決定卡方值的大小，卡方值越大咕晋，越不符合雹拄；卡方值越小，偏差越小掌呜，越趨于符合滓玖，若兩個(gè)值完全相等時(shí)，卡方值就為0站辉，表明理論值完全符合呢撞。
適用條件：
1.所有的理論數(shù)T≥5并且總樣本量n≥40,用Pearson卡方進(jìn)行檢驗(yàn).
2.如果理論數(shù)T＜5但T≥1,并且n≥40,用連續(xù)性校正的卡方進(jìn)行檢驗(yàn).
3.如果有理論數(shù)T＜1或n＜40,則用Fisher’s檢驗(yàn).
示例：
判斷5種品牌啤酒的愛好者有無(wú)顯著差異：
x<-c(210,312,170,85,223)
chisq.test(x)
x<-matrix(c(46,18,6,8),ncol=2,nrow=2)
chisq.test(x)
chisq.test(x)$expected ###查看理論值
fisher.test(x) ##進(jìn)行fisher檢驗(yàn)
秩和檢驗(yàn)
秩和檢驗(yàn)是對(duì)原假設(shè)的非參數(shù)檢驗(yàn)，在不需要假設(shè)兩個(gè)樣本空間都為正態(tài)分布的情況下饰剥，測(cè)試它們的分布是否完全相同殊霞。
library(stats)
data("mtcars")
wilcox.test(mpg~am,data = mtcars)
幾種常用的R包
目前常用差異分析的R包有edgeR、limma汰蓉、DESeq2
三種包的區(qū)別：
1.limma包做差異分析要求數(shù)據(jù)滿足正態(tài)分布或近似正態(tài)分布绷蹲，如基因芯片术幔、TPM格式的高通量測(cè)序數(shù)據(jù)躯护。
2.通常認(rèn)為Count數(shù)據(jù)不符合正態(tài)分布而服從泊松分布蒲凶。對(duì)于count數(shù)據(jù)來(lái)說浸赫，用limma包做差異分析，誤差較大
3.DESeq2裙士、和 EdgeR都是基于count蒿辙，然后兩個(gè)都是NB（negative binomial）但是在估計(jì)dispersion parameter的方法上面不一樣院峡。
4.limma测秸，edgeR疤估，DESeq2三大包基本是做轉(zhuǎn)錄組差異分析的金標(biāo)準(zhǔn)灾常，大多數(shù)轉(zhuǎn)錄組的文章都是用這三個(gè)R包進(jìn)行差異分析。
5.edgeR差異分析速度快铃拇，得到的基因數(shù)目比較多钞瀑，假陽(yáng)性高（實(shí)際不差異，結(jié)果差異）慷荔。DESeq2差異分析速度慢雕什，得到的基因數(shù)目比較少，假陰性高（實(shí)際差異显晶，結(jié)果不差異）贷岸。
6.需要注意的是制作分組信息的因子向量是，因子水平的前后順序吧碾，在R的很多模型中凰盔，默認(rèn)將因子向量的第一個(gè)水平看作對(duì)照組
如果數(shù)據(jù)量大并且要求比較conservative的話可以所有方法都用下，然后取并集倦春。
數(shù)據(jù)預(yù)處理
library("edgeR")
expr = read.csv("mRNA_exprSet.csv",sep = ',',header=T)
head(expr)
expr = avereps(expr[,-1],ID = expr$X) # 對(duì)重復(fù)基因名取平均表達(dá)量，然后將基因名作為行名
expr = expr[rowMeans(expr)>1,] # 去除低表達(dá)的基因
library(stringr)
tumor <- colnames(expr)[as.integer(substr(colnames(expr),14,15)) < 10]
normal <- colnames(expr)[as.integer(substr(colnames(expr),14,15)) >= 10]
tumor_sample <- expr[,tumor]
normal_sample <- expr[,normal]
exprSet_by_group <- cbind(tumor_sample,normal_sample)
group_list <- c(rep('tumor',ncol(tumor_sample)),rep('normal',ncol(normal_sample)))
edgeR
對(duì)于edgeR的分析流程而言落剪，我們需要輸入的數(shù)據(jù)包括：
- 表達(dá)矩陣（
counts） - 分組信息（
group） - 擬合信息（
design）：指明如何根據(jù)樣本的分組進(jìn)行建模
edgeR默認(rèn)使用 trimmed mean of M-values (TMM) 計(jì)算文庫(kù)的scale factor進(jìn)行normalization睁本，以最大程度地縮小樣本間基因表達(dá)量的log-fold change。這是因?yàn)門MM 法認(rèn)為樣本間大部分的基因都沒有發(fā)生差異表達(dá)忠怖，而那些真正差異表達(dá)的基因并不會(huì)受到normalization的嚴(yán)重影響呢堰。如此一來(lái)，便將那些由于測(cè)序引起的差異表達(dá)基因的表達(dá)量給校正了凡泣，消除了一部分的假陽(yáng)性枉疼。
data = exprSet_by_group
group_list = factor(group_list)
design <- model.matrix(~0+group_list)
rownames(design) = colnames(data)
colnames(design) <- levels(group_list)
DGElist <- DGEList( counts = data, group = group_list)
## Counts per Million or Reads per Kilobase per Million
keep_gene <- rowSums( cpm(DGElist) > 1 ) >= 2 ## 自定義
table(keep_gene)
DGElist <- DGElist[ keep_gene, , keep.lib.sizes = FALSE ]
DGElist <- calcNormFactors( DGElist )
DGElist <- estimateGLMCommonDisp(DGElist, design)##為所有基因都計(jì)算同樣的dispersion
DGElist <- estimateGLMTrendedDisp(DGElist, design)##根據(jù)不同基因的均值--方差之間的關(guān)系來(lái)計(jì)算dispersion，并擬合一個(gè)trended model
DGElist <- estimateGLMTagwiseDisp(DGElist, design)##計(jì)算每個(gè)基因的dispersion鞋拟，并利用經(jīng)驗(yàn)性貝葉斯方法（empirical bayes）壓縮至Trend Didspersion骂维。個(gè)人認(rèn)為這一項(xiàng)相當(dāng)于GLM中每個(gè)基因的beta值
#####################負(fù)二項(xiàng)式廣義對(duì)數(shù)線性模型
fit <- glmFit(DGElist, design)
results <- glmLRT(fit, contrast = c(-1, 1))
####################類似然負(fù)二項(xiàng)式廣義對(duì)數(shù)線性模型
fit <- glmQLFit(dge, design, robust = TRUE) #擬合模型
lrt <- glmQLFTest(fit)
nrDEG_edgeR <- topTags(results, n = nrow(DGElist))
nrDEG_edgeR <- as.data.frame(nrDEG_edgeR)
head(nrDEG_edgeR)
padj = 0.01 # 自定義
foldChange= 2 # 自定義
nrDEG_edgeR_signif = nrDEG_edgeR[(nrDEG_edgeR$FDR < padj &
(nrDEG_edgeR$logFC>foldChange | nrDEG_edgeR$logFC<(-foldChange))),]
nrDEG_edgeR_signif = nrDEG_edgeR_signif[order(nrDEG_edgeR_signif$logFC),]
如果沒有生物學(xué)重復(fù)
library(edgeR)
##跟DESeq2一樣，導(dǎo)入數(shù)據(jù)贺纲，預(yù)處理（用了cpm函數(shù)）
exprSet<- read.table(file = "test.txt", sep = "\t", header = TRUE, row.names = 1, stringsAsFactors = FALSE)
group_list <- factor(c(rep("Contral",1),rep("Treat",1)))
##設(shè)置分組信息航闺，并做TMM標(biāo)準(zhǔn)化
exprSet <- DGEList(counts = exprSet, group = group_list)
bcv = 0.4 #設(shè)置BCV值
et <- exactTest(exprSet, dispersion=bcv^2)
write.csv(topTags(et, n = nrow(exprSet$counts)), 'result.csv', quote = FALSE) #輸出主要結(jié)果
limma,DEseq2的代碼實(shí)現(xiàn)可以在https://blog.csdn.net/weixin_43700050/article/details/98085127找到。
參考：
1.http://www.reibang.com/p/bb0bd72bc428.
2.https://zhuanlan.zhihu.com/p/57756620.
3.https://blog.csdn.net/weixin_43700050/article/details/98085127.
