This DESeq2 teaching workflow uses a likelihood ratio test to compare a treatment-by-timepoint model with a reduced model. Data import, quality checks, testing, and visualization follow in analysis order. Data directory

/home/workshop/demo/02_DE/

Rmd

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

Teaching document: Workshop_DESeq2-LRT-demo_20231210
Output: HTML document
Author: CTLin

Likelihood ratio test (LRT)

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/LRT in the terminal

Import data

# input the tab-delimited text file of count matrix
count_matrix <- read_delim("/home/workshop/demo/03_LRT/LRT_demo_matrix.txt",
                           delim = "\\t",
                           escape_double = FALSE,
                           trim_ws = TRUE)

colnames(count_matrix) = c("Gene_id",
                           "ctrl_t1.r1", "ctrl_t1.r2", "ctrl_t1.r3",
                           "drug_t1.r1", "drug_t1.r2", "drug_t1.r3",
                           "ctrl_t2.r1", "ctrl_t2.r2", "ctrl_t2.r3",
                           "drug_t2.r1", "drug_t2.r2", "drug_t2.r3",
                           "ctrl_t3.r1", "ctrl_t3.r2", "ctrl_t3.r3",
                           "drug_t3.r1", "drug_t3.r2", "drug_t3.r3")
count_matrix <- column_to_rownames(count_matrix, "Gene_id")


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


metadata$group <- factor(metadata$group, levels = c("ctrl_t1","ctrl_t2","ctrl_t3",
                                                    "drug_t1","drug_t2","drug_t3"))
metadata$rep <- factor(metadata$rep)
metadata$treatment <- factor(metadata$treatment)
metadata$timepoint <- factor(metadata$timepoint)

Omit low quality, duplicated genes

# 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 = 9

Create dds object

dds_time <- DESeqDataSetFromMatrix(countData= count_matrix,
                                   colData= metadata,
                                   design= ~treatment+timepoint+treatment:timepoint+rep)

keep <- rowSums(counts(dds_time)>low_exp_read_count) >= minimum_exp_sample
dds_time <- dds_time[keep, ]

# 69573 -> 16393

Output cross-sample normalization counts

# Size factor
dds_time<-estimateSizeFactors(dds_time)
normalized_counts <- counts(dds_time, normalized=TRUE)
normalized_counts <- tibble::rownames_to_column(as.data.frame(normalized_counts), "Gene_id")

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

Check PCA (vst counts) & output cross-sample normalization counts with vst transformation

The transformed counts are used for PCA and visualization; the differential-expression model below uses the count-based DESeq2 object.

vst <- vst(dds_time, blind=FALSE)

plotPCA(vst,intgroup="group")
plotPCA(vst,intgroup="treatment")
plotPCA(vst,intgroup="timepoint")
plotPCA(vst,intgroup="rep")

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

DESeq2 LRT analysis

We are interested in the “treatment:timepoint” (stands for “the difference in the effect of treatment over time”). Therefore, we compare the full model with the “reduced” model.

dds_time_LRT <- DESeq(dds_time,
                      test="LRT",
                      reduced = ~treatment+timepoint+rep)


plotDispEsts(dds_time_LRT) # How to interpret the plotDispEsts

DEG analysis: Output summary table

result_time <- results(dds_time_LRT,
                        pAdjustMethod = "fdr",
                        independentFiltering = TRUE,
                        alpha = 0.05)

result_time_df <- as.data.frame(result_time)
result_time_df <- result_time_df[order(result_time_df$padj),]
write.csv(result_time_df, file = "DESeq2_output/R_output/LRT/DESeq2_LRT-DEG-result.csv")

Check outliers and genes with low mean normalized count

# Check outlier: pvalue and padj = NA
result_time_df[is.na(result_time_df$pvalue) & is.na(result_time_df$padj),]
# -> 有 0 個 gene 是 outlier

# Check how many genes were filtered by automatic independent filtering, for having a low mean normalized count: only padj = NA
result_time_df[is.na(result_time_df$padj) & !is.na(result_time_df$pvalue),]
# -> 有 0 個 gene 被 independent filtered

metadata(result_time)$filterThreshold
# DESeq2 has filtered out genes with a mean count of 7.072094 ->
# 0% of genes were filtered out of the analysis

Subset DEGs

you can also subset DEGs with the criterias you want

deg_padj01 <- subset(result_time_df, (result_time_df$padj < 1e-20))


deg_padj01 <- tibble::rownames_to_column(deg_padj01, "Gene_id")
# -> 443 DEGs




write.csv(deg_padj01, file = "DESeq2_output/R_output/LRT/DESeq2_LRT_DEG_padj<1e-20.csv")

Clustering the DEGs with similar expression profile along timepoints

# subset vst with only DEGs
vst_sig <- vst_mat[rownames(vst_mat) %in% deg_padj01$Gene_id,]


# visualization
pdf(file = "DESeq2_output/R_output/LRT/hmap_pattern_LRT_DEG_padj<1e-20.pdf")
heatmap_and_clusterPattern(vst_matrix = vst_sig,
                           metadata = metadata,
                           heatmap_title =
                             "DEGs identified by LRT analysis \\nvst_counts_scaled (scaled for each gene), break: 0 \u00b1 2,\n euclidean, ward.D2",
                           value_break = 2,
                           PatternCluster = T,
                           group_rep_delim = ".", # delim between group and replicate number
                           treatment_timepoint_delim = "_",
                           PatternCluster_height = 0.25, # 25% height
                           PatternCluster_csv.file = "DESeq2_output/R_output/LRT/cluster_results.csv", # cluster information
                           PatternCluster_txt.file = "DESeq2_output/R_output/LRT/cluster_results.txt",
                           PatternCluster_title =
                           "DEGs identified by LRT analysis \\nvst_counts_scaled (scaled for each gene), break: 0 \u00b1 2,\n euclidean, ward.D2, cutree: 25%")
dev.off()

Plot the expression profile for certain gene

eg, plot 5 genes with the smallest padj values

x=5 #  plot 5 genes' profiles with smallest padj values in result_df

for (i in 1:x) {
  assign(paste0("g_",i),
                    prepare_df_for_certainGene_pattern(normalized_counts,
                                                       deg_padj01$Gene_id[i],
                                                       group_rep_delim = ".", # delim between group and replicate
                                                       treatment_timepoint_delim = "_")
        )
  print(
  ggplot(get(paste0("g_",i)),
       aes(x = timepoint,
           y = count,
           color = treatment,
           group = treatment)) +
  geom_point() +
  stat_summary(fun.y=mean, geom="line") +
  scale_y_log10()+
  ggtitle(rownames(result_time_df)[i])
  )
}

If I want to visualize certain gene

cluster_results <- read_delim("DESeq2_output/R_output/LRT/cluster_results.txt",
                              delim = "\\t",
                              escape_double = FALSE,
                              trim_ws = TRUE)

# ex: I want to check some genes in cluster X
cluster_results[cluster_results$cluster == "cluster_2",]


gene_ID=c("XLOC_000763") # Input gene ID


df <- prepare_df_for_certainGene_pattern(normalized_counts,
                                         gene_ID,
                                         group_rep_delim = ".", # delim between group and replicate number
                                         treatment_timepoint_delim = "_")
ggplot(df,
       aes(x = timepoint,
           y = count,
           color = treatment,
           group = treatment)) +
  geom_point() +
  stat_summary(fun.y=mean, geom="line") +
  scale_y_log10()+
  ggtitle(gene_ID)




# try these non-DEGs:
#XLOC_046379
#XLOC_049384
#XLOC_028779
← Back to projects