美文网首页TCGA
DESeq双因素差异分析

DESeq双因素差异分析

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

单因素差异分析

对于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)值如下(忽略矫正的情况):

  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

相关文章

网友评论

    本文标题:DESeq双因素差异分析

    本文链接:https://www.haomeiwen.com/subject/cuvzihtx.html