在做基因差异表达分析时，经常会用DESeq2这个包，但一直没有深究其分析的统计流程。因此，在这里记录一下DESeq2校正基因表达的方法 -- 比率中值法。

开篇明义，比率中值法考虑的因素是测序深度和文库的RNA组成。

以下我们使用airway数据具体展示比率中值法的校正过程

library(airway)
# Airway smooth muscle cells lines RNA-Seq experiment
data("airway")
# randomly sample 5k genes for testing
set.seed(123)
df1 <- as.data.frame(assay(airway)[sample(1:nrow(airway), 5000), 1:4])

比率中值法具体实现如下：

Step 1: 创建假参考样本 （每个基因的几何平均值）

首先，我们为每一个基因创建一个假参考样本(pseudo-reference sample)表达量，这代表理想的基因表达量。该值通过计算每个基因的几何平均值得到

几何平均值：n个数的乘积开n次根

gm <- apply(df1, 1, function(x) exp(mean(log(x))))
df1$pseudo <- gm
knitr::kable(head(df1))

Step 2: 计算每一个样本相对于参考样本的比值

接着，对于样本的每一个基因，我们计算一个 sample/reference 的比值。

# step 2: calculates ratio of each sample to the reference
ratios <- df1[,-5]/gm
colnames(ratios) <- paste0('R_',colnames(ratios))
df2 <- merge(df1, ratios, by=0)
knitr::kable(head(df2))

以第一行基因的ratio为例，写出具体计算过程：

由于大部分基因不会发生差异表达，所以样本内的ratios应是相当的。

Step 3: 计算每个样本的校正因子

样本内（column-wise）所有比率的中值即为校正因子(normalization factor), 也就是DESeq2 的size factor

# step 3: size factor
df2 <- na.omit(df2)

sf1 <- median(df2$R_SRR1039508, na.rm = T)
sf2 <- median(df2$R_SRR1039509, na.rm = T)
sf3 <- median(df2$R_SRR1039512, na.rm = T)
sf4 <- median(df2$R_SRR1039513, na.rm = T)
sfs <- c(sf1, sf2, sf3, sf4)
> sfs
[1] 1.1241342 0.9892460 1.2740402 0.7208434

这里展示一个样本基因比率分布图，红线即为比率的中值

hist(df2$R_SRR1039508, main='sample / pseudo-reference sample', xlab='')
abline(v=median(df2$R_SRR1039508), col='red')

Step 4: 计算normalized coutns

df1_norm <- t(t(df1[,-5])/sfs)
knitr::kable(head(df1_norm))

我们可以与DESeq2的sizefactor对比一下，看看我们算的是否正确

# using deseq2
library(DESeq2)

dds <- DESeqDataSetFromMatrix(countData = df1[,-5], colData = colData(airway)[1:4,], design=~dex)
dds <- estimateSizeFactors(dds)
DESeq_SFs <- dds$sizeFactor

sfs == DESeq_SFs
SRR1039508 SRR1039509 SRR1039512 SRR1039513 
      TRUE       TRUE       TRUE       TRUE

以上就是DESeq2计算size factors，以及校正counts的方法。事实上DESeq2在差异分析时并不会真正地用normalized counts进行分析，而是使用raw counts，并将size factor作为其广义线性回归模型（generalized linear model）的一个因子进行考虑。

完。

ref:

https://hbctraining.github.io/DGE_workshop/lessons/02_DGE_count_normalization.html