Analysis of differential gene expression in DESeq in R Notebook

The data used here is from jci research article. Study is about Silencing (KO) the gene SMOC2 protects from kidney fibrosis. so the data is of this format: WT_norm, WT_fibrosis, SMOC2_overexpression_normal, SMOC2_overexpression_fibrosis

#setwd("D:/shilpa/RNAseq/datacamp")
library(DESeq2)
library(RColorBrewer)
library(pheatmap)
library(tidyverse)

raw_count_matrix <- read.csv('fibrosis_smoc2_rawcounts_unordered.csv')
ggplot(raw_count_matrix) + geom_histogram(aes(x=smoc2_normal1), stat="bin", bins= 200) + xlab("Raw expression counts") + ylab("Number of genes") + xlim(0, 25000) + ylim(0, 5000)

The gene SMOC2, bound to fate, its expression studies we now undertake, with this spirit now we shall see how our data looks like

head(raw_count_matrix)

str(raw_count_matrix)

'data.frame':   47729 obs. of  8 variables:
 $ X              : chr  "ENSMUSG00000102693" "ENSMUSG00000064842" "ENSMUSG00000051951" "ENSMUSG00000102851" ...
 $ smoc2_fibrosis1: int  0 0 72 0 0 0 0 0 0 1 ...
 $ smoc2_fibrosis4: int  0 0 30 0 0 0 0 0 0 1 ...
 $ smoc2_normal1  : int  0 0 0 0 1 0 0 0 0 1 ...
 $ smoc2_normal3  : int  0 0 3 0 0 0 0 0 0 0 ...
 $ smoc2_fibrosis3: int  0 0 36 0 0 0 0 0 0 1 ...
 $ smoc2_normal4  : int  0 0 1 0 0 0 0 0 0 0 ...
 $ smoc2_fibrosis2: int  0 0 51 0 0 0 0 0 0 1 ...

With a steady hand and keen eye, we shall forge the metadata to construct our DESeq2 object!

genotype <- c("smoc2_oe", "smoc2_oe", "smoc2_oe", "smoc2_oe", "smoc2_oe", "smoc2_oe", "smoc2_oe")
condition <- c("fibrosis", "fibrosis", "fibrosis", "fibrosis", "normal", "normal", "normal")
metadata <- data.frame(genotype, condition)

# Assign the row names of the data frame
rownames(metadata) <- c("smoc2_fibrosis1", "smoc2_fibrosis2", "smoc2_fibrosis3", "smoc2_fibrosis4", "smoc2_normal1", "smoc2_normal3", "smoc2_normal4")

metadata

NA

While there are many packages for DE analysis, we are going to use DESeq2. It has one of more popular ones and its vignetter is pretty helpful. DESeq2 models the gene expression (count data) by negetive binomial distribution.

The raw data from ncbi GEO has already been preprocessed (QC steps & counts).

What is a count matrix: count matrix represents number of reads matching the exons of each gene. Barcodes (for each cell) and UMI (for each gene) information is used to generate this matrix. A little bit about count matrix:

• The negative binomial model commonly used to model RNA-seq count data

• it is discrete and it starts with zero.

• Highest frequency near zero i.e. many genes have low number of counts

• The expression range is between zero to inf (no max limit)

• For discrete numbers i.e. counts one could also use poisson but there is too much variation in the data than poisson can handle

NOTE: DESeq2 model internally corrects for the library size so transformed or normalized values such as counts scaled by library size should not be used as input.

Bringing in data for DESeq2: we need to bring the samples from metadata and countdata in the same order

#metadata <- read.csv("metadata.csv")

#all(colnames(raw_count_matrix) == rownames(metadata))

# Use the match() function to reorder the columns of the raw counts
reorder_idx <- match(rownames(metadata), colnames(raw_count_matrix))

# Reorder the columns of the count data
reordered_raw_count_matrix <- raw_count_matrix[ , reorder_idx]
# Create a DESeq2 object (DESeqDataSet dds)
dds <- DESeqDataSetFromMatrix(countData = reordered_raw_count_matrix,
                                    colData = metadata,
                                    design = ~ condition)

Warning: some variables in design formula are characters, converting to factors

# design formula is used to estimate the dispersions and to estimate the log2 fold changes of the model

dds

class: DESeqDataSet 
dim: 47729 7 
metadata(1): version
assays(1): counts
rownames: NULL
rowData names(0):
colnames(7): smoc2_fibrosis1 smoc2_fibrosis2 ... smoc2_normal3 smoc2_normal4
colData names(2): genotype condition

There are 4 ways of constructing DESeqDataSet depending on which pipeline is used upstream. We have count matrix here so we use ‘countData’ while constructing dds. Now we can explore the data with DESeq2 functions and perform the differential expression analysis.

NORMALIZING COUNTS WITH DESeq2 -> normalized for library size while accounting for library composition (library size is the total number of gene counts per sample)

dds <- estimateSizeFactors(dds)
normalized_counts <- counts(dds, normalized=TRUE)

We will be using the normalized counts to explore similarities in gene expression between each of our samples. To do this we use clustering after using dimensional reduction. Before clustering, we need to log transform the normalized counts to improve the visualization of the clustering and hence our understanding. we can use DESeq2’s vst() function. vst stands for Variance stabilizing transformation. The transformed data should be approximated variance stabilized and also includes correction for size factors or normalization factors. The transformed data is on the log2 scale for large counts. The blind=TRUE argument specifies that the transformation should be blind to the sample information given in the design formula; this argument should be specified when performing quality assessment.

vsd_trans <- vst(dds, blind = TRUE)
vsd_trans

class: DESeqTransform 
dim: 47729 7 
metadata(1): version
assays(1): ''
rownames: NULL
rowData names(4): baseMean baseVar allZero dispFit
colnames(7): smoc2_fibrosis1 smoc2_fibrosis2 ... smoc2_normal3 smoc2_normal4
colData names(3): genotype condition sizeFactor

From vsd, we need to extract the VST-trasformed normalized counts as a matrix

vsd_mat <- assay(vsd_trans)

Then, we need to compute the pairwise correlation values between each pair of samples

vsd_cor <- cor(vsd_mat)

Now we do clustering. We choose Hierarchical clustering here.

pheatmap(vsd_cor, annotation = select(metadata, condition)) 

This heatmap shows us that the biological replicates cluster together (as expected, phew…). This is a good sign because, very likely our differentially expressed genes must be driving this separation. Also we do not see any outliers.

We will crosscheck this with PCA (QC step 2)

plotPCA(vsd_trans, intgroup="condition")

Now that we have explored the quality of the samples, checked for outliers, we can go ahead and find which genes have significant differences in expression between the normal and fibrosis sample groups.

At this stage we do not need to re-create the object since we have not removed samples or found additional sources of variation during QC steps.

dds_run <- DESeq(dds)

using pre-existing size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

dds_run

class: DESeqDataSet 
dim: 47729 7 
metadata(1): version
assays(4): counts mu H cooks
rownames: NULL
rowData names(22): baseMean baseVar ... deviance maxCooks
colnames(7): smoc2_fibrosis1 smoc2_fibrosis2 ... smoc2_normal3 smoc2_normal4
colData names(3): genotype condition sizeFactor

mean_counts <- apply(raw_count_matrix[,1:3], 1, mean)

Dispersion estimates are used to model the raw counts; if the dispersions don’t follow the assumptions made by DESeq2, then the variation in the data could be poorly estimated and outout of which genes are differentially expressed could be wrong. The assumptions made by DESeq2 are:

• The counts data fit to the negative binomial distribution
• The dispersions should generally decrease with increase in the mean

plotDispEsts(dds_run)

From the dispersion plot we can see that fit of our data to the model is pretty good. Gene-est is estimated dispersion of the genes.

Plotting the dispersion estimates is a useful diagnostic. The dispersion plot which we see here is a typical one with the final estimates shrunk from the gene-wise estimates towards the fitted estimates. Some gene-wise estimates are flagged (open circles) as outliers and not shrunk towards the fitted value.

res0.01 <- results(dds_run,alpha = 0.01)
summary(res0.01)

out of 29556 with nonzero total read count
adjusted p-value < 0.01
LFC > 0 (up)       : 4534, 15%
LFC < 0 (down)     : 4798, 16%
outliers [1]       : 15, 0.051%
low counts [2]     : 8317, 28%
(mean count < 2)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

# contrasts
resultsNames(dds_run)

[1] "Intercept"                    "condition_normal_vs_fibrosis"

# MA plot
plotMA(res0.01)