2 min read

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 -

  1. Convert raw counts to loge(raw counts)
  2. Average each row of the log-transformed counts matrix
  3. Remove genes with infinity values from the matrix and averages
  4. For each sample, subtract the averages from the loge(raw counts)
  5. For each sample, compute median of the resulting values
  6. Raise e to the power of -(median values) to obtain size factors for each sample
  7. 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