# Normalize Counts Matrix in DESeq2

Aarthi Ramakrishnan

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)

Samples=c("Sample_1", "Sample_2", "Sample_3"),
group=c("ctrl", "ctrl", "treat")
)

#           Samples group
# Sample_1 Sample_1  ctrl
# Sample_2 Sample_2  ctrl
# Sample_3 Sample_3 treat

dds <- DESeqDataSetFromMatrix(countData = counts_mat,
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``````

Bioinformatics

## Interesting Bioinformatics Articles

Following is a collection of articles which I feel every Bioinformatician must be aware of. I will keep updating this list from time to time - 1. All biology is computational biology 2. Core services: Reward bioinformaticians 3. Importance of stupidity in scientific research

## Screen command in UNIX

Screen is a very useful command to have in your toolbox if you frequently use interactive sessions on your supercomputer logged in through a VPN. A VPN typically has a time limit, and you may get disconnected from it without any warning when you have poor internet connection. Screen program

## Types of models in DESeq2

There are 2 major types of regression models one can specify in DESeq2 to explore the raw count matrices from an RNA-seq experiment - * Mean-reference model for Factors * Regression model for Covariates Mean-reference model for Factors - Factors typically represent categorical variable such as Gender, Ethnicity, Race etc. The mean-reference