Normalize Counts Matrix in DESeq2
I followed Josh Starmer's YouTube video on how DESeq2 normalizes the raw counts matrix, and attempted to reproduce the steps in R. This normalization procedure accounts for both the library size (total no. of reads in each sample) as well as the library composition (enables us to compare between different tissues). Following are the steps -
- Convert raw counts to loge(raw counts)
- Average each row of the log-transformed counts matrix
- Remove genes with infinity values from the matrix and averages
- For each sample, subtract the averages from the loge(raw counts)
- For each sample, compute median of the resulting values
- Raise e to the power of -(median values) to obtain size factors for each sample
- Divide the raw counts in step 1 by size factors
R implementation -
library(matrixStats)
# Set raw counts matrix
counts_mat <- data.frame (
Sample_1 = c(0, 10, 35, 45, 5),
Sample_2 = c(24, 20, 70, 90, 10),
Sample_3 = c(12, 30, 42, 5, 70)
)
row.names(counts_mat) <- c("Fam189a2", "RGD1562699",
"Rnf115", "Mrnip", "AABR070585111")
counts_mat
# Sample_1 Sample_2 Sample_3
# Fam189a2 0 24 12
# RGD1562699 10 20 30
# Rnf115 35 70 42
# Mrnip 45 90 5
# AABR070585111 5 10 70
# log transform the raw counts
dat <- log(counts_mat)
dat
# Sample_1 Sample_2 Sample_3
# Fam189a2 -Inf 3.178054 2.484907
# RGD1562699 2.302585 2.995732 3.401197
# Rnf115 3.555348 4.248495 3.737670
# Mrnip 3.806662 4.499810 1.609438
# AABR070585111 1.609438 2.302585 4.248495
# Compute row averages
rowsum_values <- rowMeans(dat)
rowsum_values
# Fam189a2 RGD1562699 Rnf115 Mrnip AABR070585111
# -Inf 2.899838 3.847171 3.305303 2.720173
# Remove genes with infinity values
dat <- dat[!is.infinite(rowsum_values),]
dat
# Sample_1 Sample_2 Sample_3
# RGD1562699 2.302585 2.995732 3.401197
# Rnf115 3.555348 4.248495 3.737670
# Mrnip 3.806662 4.499810 1.609438
# AABR070585111 1.609438 2.302585 4.248495
# Remove genes with infinity values
rowsum_values <- rowsum_values[!is.infinite(rowsum_values)]
rowsum_values
# RGD1562699 Rnf115 Mrnip AABR070585111
# 2.899838 3.847171 3.305303 2.720173
# Subtract the row averages from each column of counts
dat <- apply(dat, 2, function(x) x - rowsum_values)
dat
# Sample_1 Sample_2 Sample_3
# RGD1562699 -0.5972532 0.09589402 0.5013591
# Rnf115 -0.2918229 0.40132427 -0.1095014
# Mrnip 0.5013591 1.19450631 -1.6958654
# AABR070585111 -1.1107348 -0.41758766 1.5283225
# Compute median of resulting values for each sample
median_vals <- colMedians(as.matrix(dat))
median_vals
# [1] -0.4445380 0.2486091 0.1959289
# Calculate size factors using median
size_factor <- exp(median_vals)
size_factor
# [1] 0.6411204 1.2822408 1.2164404
# Divide raw counts by size factors
normalized_counts <- t(apply(counts_mat, 1,
function(x) x/size_factor))
normalized_counts
# Sample_1 Sample_2 Sample_3
# Fam189a2 0.000000 18.717234 9.864848
# RGD1562699 15.597695 15.597695 24.662121
# Rnf115 54.591931 54.591931 34.526969
# Mrnip 70.189626 70.189626 4.110353
# AABR070585111 7.798847 7.798847 57.544948
In order to confirm whether the steps are correct, I used the in-built function of DESeq2 -
library(DESeq2)
metadata <- data.frame(
Samples=c("Sample_1", "Sample_2", "Sample_3"),
group=c("ctrl", "ctrl", "treat")
)
row.names(metadata) <- metadata$Samples
metadata
# Samples group
# Sample_1 Sample_1 ctrl
# Sample_2 Sample_2 ctrl
# Sample_3 Sample_3 treat
dds <- DESeqDataSetFromMatrix(countData = counts_mat,
colData = metadata,
design = ~ group)
dds <- DESeq(dds)
normalized_counts <- counts(dds, normalized=T)
normalized_counts
# Sample_1 Sample_2 Sample_3
# Fam189a2 0.000000 18.717234 9.864848
# RGD1562699 15.597695 15.597695 24.662121
# Rnf115 54.591931 54.591931 34.526969
# Mrnip 70.189626 70.189626 4.110353
# AABR070585111 7.798847 7.798847 57.544948