DEseq2

new_rnaCounts 为count文件
new_metaMatrix2 为样本分组文件
group为new_metaMatrix2 中的一列

library(DESeq2)
# 构建DESeq2中的对象
dds <- DESeqDataSetFromMatrix(
  countData = new_rnaCounts,
  colData = new_metaMatrix2,
  design = ~ group)

# 指定哪一组作为control
dds$group <- relevel(dds$group, ref = "long")


dds <- DESeq(dds)
res <- results(dds)

在读取数据的过程中，有两点需要注意，第一个就是根据表达量对基因进行过滤，通常是过滤低表达量的基因，这一步是可选的，阈值可以自己定义；另外一个就是指定哪一组作为control组，在计算log2FD时 ，需要明确control组，默认会字符串顺序对分组的名字进行排序，排在前面的作为control组，这种默认行为选出的control可能与我们的实验设计不同，所以必须明确指定control组。

2. 归一化

计算每个样本的归一化系数，代码如下

dds <- estimateSizeFactors(dds)

将原始的表达量除以每个样本的归一化系数，就得到了归一化之后的表达量。

3. 估计基因的离散程度

DESeq2假定基因的表达量符合负二项分布，有两个关键参数，总体均值和离散程度α值， 如下图所示

image

这个α值衡量的是均值和方差之间的关系，表达式如下

image

通过下列代码估算每个基因的α值

dds <- estimateDispersions(dds)

4. 差异分析

代码如下

dds <- nbinomWaldTest(dds)
res <- results(dds)

为了简化调用，将第二部到第四部封装到了DESeq这个函数中，代码如下

dds <- DESeq(dds)
res <- results(dds)
write.table(
  res,
  "DESeq2.diff.tsv",
  sep="\t",
  quote=F,
  col.names = NA)

在DESeq2差异分析的过程中，已经考虑到了样本之间已有的差异，所以可以发现，最终结果里的log2FD值和我们拿归一化之后的表达量计算出来的不同, 示意如下

 head(results(dds)[, 1:2])
log2 fold change (MLE): condition B vs A
DataFrame with 6 rows and 2 columns
        baseMean log2FoldChange
       <numeric>      <numeric>
gene1   7.471250     -0.8961954
gene2  18.145279      0.4222174
gene3   2.329461     -2.3216915
gene4 165.634256     -0.1974001
gene5  38.300621      1.3573162
gene6   7.904819      1.8384322

提取归一化之后的表达量，自己计算baseMean和logFoldChange, 示例数据包含了12个样本，其中前6个样本为control, 后6个样本为case , 代码如下

> nor_count <- counts(dds, nor = T)
> res <- data.frame(
   baseMean = apply(nor_count, 1, mean),
   log2FoldChange = apply(nor_count, 1, function(t){ mean(t[7:12]) / mean(t[1:6])})
  )
> head(res)
        baseMean log2FoldChange
gene1   7.471250      0.5380191
gene2  18.145279      1.3404422
gene3   2.329461      0.1991979
gene4 165.634256      0.8719078
gene5  38.300621      2.5621035
gene6   7.904819      3.5365201

对比DESeq2和自己计算的结果，可以发现，baseMeans是一致的，而log2Foldchange 差异很大，甚至连数值的正负都发生了变化。

log2FD 反映的是不同分组间表达量的差异，这个差异由两部分构成，一种是样本间本身的差异，比如生物学重复样本间基因的表达量就有一定程度的差异，另外一部分就是我们真正感兴趣的，由于分组不同或者实验条件不同造成的差异。

用归一化之后的数值直接计算出的log2FD包含了以上两种差异，而我们真正感兴趣的只有分组不同造成的差异，DESeq2在差异分析的过程中已经考虑到了样本本身的差异，其最终提供的log2FD只包含了分组间的差异，所以会与手动计算的不同。