Workshop: DESeq2 pairwise comparisons

A DESeq2 teaching workflow for pairwise treatment comparisons, with data import, gene filtering, VST visualization, differential-expression testing, and plotting.

Teaching document: DESeq2_pairwise_demo-20231212
Output: HTML document
Author: CTLin, ICOB Data directory

/home/workshop/demo/02_DE/

Rmd

~/Shared/518_projects/workshop_20231213/DESeq2_pairwise_demo.Rmd

DESeq2 pairwise comparison

Library

library(DESeq2)
library(readr)
library(tibble)
library(pheatmap)
library(dplyr)
library(ggplot2)
source("/home/workshop/demo/02_DE/Hmap_PatternCluster_func.R")
source("/home/workshop/demo/02_DE/prepare_df_for_certainGene_pattern_func.R")

# run mkdir DESeq2_output/R_output in the terminal

Import data

# import the tab-delimited text file of count matrix
count_matrix <- read_delim("~/count_matrix.txt",
                           delim = "\t",
                           escape_double = FALSE,
                           trim_ws = TRUE)

count_matrix <- column_to_rownames(count_matrix, "Geneid")

# import the metadata
metadata <- read_delim("/home/workshop/demo/02_DE/metadata_pairwise.txt",
                       delim = "\t",
                       escape_double = FALSE,
                       trim_ws = TRUE)


metadata$rep <- factor(metadata$rep, levels = c("r1","r2","r3"))
metadata$treatment <- factor(metadata$treatment, levels = c("ctrl","drug"))

# eliminate rows with NA
table(is.na(count_matrix))
# no NA

# check if there are duplicated gene id?
length(unique(rownames(count_matrix)))
# no

# the read count cut-off of the low-expression genes
low_exp_read_count = 10
# minimum samples that express this gene (總共有24個 sample)
minimum_exp_sample = 3

dds object & gene filtering

dds_pw <- DESeqDataSetFromMatrix(countData= count_matrix,
                                colData= metadata,
                                design= ~treatment+rep) # batch effect

keep <- rowSums(counts(dds_pw)>low_exp_read_count) >= minimum_exp_sample
dds_pw <- dds_pw[keep, ] # 17539 genes

VST transformation

# Output cross-sample normalization counts
# Size factor
dds_pw<-estimateSizeFactors(dds_pw)

normalized_counts <- counts(dds_pw, normalized=TRUE)
normalized_counts <- tibble::rownames_to_column(as.data.frame(normalized_counts), "Gene_id")

write.csv(normalized_counts, file = "DESeq2_output/R_output/DESeq2_cross-sample_normalized_counts.csv")





# vst transformation
vst <- vst(dds_pw, blind=FALSE)
plotPCA(vst,intgroup="treatment")
plotPCA(vst,intgroup="rep")
vst_mat <- assay(vst)

write.csv(vst_mat, file = "DESeq2_output/R_output/DESeq2_cross-sample_normalized_vst-counts.csv")

dds_wald <- DESeq(dds_pw,
              test = "Wald",
              fitType = "parametric",
              betaPrior = F) #  whether or not to put a zero-mean normal prior on the non-intercept coefficients
plotDispEsts(dds_wald)

DEG between treatments

result_drug_ctrl <-results(dds_wald,
                      contrast=c("treatment","drug","ctrl"),      # base factor level 放後面
                      alpha = 0.05,
                      independentFiltering = TRUE,
                      pAdjustMethod = "fdr")
# before shrinkage
result_drug_ctrl

# afger shrinkage
result_drug_ctrl_shrink_ashr <- lfcShrink(dds_wald,
                                      type="ashr",
                                      res = result_drug_ctrl)

plotMA(result_drug_ctrl)
plotMA(result_drug_ctrl_shrink_ashr)

subset DEG

result_drug_ctrl_shrink_ashr_df <- as.data.frame(result_drug_ctrl_shrink_ashr)
result_drug_ctrl_shrink_ashr_df <- result_drug_ctrl_shrink_ashr_df[order(result_drug_ctrl_shrink_ashr_df$padj),]


deg_drug_ctrl_padj01 <- subset(result_drug_ctrl_shrink_ashr_df,
                                 (result_drug_ctrl_shrink_ashr_df$padj < 0.05 &
                                  abs(result_drug_ctrl_shrink_ashr_df$log2FoldChange) > 1))




deg_drug_ctrl_padj01 <- tibble::rownames_to_column(deg_drug_ctrl_padj01, "Gene_id")
# -> 166 DEGs

write.csv(deg_drug_ctrl_padj01, file = "DESeq2_output/R_output/DEGs_drug_ctrl_padj<0.01.csv")

Visualization

result_drug_ctrl_shrink_ashr_df <- na.omit(result_drug_ctrl_shrink_ashr_df)
write.csv(result_drug_ctrl_shrink_ashr_df, file = "DESeq2_output/R_output/DESeq2_drug_ctrl_Shrinklog2FC.csv")





# plot
result_drug_ctrl_shrink_ashr_df <- as.data.frame(result_drug_ctrl_shrink_ashr)
result_drug_ctrl_shrink_ashr_df$shape <- ifelse(result_drug_ctrl_shrink_ashr_df$padj<10^-30,2,20)
result_drug_ctrl_shrink_ashr_df$padj[result_drug_ctrl_shrink_ashr_df$padj<10^-30] <- 10^-30
pdf(file = "DESeq2_output/R_output/volcano_drug.vs.ctrl.pdf", height = 5, width = 7)
with(result_drug_ctrl_shrink_ashr_df,
     plot(log2FoldChange,
          -log10(padj),
          pch=shape,
          main=""))

# Add colored points: blue if log2FC < -1 and padj < 0.05
# red if log2FC > 1 and padj<0.05)
with(subset(result_drug_ctrl_shrink_ashr_df,
            padj<.05 & log2FoldChange < -1),
     points(log2FoldChange,
            -log10(padj),
            pch=shape,
            col="blue"))

with(subset(result_drug_ctrl_shrink_ashr_df,
            padj<.05 & log2FoldChange > 1),
     points(log2FoldChange,
            -log10(padj),
            pch=shape,
            col="red"))

legend("topright",
       legend=c("log2FC < -1 and padj < 0.05",
                "log2FC > 1 and padj < 0.05"),
       col=c("blue", "red"),
       cex = 0.6,
       pch = 20)
dev.off()

← Back to projects