DESeq双因素差异分析

单因素差异分析

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



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

双因素差异分析

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



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

双因素差异分析的交互项：
如果不考虑交互项，那么意味着在 genotype I ，genotype II 和 genotype III 中 treat 比较于 control 的差异值都是相等的；而事实并非如此，在 genotype I ，genotype II 和 genotype III 中 treat 比较于 control 的差异值很可能是不一样的，因此需要引入交互项来评价在 genotype I ，genotype II 和 genotype III 中 treat 比较于 control 差异值的不同


交互项的意义是当控制第二个变量时，我们讨论 control 与treat 的差异，那么效应值为 a2，这意味着在第二个变量 genotype I ，genotype II 和 genotype III 中 control 与treat 的差异都是 a2。可是不是在第二个变量 genotype I ，genotype II 和 genotype III 类中的差异都是 a2？需要引入交互项来区分这种差异
因此引入交互项来区分在第二个变量 genotype I ，genotype II 和 genotype III 中 control 与treat 的差异，如上图 y2 和 y3，同样描述 control 与treat 的差异（效应值为 a2）但在第二个变量 genotype II 和 genotype III 中是不同的，在 genotype II 类中 control 与 treat 的差异为 a2+d5，在 genotype III 中 control 与 treat 的差异为a2+d6


交互项显著代表 control 与 treat 的差异在不同的 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)值如下（忽略矫正的情况）：

  1. Gene 1: log2(FC) 为 -0.495586066
  2. Gene 2: log2(FC) 为 -0.103780171
  3. 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)值如下（忽略矫正的情况）：

  1. Gene 1: log2(FC) 为 -0.3179995795
  2. Gene 2: log2(FC) 为 0.4551878042
  3. 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)值如下（忽略矫正的情况）：

  1. Gene 1: log2(FC) 为 -0.5515813451
  2. Gene 2: log2(FC) 为 0.7185941762
  3. 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)值如下（忽略矫正的情况）：

  1. Gene 1: log2(FC) 为 -0.23358177
  2. Gene 2: log2(FC) 为 0.26340637
  3. 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")
  1. genotype_II_vs_I 表示的是在condition_A的条件下，genotype II 与genotype I 相比较下的差异
  2. genotype_III_vs_I 表示的是在condition_A的条件下，genotype III 与genotype I 相比较下的差异
  3. condition_B_vs_A 表示对于genotype I 来说，在condition_B与condition_A相比较下的差异
  4. genotypeII.conditionB 表示双因素交互项，目的是检测在不同genotype中，condition的效果是否会不同，即先计算condition_B的genotype II 与condition_A的genotype II 相比较下的差异和condition_B的genotype I 与condition_A的genotype I 的差异，两者再进行差异比较（前者比后者）
  5. 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)值如下（忽略矫正的情况）：

  1. Gene 1: log2(FC) 为 0.367333
  2. Gene 2: log2(FC) 为 0.313593
  3. 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)值如下（忽略矫正的情况）：

  1. Gene 1: log2(FC) 为 -0.364814
  2. Gene 2: log2(FC) 为 -0.029441
  3. 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)值如下（忽略矫正的情况）：

  1. Gene 1: log2(FC) 为 -0.633304150
  2. Gene 2: log2(FC) 为 -0.129958993
  3. 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)值如下（忽略矫正的情况）：

  1. Gene 1: log2(FC) 为 -0.46013841
  2. Gene 2: log2(FC) 为 -0.08729557
  3. 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)值如下（忽略矫正的情况）：

  1. Gene 1: log2(FC) 为 1.3656608
  2. Gene 2: log2(FC) 为 -0.3305137
  3. 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)值如下（忽略矫正的情况）：

  1. Gene 1: log2(FC) 为 0.17316574
  2. Gene 2: log2(FC) 为 0.04266342
  3. 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)值如下（忽略矫正的情况）：

  1. Gene 1: log2(FC) 为 -0.072876535
  2. Gene 2: log2(FC) 为 0.007112289
  3. 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)值如下（忽略矫正的情况）：

  1. Gene 1: log2(FC) 为 -0.46202782
  2. Gene 2: log2(FC) 为 0.73984754
  3. 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)值如下（忽略矫正的情况）：

  1. Gene 1: log2(FC) 为 -0.972029
  2. Gene 2: log2(FC) 为 -0.412843
  3. 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)值如下（忽略矫正的情况）：

  1. Gene 1: log2(FC) 为 0.3033151
  2. Gene 2: log2(FC) 为 -0.1632627
  3. 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)值如下（忽略矫正的情况）：

  1. Gene 1: log2(FC) 为 -0.267888
  2. Gene 2: log2(FC) 为 0.638842
  3. 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)值如下（忽略矫正的情况）：

  1. Gene 1: log2(FC) 为 0.0778375
  2. Gene 2: log2(FC) 为 0.9659273
  3. 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)值如下（忽略矫正的情况）：

  1. Gene 1: log2(FC) 为 1.812546
  2. Gene 2: log2(FC) 为 0.959211
  3. 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

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 214,444评论 6赞 496
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 91,421评论 3赞 389
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 160,036评论 0赞 349
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 57,363评论 1赞 288
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 66,460评论 6赞 386
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 50,502评论 1赞 292
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 39,511评论 3赞 412
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 38,280评论 0赞 270
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 44,736评论 1赞 307
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,014评论 2赞 328
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,190评论 1赞 342
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 34,848评论 5赞 338
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 40,531评论 3赞 322
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,159评论 0赞 21
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,411评论 1赞 268
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 47,067评论 2赞 365
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 44,078评论 2赞 352

推荐阅读更多精彩内容

  • 2018-07-14
    选择题部分 1.(),只有在发生短路事故时或者在负荷电流较大时,变流器中才会有足够的二次电流作为继电保护跳闸之用。...
    skystarwuwei阅读 12,876评论 0赞 7
  • 日本食管学会之食管癌实践指南2017：（第1部分）
    引言指南目的这些指南的主要目的是为一般临床医生提供信息，指导他们对食管癌的可用诊断/治疗策略做出恰当的选择（针对食...
    医博云天阅读 4,310评论 0赞 2
  • 临床2
    ★681★触觉语颤增强主要见于 E.肺泡内炎症性浸润 ★682★女，12岁，白血病，近日左眼眶出现2cmX2cm无...
    捷龙阅读 4,516评论 0赞 2
  • 配置GLEW和GLFW
    首先在mac中要安装glew和glfw，这里我安装这两个工具使用的是homebrew包管理，这东西超级好用，安装命...
    yyyqy阅读 4,872评论 0赞 0
  • 2019-01-21
    20190121学习笔记： 什么才是受尊敬的企业？ 一家受到尊敬的公司， 理所应当为股东创造物质财富， 与此同时，...
    马唐阅读 321评论 0赞 0