Bulk RNAseq analysis

set.seed(12345)

pacman::p_load(tidyverse, SummarizedExperiment, airway, org.Hs.eg.db,
               RColorBrewer,EnhancedVolcano, pheatmap, AnnotationDbi, DESeq2, styler, clusterProfiler, remotes)

Installing airway

if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install("airway")

data(airway)

class(airway)

str(airway)

## Access metadata
meta <- airway %>% 
  colData() %>% ## To access metadata information
  as.data.frame()
head(meta)

meta <- meta[,c(2,3)]
## Suggest a different approach

meta$dex <- gsub('trt', 'treated', meta$dex)
## Suggest a different approach

meta$dex <- gsub('untrt', 'untreated', meta$dex)

## Suggest a different approach
names(meta) <- c('cellLine', 'dexamethasone')

write.table(meta, file = "meta.csv", sep = ',', 
            col.names = T, row.names = T, 
            quote = F)

counts <- assay(airway)
write.table(counts, file = "counts_data.csv", sep = ',', col.names = T, row.names = T, quote = F)

counts <- read.csv("counts_data.csv")

counts[1:6, 1:6]

meta$dexamethasone <- relevel(factor(meta$dexamethasone), ref = "untreated")

x <- SummarizedExperiment(list(counts=as.matrix(counts)), 
                                colData = meta)
x

symbols <- mapIds(
    org.Hs.eg.db,
    keys = rownames(x),
    column = "SYMBOL",
    keytype = "ENSEMBL",
    multiVals = "first"
)

rowData(x)$symbol <- symbols

x

ens_ids <- rownames(x)
ens_ids_clean <- sub("\\..*$", "", ens_ids)

gene_symbols <- mapIds(
  org.Hs.eg.db,
  keys = ens_ids_clean,
  keytype = "ENSEMBL",
  column = "SYMBOL",
  multiVals = "first"
)

gene_names <- mapIds(
  org.Hs.eg.db,
  keys = ens_ids_clean,
  keytype = "ENSEMBL",
  column = "GENENAME",
  multiVals = "first"
)

rowData(x)$ensembl_id <- ens_ids_clean
rowData(x)$symbol <- gene_symbols
rowData(x)$gene_name <- gene_names

x

# making sure the row names in colData matches to column names in counts_data
all(colnames(counts) %in% rownames(meta))
# are they in the same order?
all(colnames(counts) == rownames(meta))

# Step 2: construct a DESeqDataSet object ----------
dds <- DESeqDataSetFromMatrix(countData = counts,
                       colData = meta,
                       design = ~ dexamethasone)
dds
# pre-filtering: removing rows with low gene counts, keeping rows that have at least 10 reads total
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]

dds
# set the factor level
dds$dexamethasone <- relevel(dds$dexamethasone, ref = "untreated")

# Step 3: Run DESeq ----------------------
dds <- DESeq(dds)

res <- results(dds)

str(res)

summary(res)

# contrasts
resultsNames(dds)
# e.g.: treated_4hrs, treated_8hrs, untreated

results(dds, contrast = c("dexamethasone", "treated", "untreated"))
# MA plot
plotMA(res)

vsd <- vst(dds, blind = FALSE)
mat <- assay(vsd)

plotPCA(vsd, intgroup = "dexamethasone")

pca_data <- plotPCA(vsd, intgroup = "dexamethasone", returnData = TRUE)

percent_var <- round(100 * attr(pca_data, "percentVar"))
ggplot(pca_data, aes(PC1, PC2, color = dexamethasone, label = name)) +
geom_point(size = 4) +
geom_text(vjust = -1) +
xlab(paste0("PC1: ", percent_var[1], "% variance")) +
ylab(paste0("PC2: ", percent_var[2], "% variance")) +
theme_minimal()

annot <- data.frame(ensembl_id = rowData(x)$ensembl_id,
                    symbol = rowData(x)$symbol,
                    gene_name = rowData(x)$gene_name)

head(annot)

head(res, 4)

res_df <- as.data.frame(res) %>%
  rownames_to_column("ensembl_id") %>%
  mutate(ensembl_id = sub("\\..*$", "", ensembl_id)) %>%
  left_join(annot, by = "ensembl_id")

VOLCANO_Plot <- function(df, units = "in", height = 16, width = 16,
filename, dpi = 600,
cols = c("UP" = "firebrick3", "DOWN" = "steelblue", "ns" = "grey"),
ylim = c("",""),
xlim = c(-2.5,2.5)){
vp1 <- EnhancedVolcano(df,
lab = df$symbol,
x = 'log2FoldChange',
y = 'padj', pCutoff = 0.05, FCcutoff = 0.5)
print(vp1)
# ggsave(filename, units = units, height = height, width = width, dpi = dpi)
return(vp1)
}

VOLCANO_Plot(df = res_df)

sig_res <- subset(res_df, 
                  padj < 0.05 & 
                    abs(log2FoldChange) > 1)
head(sig_res)

nrow(sig_res)

head(sig_res[order(-sig_res$log2FoldChange), ])

head(sig_res[order(sig_res$log2FoldChange), ])

res_df$significant <- ifelse(res_df$padj < 0.05 & 
                               abs(res_df$log2FoldChange) > 1, "yes", "no")

ggplot(res_df, aes(x = log2FoldChange, y = -log10(padj), color = significant)) +
geom_point(alpha = 0.7) +
theme_minimal() +
xlab("log2 fold change") +
ylab("-log10 adjusted p-value")

res_df <- res %>% 
  as.data.frame() %>% 
  tibble::rownames_to_column("ensembl_id") %>% 
  mutate(ensembl_id = sub("\\..*$", "", ensembl_id)) %>% 
  left_join(annot, by = "ensembl_id") %>% 
  filter(padj < 0.05, 
         abs(log2FoldChange) > 1) %>% 
  mutate(
    DE = case_when(
      log2FoldChange > 1 ~ "UP",
      log2FoldChange < -1 ~ "DOWN"
    )
  )

head(res_df)

symbol_to_ens <- setNames(
  rowData(x)$ensembl_id,
  rowData(x)$symbol
)

plotCounts(dds, gene = "PDK4", intgroup = "dexamethasone", main = "PDK4")

plotCounts(dds, gene = symbol_to_ens["PDK4"], intgroup = "dexamethasone", main = "PDK4")

plotCounts(dds, gene = symbol_to_ens["TP53"], intgroup = "dexamethasone", main = "TP53")

Construct a loop to plot 10 genes

library(dplyr)

top_up <- res_df %>%
  filter(DE == "UP") %>%
  arrange(padj) %>%
  slice_head(n = 20)

top_down <- res_df %>%
  filter(DE == "DOWN") %>%
  arrange(padj) %>%
  slice_head(n = 20)

heatmap_df <- bind_rows(top_up, top_down)

mat_sub <- mat[heatmap_df$ensembl_id, , drop = FALSE]

gene_labels <- heatmap_df$symbol
gene_labels[is.na(gene_labels) | gene_labels == ""] <- heatmap_df$ensembl_id[is.na(gene_labels) | gene_labels == ""]
gene_labels <- make.unique(gene_labels)

rownames(mat_sub) <- gene_labels

mat_scaled <- t(scale(t(mat_sub)))
mat_scaled[is.na(mat_scaled)] <- 0

annotation_col <- data.frame(
  dexamethasone = colData(dds)$dexamethasone
)

rownames(annotation_col) <- colnames(mat_scaled)

pheatmap(
  mat_scaled,
  annotation_col = annotation_col,
  cluster_rows = TRUE,
  cluster_cols = TRUE,
  show_rownames = TRUE,
  show_colnames = TRUE,
  fontsize_row = 8,
  fontsize_col = 10,
  border_color = NA,
  color = colorRampPalette(rev(brewer.pal(n = 11, name = "RdBu")))(100),
  main = "Top differentially expressed genes"
)

GO terms

GO_res <- enrichGO(
  gene = res_df$symbol,
  OrgDb = "org.Hs.eg.db",
  keyType = "SYMBOL",
  ont = "ALL",
  pAdjustMethod = "BH",
  qvalueCutoff = 0.05,
  pvalueCutoff = 0.05,
  readable = TRUE
)

str(GO_res)

dotplot(GO_res)

Up_reg <- res_df %>% filter(DE == "UP")
Down_reg <- res_df %>% filter(DE == "DOWN")

GO_res_up <- enrichGO(
  gene = Up_reg$symbol,
  OrgDb = "org.Hs.eg.db",
  keyType = "SYMBOL",
  ont = "ALL",
  pAdjustMethod = "BH",
  qvalueCutoff = 0.05,
  pvalueCutoff = 0.05,
  readable = TRUE
)

dotplot(GO_res_up)

GO_res_down <- enrichGO(
  gene = Down_reg$symbol,
  OrgDb = "org.Hs.eg.db",
  keyType = "SYMBOL",
  ont = "ALL",
  pAdjustMethod = "BH",
  qvalueCutoff = 0.05,
  pvalueCutoff = 0.05,
  readable = TRUE
)

dotplot(GO_res_down)

More visualization tools

enrichplot::cnetplot(GO_res_up)

GO_res_up <- enrichplot::pairwise_termsim(GO_res_up)
enrichplot::emapplot(GO_res_up)

enrichplot::emapplot_cluster(GO_res_up)

enrichplot::upsetplot(GO_res_up)
1

enrichplot::ssplot(GO_res_up)
ggupset, 1
1

enrichplot::treeplot(GO_res_up, showCategory = 10)