DESeq双因素差异分析

作者: 小潤澤 | 来源:发表于2020-04-24 13:16 被阅读0次

DESeq双因素差异分析
你最关心的差异基因是怎么挑出来的？！
转录组数据差异分析
DESeq差异分析
转录组差异分析流程三大R包比较
RNASEQ分析入门笔记8-使用DESeq2进行表达差异分析
R|FPKM、RPKM差异分析
RNA-seq 详细教程： `DESeq2` 差异表达分析（7）
用DESeq2包来对RNA-seq数据进行差异分析
171025 差异分析和KEGG GO

单因素差异分析

对于DESeq来说，常用的单因素的差异分析已经比较多了，

一般来说，对于基因的表达矩阵，放在前面的为对照组，后面的为处理组，所以结果为后面的比较前面的
并且对于control组或者treat组这样的组来说，组内各样本的表达谱姑且认为是相同的

双因素差异分析

那么，对于双因素差异分析

所谓双因素，是指一个为主因素，另一个为次要因素，我们除了考虑treat与control的差异以外，分别在这两个组内，不同genotype也会造成彼此之间的差异

双因素差异分析结构

这就是各个因素之间的结构和内在联系，一般我们把设计矩阵的前面那几列作为对照

双因素实例

若不带交互项

dds <- makeExampleDESeqDataSet(n=100,m=18)
dds$genotype <- factor(rep(rep(c("I","II","III"),each=3),2))
#我们把主要因素condition写在后面，genotype为次要因素写在前面，“：”表示交互项
design(dds) <- ~ genotype + condition
dds <- DESeq(dds)
resultsNames(dds)

我们可以发现，若没有交互项，各组的比较仅停留在单因素两两进行比较。
设计矩阵的因子水平如下：

设计矩阵的因子水平

第一种情况：

dat = counts(dds, normalized=TRUE)
dds_result = data.frame(results(dds, contrast=c("condition","B","A")))

这样的差异表达方式只考虑condition B和A的差异，而不区分是哪一种genotype
我们挑选前三个基因手动计算它们的log2(FC)值：

# dat 的 10-18列包括了所有 genotype I-III 的 conditionB 的 sample
# dat 的 1-9列包括了所有 genotype I-III 的 conditionA 的 sample

## Gene 1
log2(mean(dat[1,10:18])/mean(dat[1,1:9]))
[1] -0.5054413

## Gene 2
log2(mean(dat[2,10:18])/mean(dat[2,1:9]))
[1] -0.09992514

## Gene 3
log2(mean(dat[3,10:18])/mean(dat[3,1:9]))
[1] -0.2876421

而软件计算的log2(FC)值如下（忽略矫正的情况）：

Gene 1: log2(FC) 为 -0.495586066

Gene 2: log2(FC) 为 -0.103780171

Gene 3: log2(FC) 为 -0.322288959

第二种情况：

dat = counts(dds, normalized=TRUE)
dds_result = data.frame(results(dds, contrast=c("genotype","II","I")))

这样的差异表达方式只考虑genotype II 和 I 之间的差异，而不区分是哪一种condition
我们挑选前三个基因手动计算它们的log2(FC)值：

# dat 的 4-6列和13-15列包括了所有 condition 的 genotype II 的 sample
# dat 的 1-3列和10-12列包括了所有 condition 的 genotype I 的 sample

## Gene 1
log2(mean(dat[1,c(4:6,13:15)])/mean(dat[1,c(1:3,10:12)]))
[1] -0.3342552

# Gene 2
log2(mean(dat[2,c(4:6,13:15)])/mean(dat[2,c(1:3,10:12)]))
[1] 0.4529269

# Gene 3
log2(mean(dat[3,c(4:6,13:15)])/mean(dat[3,c(1:3,10:12)]))
[1] -0.1632169

而软件计算的log2(FC)值如下（忽略矫正的情况）：

Gene 1: log2(FC) 为 -0.3179995795

Gene 2: log2(FC) 为 0.4551878042

Gene 3: log2(FC) 为 -0.2035258259

第三种情况：

dat = counts(dds, normalized=TRUE)
dds_result = data.frame(results(dds, contrast=c("genotype","III","I")))

这样的差异表达方式只考虑genotype III 和 I 之间的差异，而不区分是哪一种condition
我们挑选前三个基因手动计算它们的log2(FC)值：

# dat 的 7-9列和16-18列包括了所有 condition 的 genotype III 的 sample
# dat 的 1-3列和10-12列包括了所有 condition 的 genotype I 的 sample

## Gene 1
log2(mean(dat[1,c(7:9,16:18)])/mean(dat[1,c(1:3,10:12)]))
[1] -0.5673898

# Gene 2
log2(mean(dat[2,c(7:9,16:18)])/mean(dat[2,c(1:3,10:12)]))
[1] 0.7163037

# Gene 3
log2(mean(dat[3,c(7:9,16:18)])/mean(dat[3,c(1:3,10:12)]))
[1] -0.2464887

而软件计算的log2(FC)值如下（忽略矫正的情况）：

Gene 1: log2(FC) 为 -0.5515813451

Gene 2: log2(FC) 为 0.7185941762

Gene 3: log2(FC) 为 -0.3108995692

第四种情况：

dat = counts(dds, normalized=TRUE)
dds_result = data.frame(results(dds, contrast=c("genotype","III","II")))

这样的差异表达方式只考虑genotype III 和 II 之间的差异，而不区分是哪一种condition
我们挑选前三个基因手动计算它们的log2(FC)值：

# dat 的 7-9列和16-18列包括了所有 condition 的 genotype III 的 sample
# dat 的 4-6列和13-15列包括了所有 condition 的 genotype II 的 sample

## Gene 1
log2(mean(dat[1,c(7:9,16:18)])/mean(dat[1,c(4:6,13:15)]))
[1] -0.2331346

# Gene 2
log2(mean(dat[2,c(7:9,16:18)])/mean(dat[2,c(4:6,13:15)]))
[1] 0.2633768

# Gene 3
log2(mean(dat[3,c(7:9,16:18)])/mean(dat[3,c(4:6,13:15)]))
[1] -0.08327181

而软件计算的log2(FC)值如下（忽略矫正的情况）：

Gene 1: log2(FC) 为 -0.23358177

Gene 2: log2(FC) 为 0.26340637

Gene 3: log2(FC) 为 -0.10737374

若带有交互项

我们用帮助文档中的例子作为说明

dds <- makeExampleDESeqDataSet(n=100,m=18)
dds$genotype <- factor(rep(rep(c("I","II","III"),each=3),2))
#我们把主要因素condition写在后面，genotype为次要因素写在前面，“：”表示交互项
design(dds) <- ~ genotype + condition + genotype:condition
dds <- DESeq(dds)
resultsNames(dds)

如果带有交互项，那么我们就可以比较双因素共同作用下两组的差异了；也可以限定一个因素，比较某单因素作用下两组的差异

注意事项

DESeq2的建模原理是利用广义线性模型去拟合高维数据，然后再做差异比较
值得注意的是，不带有交互项（design(dds) <- ~ genotype + condition）对高维数据的拟合的方程是线性的（两个维度，其中一个维度是genotype，另一个维度是condition）；带有交互项（design(dds) <- ~ genotype + condition + genotype:condition）对高维数据的拟合的方程是非线性的（两个维度，其中一个维度是genotype，另一个维度是condition ，再带有一个非线性交互项 genotype:condition）；因此两个模型计算出来的基因表达估计值是不一样的，因此计算出来的log2(FC)是不一样的，但是不会影响上下调的方向

上面是没有交互项的；下面是有交互项的

那么DESeq的原理是哪一个sample放前面，哪一个就为对照组，在该例子里面，对于condition，condition_A为对照，对于genotype，genotype I为对照组
当然我们也可以自行设定：

# to compare A vs B, make B the reference level,
# and select the last coefficient
condition <- relevel(condition, "B")

genotype_II_vs_I 表示的是在condition_A的条件下，genotype II 与genotype I 相比较下的差异

genotype_III_vs_I 表示的是在condition_A的条件下，genotype III 与genotype I 相比较下的差异

condition_B_vs_A 表示对于genotype I 来说，在condition_B与condition_A相比较下的差异

genotypeII.conditionB 表示双因素交互项，目的是检测在不同genotype中，condition的效果是否会不同，即先计算condition_B的genotype II 与condition_A的genotype II 相比较下的差异和condition_B的genotype I 与condition_A的genotype I 的差异，两者再进行差异比较（前者比后者）

genotypeIII.conditionB 表示双因素交互项，目的是检测在不同genotype中，condition的效果是否会不同，即先计算condition_B的genotype III 与condition_A的genotype III 相比较下的差异和condition_B的genotype I 与condition_A的genotype I 的差异，两者再进行差异比较（前者比后者）

接下来我们介绍下参数 contrast：

这个参数目的是指定要从对象中提取哪些比较以构建结果表
简而言之，这个参数控制着提取谁与谁比较的结果

第一种情况：
那么：

dat = counts(dds, normalized=TRUE)
res = results(dds, list( c("genotype_II_vs_I")))

这个表示对于condition_A 来说，在genotype II 与genotype I 相比较下的差异
FC＝(genotype II in condition_A) / (genotype I in condition_A)
我们用代码计算下：

# dat 的 4-6列代表 genotype II in condition_A 的 sample
# dat 的 1-3列代表 genotype I in condition_A 的 sample

## Gene 1
log2(mean(dat[1,4:6])/mean(dat[1,1:3]))
[1] 0.3768983

## Gene 2
log2(mean(dat[2,4:6])/mean(dat[2,1:3]))
[1] 0.3134871

## Gene 3
log2(mean(dat[3,4:6])/mean(dat[3,1:3]))
[1] -0.5114681

而软件计算的log2(FC)值如下（忽略矫正的情况）：

Gene 1: log2(FC) 为 0.367333

Gene 2: log2(FC) 为 0.313593

Gene 3: log2(FC) 为 -0.506202

第二种情况：
那么：

dat = counts(dds, normalized=TRUE)
res = results(dds, list( c("genotype_III_vs_I")))

这个表示对于condition_A 来说，在genotype III 与genotype I 相比较下的差异
FC＝(genotype III in condition_A) / (genotype I in condition_A)
我们用代码计算下：

# dat 的 7-9列代表 genotype III in condition_A 的 sample
# dat 的 1-3列代表 genotype I in condition_A 的 sample

## Gene 1
log2(mean(dat[1,7:9])/mean(dat[1,1:3]))
[1] -0.3639434

## Gene 2
log2(mean(dat[2,7:9])/mean(dat[2,1:3]))
[1] -0.03030897

## Gene 3
log2(mean(dat[3,7:9])/mean(dat[3,1:3]))
[1] -0.6220792

而软件计算的log2(FC)值如下（忽略矫正的情况）：

Gene 1: log2(FC) 为 -0.364814

Gene 2: log2(FC) 为 -0.029441

Gene 3: log2(FC) 为 -0.610125

第三种情况：
那么：

dat = counts(dds, normalized=TRUE)
# the condition effect for genotype I (the main effect)
results(dds, contrast=c("condition","B","A"))

这个表示对于genotype I 来说，在condition_B与condition_A相比较下的差异
FC＝(genotype I in condition_B) / (genotype I in condition_A)
我们用代码计算下：

# dat 的 10-12列代表 genotype I in condition_B 的 sample
# dat 的 1-3列代表 genotype I in condition_A 的 sample

## Gene 1
log2(mean(dat[1,10:12])/mean(dat[1,1:3]))
[1] -0.6302513

## Gene 2
log2(mean(dat[2,10:12])/mean(dat[2,1:3]))
[1] -0.1265862

## Gene 3
log2(mean(dat[3,10:12])/mean(dat[3,1:3]))
[1] 0.3226792

而软件计算的log2(FC)值如下（忽略矫正的情况）：

Gene 1: log2(FC) 为 -0.633304150

Gene 2: log2(FC) 为 -0.129958993

Gene 3: log2(FC) 为 0.321771376

第四种情况：

dat = counts(dds, normalized=TRUE)
# the condition effect for genotype III.
# this is the main effect *plus* the interaction term
# (the extra condition effect in genotype III compared to genotype I).
results(dds, list( c("condition_B_vs_A","genotypeIII.conditionB")))

表示主要考察在genotype III 这个因素下，condition_B 与condition_A 相比带来的差异。
也就是说，condition_B的genotype III 与 condition_A的genotype III 相比较所带来的差异
这里的extra condition effect应该理解为在genotype III 的条件下，某两个condition相比较下的差异
FC ＝ (genotype III in condition_B) / (genotype III in condition_A)

# dat 的 16-18列代表 genotype III in condition_B 的 sample
# dat 的 7-9列代表 genotype III in condition_A 的 sample

## Gene 1
log2(mean(dat[1,16:18])/mean(dat[1,7:9]))
[1] -0.4690318

## Gene 2
log2(mean(dat[2,16:18])/mean(dat[2,7:9]))
[1] -0.0880866

## Gene 3
log2(mean(dat[3,16:18])/mean(dat[3,7:9]))
[1] -0.8147309

而软件计算的log2(FC)值如下（忽略矫正的情况）：

Gene 1: log2(FC) 为 -0.46013841

Gene 2: log2(FC) 为 -0.08729557

Gene 3: log2(FC) 为 -0.81238535

第五种情况：

dat = counts(dds, normalized=TRUE)
# the condition effect for genotype II.
# this is the main effect *plus* the interaction term
# (the extra condition effect in genotype II compared to genotype I).
results(dds, list( c("condition_B_vs_A","genotypeII.conditionB")))

表示主要考察在genotype II 这个因素下，condition_B 与condition_A 相比带来的差异。
也就是说，condition_B的genotype II 与 condition_A的genotype II 相比较所带来的差异
这里的extra condition effect应该理解为在genotype II 的条件下，某两个condition相比较下的差异
FC ＝ (genotype II in condition_B) / (genotype II in condition_A)

# dat 的 13-15列代表 genotype II in condition_B 的 sample
# dat 的 4-6列代表 genotype II in condition_A 的 sample

## Gene 1
log2(mean(dat[1,13:15])/mean(dat[1,4:6]))
[1] 1.375654

## Gene 2
log2(mean(dat[2,13:15])/mean(dat[2,4:6]))
[1] -0.3227521

## Gene 3
log2(mean(dat[3,13:15])/mean(dat[3,4:6]))
[1] 0.08036955

而软件计算的log2(FC)值如下（忽略矫正的情况）：

Gene 1: log2(FC) 为 1.3656608

Gene 2: log2(FC) 为 -0.3305137

Gene 3: log2(FC) 为 0.0780536

第六种情况：

dat = counts(dds, normalized=TRUE)
# the interaction term for condition effect in genotype III vs genotype I.
# this tests if the condition effect is different in III compared to I
results(dds, name="genotypeIII.conditionB")

这一项是交互项，即目的是检测在不同genotype中，condition的效果是否会不同，即先计算condition_B的genotype III 与condition_A的genotype III 相比较下的差异和condition_B的genotype I 与condition_A的genotype I 的差异，两者再进行差异比较（前者比后者），英文的解释为：This tests if the condition effect is different in III compared to I:

FC ＝ (genotype III in condition_B / genotype III in condition_A) / (genotype I in condition_B / genotype I in condition_A)

这个 FC 主要考察的是 condition 的改变对 genotype III 产生的变化倍数是否和 condition 的改变对 genotype I 产生的变化倍数相同；如果相同 FC的值为1；如果前者变化倍数比后者小，那么 FC 介于0到1，如果前者变化倍数比后者大，那么 FC 大于1
代码如下：

# dat 的 16-18列代表 genotype III in condition_B 的 sample
# dat 的 7-9列代表 genotype III in condition_A 的 sample
# dat 的 10-12列代表 genotype I in condition_B 的 sample
# dat 的 1-3列代表 genotype I in condition_A 的 sample

## Gene 1
log2((mean(dat[1,16:18])/mean(dat[1,7:9]))/(mean(dat[1,10:12])/mean(dat[1,1:3])))
[1] 0.1612196

## Gene 2
log2((mean(dat[2,16:18])/mean(dat[2,7:9]))/(mean(dat[2,10:12])/mean(dat[2,1:3])))
[1] 0.03849963

## Gene 3
log2((mean(dat[3,16:18])/mean(dat[3,7:9]))/(mean(dat[3,10:12])/mean(dat[3,1:3])))
[1] -1.13741

而软件计算的log2(FC)值如下（忽略矫正的情况）：

Gene 1: log2(FC) 为 0.17316574

Gene 2: log2(FC) 为 0.04266342

Gene 3: log2(FC) 为 -1.13415673

这里在插播一下name这个参数：

通常只提取一项的话，用 name 比较好

第七种情况：

dat = counts(dds, normalized=TRUE)
# the interaction term for condition effect in genotype III vs genotype II.
# this tests if the condition effect is different in III compared to II
results(dds, contrast=list("genotypeIII.conditionB", "genotypeII.conditionB"))

这个是两个交互项的比较，其英文解释为：This tests if the condition effect is different in III compared to II
目的是检测在不同genotype中，condition的效果是否会不同，即先计算condition_B的genotype III 与condition_A的genotype III 相比较下的差异和condition_B的genotype II 与condition_A的genotype II 的差异，两者再进行差异比较（前者比后者）
FC ＝ (genotype III in condition_B / genotype III in condition_A) / (genotype II in condition_B / genotype II in condition_A)

这个 FC 主要考察的是 condition 的改变对 genotype III 产生的变化倍数是否和 condition 的改变对 genotype II 产生的变化倍数相同；如果相同 FC的值为1；如果前者变化倍数比后者小，那么 FC 介于0到1，如果前者变化倍数比后者大，那么 FC 大于1

# dat 的 16-18列代表 genotype III in condition_B 的 sample
# dat 的 7-9列代表 genotype III in condition_A 的 sample
# dat 的 13-15列代表 genotype II in condition_B 的 sample
# dat 的 4-6列代表 genotype II in condition_A 的 sample

## Gene 1
log2((mean(dat[1,16:18])/mean(dat[1,7:9]))/(mean(dat[1,13:15])/mean(dat[1,4:6])))
[1] -0.08777481

## Gene 2
log2((mean(dat[2,16:18])/mean(dat[2,7:9]))/(mean(dat[2,13:15])/mean(dat[2,4:6])))
[1] 0.006579703

## Gene 3
log2((mean(dat[3,16:18])/mean(dat[3,7:9]))/(mean(dat[3,13:15])/mean(dat[3,4:6])))
[1] -0.3210265

而软件计算的log2(FC)值如下（忽略矫正的情况）：

Gene 1: log2(FC) 为 -0.072876535

Gene 2: log2(FC) 为 0.007112289

Gene 3: log2(FC) 为 -0.321480534

第八种情况：

dat = counts(dds, normalized=TRUE)
res = results(dds, list( c("genotype_III_vs_I","genotypeIII.condition_B")))

这个比较的是在condition_B的条件下，genotype III 和genotype I 相比较下的差异
即condition_B的genotype III 与condition_B的genotype I 相比较下的差异
FC = (genotype III in condition_B) / (genotype I in condition_B)
手动计算代码：

# dat 的 16-18列代表 genotype III in condition_B 的 sample
# dat 的 10-12列代表 genotype I in condition_B 的 sample

## Gene 1
log2(mean(dat[1,16:18])/mean(dat[1,10:12]))
[1] -0.4716114

## Gene 2
log2(mean(dat[2,16:18])/mean(dat[2,10:12]))
[1] 0.7362693

## Gene 3
log2(mean(dat[3,16:18])/mean(dat[3,10:12]))
[1] -0.8629554

而软件计算的log2(FC)值如下（忽略矫正的情况）：

Gene 1: log2(FC) 为 -0.46202782

Gene 2: log2(FC) 为 0.73984754

Gene 3: log2(FC) 为 -0.86640184

第九种情况：

dat = counts(dds, normalized=TRUE)
res = results(dds, list( c("genotype_II_vs_I","genotypeII.conditionB")))

这个比较的是在condition_B的条件下，genotype II 和genotype I 相比较下的差异
即condition_B的genotype II 与condition_B的genotype I 相比较下的差异
FC = (genotype II in condition_B) / (genotype I in condition_B)
手动计算代码：

# dat 的 13-15列代表 genotype II in condition_B 的 sample
# dat 的 10-12列代表 genotype I in condition_B 的 sample

## Gene 1
log2(mean(dat[1,13:15])/mean(dat[1,10:12]))
[1] -0.9762591

## Gene 2
log2(mean(dat[2,13:15])/mean(dat[2,10:12]))
[1] -0.4114295

## Gene 3
log2(mean(dat[3,13:15])/mean(dat[3,10:12]))
[1] 1.126437

而软件计算的log2(FC)值如下（忽略矫正的情况）：

Gene 1: log2(FC) 为 -0.972029

Gene 2: log2(FC) 为 -0.412843

Gene 3: log2(FC) 为 1.132502

第十种情况：

dat = counts(dds, normalized=TRUE)
# 注意要分为两个 c() 来写
res = results(dds, list( c("genotype_III_vs_I","genotypeIII.conditionB"),c("genotype_II_vs_I","genotypeII.conditionB")))

这个比较的是在condition_B的条件下，genotype III 和genotype II 相比较下的差异
即condition_B的genotype III 与condition_B的genotype II 相比较下的差异
FC = (genotype III in condition_B) / (genotype II in condition_B)
手动计算代码：

# dat 的 16-18列代表 genotype III in condition_B 的 sample
# dat 的 13-15列代表 genotype II in condition_B 的 sample

## Gene 1
log2(mean(dat[1,16:18])/mean(dat[1,13:15]))
[1] 0.2984419

## Gene 2
log2(mean(dat[2,16:18])/mean(dat[2,13:15]))
[1] -0.1667908

## Gene 3
log2(mean(dat[3,16:18])/mean(dat[3,13:15]))
[1] -0.02535796

而软件计算的log2(FC)值如下（忽略矫正的情况）：

Gene 1: log2(FC) 为 0.3033151

Gene 2: log2(FC) 为 -0.1632627

Gene 3: log2(FC) 为 -0.0172830

如果我们想要计算跨condition和跨genotype的差异

第十一种情况：

dds <- makeExampleDESeqDataSet(n=100,m=18)
dds$genotype <- factor(rep(rep(c("I","II","III"),each=3),2))
dds$genotype <- relevel(dds$genotype, "III")
#我们把主要因素condition写在后面，genotype为次要因素写在前面，“：”表示交互项
design(dds) <- ~ genotype + condition + genotype:condition
dds <- DESeq(dds)
resultsNames(dds)


dat = counts(dds, normalized=TRUE)
res = results(dds, list(c("condition_B_vs_A"),c("genotype_I_vs_III")))

这里计算的是condition_B的genotype III 与condition_A的genotype I 相比较下的差异
FC = (genotype III in condition_B) / (genotype I in condition_A)

手动计算代码：

# dat 的 16-18列代表 genotype III in condition_B 的 sample
# dat 的 1-3列代表 genotype I in condition_A 的 sample

## Gene 1
log2(mean(dat[1,16:18])/mean(dat[1,1:3]))
[1] -0.2684356

## Gene 2
log2(mean(dat[2,16:18])/mean(dat[2,1:3]))
[1] 0.6341815

## Gene 3
log2(mean(dat[3,16:18])/mean(dat[3,1:3]))
[1] 0.111087

而软件计算的log2(FC)值如下（忽略矫正的情况）：

Gene 1: log2(FC) 为 -0.267888

Gene 2: log2(FC) 为 0.638842

Gene 3: log2(FC) 为 0.108176

第十二种情况：

dds <- makeExampleDESeqDataSet(n=100,m=18)
dds$genotype <- factor(rep(rep(c("I","II","III"),each=3),2))
dds$genotype <- relevel(dds$genotype, "III")
#我们把主要因素condition写在后面，genotype为次要因素写在前面，“：”表示交互项
design(dds) <- ~ genotype + condition + genotype:condition
dds <- DESeq(dds)
resultsNames(dds)


dat = counts(dds, normalized=TRUE)
res = results(dds, list(c("condition_B_vs_A"),c("genotype_II_vs_III")))

这里计算的是condition_B的genotype III 与condition_A的genotype II 相比较下的差异
FC = (genotype III in condition_B) / (genotype II in condition_A)

手动计算代码：

# dat 的 16-18列代表 genotype III in condition_B 的 sample
# dat 的 4-6列代表 genotype II in condition_A 的 sample

## Gene 1
log2(mean(dat[1,16:18])/mean(dat[1,4:6]))
[1] 0.07651671

## Gene 2
log2(mean(dat[2,16:18])/mean(dat[2,4:6]))
[1] 0.966533

## Gene 3
log2(mean(dat[3,16:18])/mean(dat[3,4:6]))
[1] -0.1500791

而软件计算的log2(FC)值如下（忽略矫正的情况）：

Gene 1: log2(FC) 为 0.0778375

Gene 2: log2(FC) 为 0.9659273

Gene 3: log2(FC) 为 -0.1501982

第十三种情况：

dds <- makeExampleDESeqDataSet(n=100,m=18)
dds$genotype <- factor(rep(rep(c("I","II","III"),each=3),2))
dds$genotype <- relevel(dds$genotype, "II")
#我们把主要因素condition写在后面，genotype为次要因素写在前面，“：”表示交互项
design(dds) <- ~ genotype + condition + genotype:condition
dds <- DESeq(dds)
resultsNames(dds)


dat = counts(dds, normalized=TRUE)
res = results(dds, list(c("condition_B_vs_A"),c("genotype_I_vs_II")))

这里计算的是condition_B的genotype II 与condition_A的genotype I 相比较下的差异
FC = (genotype II in condition_B) / (genotype I in condition_A)

手动计算代码：

# dat 的 13-15列代表 genotype II in condition_B 的 sample
# dat 的 1-3列代表 genotype I in condition_A 的 sample

## Gene 1
log2(mean(dat[1,13:15])/mean(dat[1,1:3]))
[1] 1.852217

## Gene 2
log2(mean(dat[2,13:15])/mean(dat[2,1:3]))
[1] 0.9619238

## Gene 3
log2(mean(dat[3,13:15])/mean(dat[3,1:3]))
[1] -0.5547154

而软件计算的log2(FC)值如下（忽略矫正的情况）：

Gene 1: log2(FC) 为 1.812546

Gene 2: log2(FC) 为 0.959211

Gene 3: log2(FC) 为 -0.553264

关于contract的注意事项

contract这个参数后面接contract = c("..." , "...")，contract = list(c("..." , "..."))和contract = list(c("..."),c("..."))这几个所对比的是不一样的
第一个要求的向量长度为3，那么只适用于某一个因素不同的两组进行对比
例如：

对于上图的condition_B_vs_A，这是一个条件，所以有：

results(dds, contrast=c("condition","B","A"))

那么对于多因素的情况，比如下面的情况

results(dds, contrast=list( c("condition_B_vs_A","genotypeIII.conditionB") ))

主要考察在genotype III 这个因素下，condition_B 与condition_A 相比带来的差异。
也就是说，condition_B的genotype III 与 condition_A的genotype III 相比较所带来的差异

那么下面的这种情况：

results(dds, contrast=list( c("genotype_III_vs_I","genotypeIII.conditionB"),c("genotype_II_vs_I","genotypeII.conditionB")))

前一个c()作分子，后一个c()作分母；
表示的是在condition_B条件下，genotype III 与genotype I 相比较下的差异（分子）和condition_B条件下，genotype II 与genotype I 相比较下的差异（分母），两者再进行比较下的差异（再相除）

FC = [(genotype III in condition_B) / (genotype I in condition_B)] / (genotype II in condition_B) / (genotype I in condition_B)] = (genotype III in condition_B) / (genotype II in condition_B)，相当于约分

参考：《DESeq2说明文档》
https://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
Bioconductor issue：
https://support.bioconductor.org/p/121674/
主要解释：
https://rstudio-pubs-static.s3.amazonaws.com/329027_593046fb6d7a427da6b2c538caf601e1.html

DESeq双因素差异分析
单因素差异分析对于DESeq来说，常用的单因素的差异分析已经比较多了，一般来说，对于基因的表达矩阵，放在前面的...
你最关心的差异基因是怎么挑出来的？！
差异基因分析做基因表达分析时必然要做差异基因分析，做差异基因分析最常用的软件就是DESeq2，使用DESeq2对...
转录组数据差异分析
差异分析• limma• edgeR• DESeq2 title: "三大R包差异分析"output: html_...
DESeq差异分析
基因表达矩阵,readcounts
转录组差异分析流程三大R包比较
总目录：数据预处理 limma包进行差异分析 edgeR包进行差异分析 DESeq2包进行差异分析三大R包比较...
RNASEQ分析入门笔记8-使用DESeq2进行表达差异分析
基因的差异表达分析，通常使用R中的软件包，包括：DESeq2，edgeR，limma等，今天介绍DESeq2的分析...
R|FPKM、RPKM差异分析
芯片数据差异分析，常规用limma进行差异分析，而对于RNA-seq数据，常用edgeR、DEseq2和limma...
RNA-seq 详细教程： `DESeq2` 差异表达分析（7）
学习目标了解如何设计公式了解如何使用 DESeq2 执行差异表达分析 1. DE 分析差异表达分析工作流程的...
用DESeq2包来对RNA-seq数据进行差异分析
转自jmzeng 用DESeq2包来对RNA-seq数据进行差异分析 jmzeng 2016年4月11日差异分析...
171025 差异分析和KEGG GO
物种Verticillium dahliae Kleb 1.deseq2差异分析只有count数据。按照coun...