理论 | DESeq2 normalized counts --
2021-10-16 本文已影响0人
尘世中一个迷途小书僮
在做基因差异表达分析时，经常会用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))
| SRR1039508 | SRR1039509 | SRR1039512 | SRR1039513 | pseudo | |
|---|---|---|---|---|---|
| ENSG00000260166 | 0 | 0 | 0 | 0 | 0.000 |
| ENSG00000266931 | 0 | 0 | 0 | 0 | 0.000 |
| ENSG00000104774 | 1939 | 1906 | 2294 | 1205 | 1787.806 |
| ENSG00000267583 | 0 | 1 | 0 | 0 | 0.000 |
| ENSG00000227581 | 1 | 0 | 0 | 0 | 0.000 |
| ENSG00000227317 | 0 | 0 | 0 | 0 | 0.000 |
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为例，写出具体计算过程：
| Row.names | SRR1039508 | SRR1039509 | SRR1039512 | SRR1039513 | pseudo | R_SRR1039508 | R_SRR1039509 | R_SRR1039512 | R_SRR1039513 |
|---|---|---|---|---|---|---|---|---|---|
| ENSG00000000460 | 60 | 55 | 40 | 35 | 46.36182 | 60/46.36182=1.2941683 | 55/46.36182=1.1863209 | 40/46.36182=0.8627789 | 35/46.36182=0.7549315 |
| ENSG00000000971 | 3251 | 3679 | 6177 | 4252 | 4209.97390 | 0.7722138 | 0.8738772 | 1.4672300 | 1.0099825 |
| ENSG00000002919 | 234 | 246 | 298 | 211 | 245.28016 | 0.9540111 | 1.0029348 | 1.2149372 | 0.8602408 |
| ENSG00000003436 | 4373 | 4183 | 5794 | 3230 | 4301.42236 | 1.0166405 | 0.9724690 | 1.3469963 | 0.7509144 |
| ENSG00000004660 | 95 | 55 | 59 | 35 | 57.31279 | 1.6575705 | 0.9596461 | 1.0294385 | 0.6106839 |
| ENSG00000004766 | 492 | 413 | 546 | 337 | 439.72831 | 1.1188727 | 0.9392163 | 1.2416758 | 0.7663823 |
由于大部分基因不会发生差异表达，所以样本内的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))
| SRR1039508 | SRR1039509 | SRR1039512 | SRR1039513 | |
|---|---|---|---|---|
| ENSG00000260166 | 0.0000000 | 0.000000 | 0.000 | 0.000 |
| ENSG00000266931 | 0.0000000 | 0.000000 | 0.000 | 0.000 |
| ENSG00000104774 | 1724.8830895 | 1926.719981 | 1800.571 | 1671.653 |
| ENSG00000267583 | 0.0000000 | 1.010871 | 0.000 | 0.000 |
| ENSG00000227581 | 0.8895735 | 0.000000 | 0.000 | 0.000 |
| ENSG00000227317 | 0.0000000 | 0.000000 | 0.000 | 0.000 |
我们可以与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
