====== 方差分析(ANOVA) ====== 方差分析是统计分析的一种常见的分析方法,对于老师实验设计比较复杂,比如平常农口实验设计多时间点的处理样本,我们常用的两两比较在处理这种时间序列的实验设计就感觉比较单一,不能满足老师的分析要求,因而此时可以使用的方差分析(ANOVA)完成分析。 ===== 2.1 one-way-ANOVA(单因素方差分析) ===== 还是以之前的一个项目为例 ^ 0h ^ 1h ^ 6h ^ 24h ^ 7d ^ 18d ^ 检验方法 ^ | eSR0h | eSR1h | eSR6h | eSR24h | eSR7d | eSR18d | one-way-anova | | eSR0h | eSR1h | eSR6h | eSR24h | | | short_term(one-way-anova) | | eSR0h | eSR7d | eSR18d | | | | long_term (one way anova) | 这个项目是老师采取的一系列时间点的处理的胡杨的根组织,老师的分析要求是探究在这连续的时间序列内的显著变化的差异基因,这是一个典型的单因素的方差分析。以处理和非处理的为单一因素,探究时间序列上差异基因变化。这种检验方法使用R里面的edgeR实现起来比较容易。 针对单一因素方差分析,其实在edgeR包内也不是一般意义下的单因素anova分析,在edgeR内该分析被称作An ANOVA-like test for any di erences。在edgeR中,我们设计design的时候可以查看design的情况,在差异基因分析是,有一个重要的参数来coef(coefficients ),选择共同表达的组别。coef表示的时你选择的共同表达分析的比较组状况。 libary(edgeR) libary(limma) setwd('/TJPROJ1/RNA/project/201511/Regcon/NHM150353_huyang_lnc/aftersale/ANOVA_multicompare/1st_dataset/') #设置工作目录 rc <- read.delim('first_dataset.xls', check.names=FALSE, stringsAsFactors=FALSE) #读入文件 group <- factor( c("control","control","control","treat1","treat1","treat1","treat2","treat2","treat2", "treat3","treat3","treat3","treat4","treat4","treat4","treat5","treat5","treat5"), levels <- c('control','treat1','treat2','treat3','treat4','treat5')) #设置分组及处理因素 y <- DGEList(counts=rc[, 2:19],genes=rc[, 1], group=group ) #edgeR读入数据数据 y <- calcNormFactors(y) #标准化 design <- model.matrix(~group) #设置实验设计 rownames(design) <- colnames(y) #导入列名称 y <- estimateGLMCommonDisp(y, design, verbose=TRUE) #计算离差 y <- estimateGLMTrendedDisp(y, design) y <- estimateGLMTagwiseDisp(y, design) fit <- glmFit(y, design) #线性拟合 lrt <- glmLRT(fit, coef=2:6) # 找出与实验设计变化相关的基因 ttg <- topTags( lrt, adjust.method="BH", sort.by="p.value", n=90587) #对差异分析结果进行FDR校验 write.table(ttg, file='./first_dataset_test_out.xls', sep='\t', quote=F, row.names=F, col.names=T) #导出运算结果 ===== 2.2 two-way-anova(双因素方差分析) ===== 以在做的一个项目为例,跟之前的单因素方差是一个项目,如下图。 ^ ^ salt ^control ^ | root | eSR24h | eCR24h | | stem | eSS24h | eCS24h | | leaf | eSL24h | eCL24h | library(edgeR) library(multcomp) fpkm <- read.delim("dataset.xls", sep="\t", header=T, quote="") treat <- factor(c(rep(c('control', 'control', 'control', 'salt', 'salt', 'salt'), 3))) tissue <-factor(c('root','root','root','root','root','root','stem','stem','stem','stem','stem','stem','leaf','leaf','leaf','leaf','leaf','leaf')) design < -model.matrix(~ treat+tissue) #实验设计(重要) dge<-DGEList(counts=fpkm[, 2:19], genes=fpkm[, 1], group=rep(c('R_control','R_salt','S_control','S_salt','L_control','L_salt'),each=3)) dge <- calcNormFactors(dge) dge <- estimateGLMCommonDisp(dge, design) dge <- estimateGLMTagwiseDisp(dge, design) fit <- glmFit(dge, design ) lrt.dge <- glmLRT(fit, coef=2:6 ) #设置共分析的组别(重要) xx<- topTags(lrt.dge, adjust.method="BH", sort.by="logFC", n=90587) write.table(xx, "two_way_anova_logfc.xls", sep="\t") 老师实验设计比较复杂,需要分析 在处理和非处理情况下,不同组织的基因表达差异分析,这个就包含了两对外在因素。处理和非处理,组织的不同(根,茎,叶)。要在这两个因素下进行差异分析的分析。 同时对于两个因素的处理也会有所考虑。针对双因素的方差分析,在design的设置上也有一些其他的注意事项。在本次试验中,双因素包括为两个因素:tissue 和 treat ,这两个因素之间没有交叉影响,因而在设置design的时候两因素之间使用的是“+”,但是如果两因素之间有交叉作用的话,实验设计应该使用tissue*treat。 针对与双因素方差分析的实验设计,同样查阅了《R语言实战》一书,书内对于常见研究设计有表达式,下表表达式针对于R包里普通的aov函数。 ^ 设计 ^表达式 ^ | 单因素ANOVA | y~A | | 含单个协变量的单因素ANCOVA | y~x+A | | 双因素ANOVA | y~A*B | | 含两个协变量的双因素ANCVOA | y~x1+x2+A*B | ======参考资料====== 生物统计学(第五版)科学出版社 李春喜、姜丽娜、邵云 等 R语言实战 (第三板) {{:知识:r语言实战_中文完整版_.pdf|}} edgeRUsersGuide {{:知识:edgerusersguide.pdf|}} T-test {{:知识:t检验.ppt|}}