Final Data Analysis

knitr::opts_chunk$set(echo = TRUE)

Load Data

raw <- read.csv(
  "Count Table RNA Seq Data for 6703 presentation.csv",
  header = TRUE,
  stringsAsFactors = FALSE
)

# Select only count columns
counts <- raw[, c("NBF1","NBF2","NBF3",
                  "BF3hr1","BF3hr2","BF3hr3",
                  "BF24hr1","BF24hr3","BF24hr4")]

rownames(counts) <- raw$id

counts <- as.matrix(counts)
mode(counts) <- "numeric"

condition <- factor(c(rep("NBF",3),
                      rep("BF3hr",3),
                      rep("BF24hr",3)))

coldata <- data.frame(
  row.names = colnames(counts),
  condition = condition
)

DESeq2 Analysis

dds <- DESeqDataSetFromMatrix(
  countData = counts,
  colData = coldata,
  design = ~ condition
)

## converting counts to integer mode

dds <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

dds_vst1 <- varianceStabilizingTransformation(dds)

PCA

pca_df <- plotPCA(
  dds_vst1,
  ntop = 500,
  intgroup = "condition",
  returnData = TRUE
)

## using ntop=500 top features by variance

pct_var <- round(100 * attr(pca_df, "percentVar"), 1)

ggplot(pca_df, aes(x = PC1, y = PC2, color = condition)) +
  geom_point(size = 8) +
  scale_color_brewer(palette = "Dark2", name = "Feeding Treatment") +
  theme_bw(base_size = 20) +
  labs(
    x = paste0("PC1 (", pct_var[1], "%)"),
    y = paste0("PC2 (", pct_var[2], "%)")
  )

Differential Expression

NBF vs BF3hr

res_NBF.BF3 <- results(dds, contrast = c("condition", "NBF", "BF3hr"))
head(res_NBF.BF3)

## log2 fold change (MLE): condition NBF vs BF3hr 
## Wald test p-value: condition NBF vs BF3hr 
## DataFrame with 6 rows and 6 columns
##                        baseMean log2FoldChange     lfcSE       stat      pvalue
##                       <numeric>      <numeric> <numeric>  <numeric>   <numeric>
## Locus10002v1rpkm9.17   155.2077      -3.042106 0.2926413 -10.395343 2.60348e-25
## Locus10003v1rpkm9.17    37.9684       0.556711 0.4656658   1.195517 2.31885e-01
## Locus1000v1rpkm136.27 2187.7091       0.569815 0.0922209   6.178813 6.45854e-10
## Locus1000v2rpkm50.24   590.7223      -0.121061 0.1536421  -0.787939 4.30732e-01
## Locus10016v1rpkm9.16    62.6007       0.277192 0.2866525   0.966997 3.33546e-01
## Locus10018v1rpkm9.16    67.4962      -0.780650 0.3325833  -2.347230 1.89136e-02
##                              padj
##                         <numeric>
## Locus10002v1rpkm9.17  1.96121e-23
## Locus10003v1rpkm9.17  3.59586e-01
## Locus1000v1rpkm136.27 1.04455e-08
## Locus1000v2rpkm50.24  5.66604e-01
## Locus10016v1rpkm9.16  4.71060e-01
## Locus10018v1rpkm9.16  4.99676e-02

sum(res_NBF.BF3$padj < 0.05, na.rm = TRUE)

## [1] 3108

NBF vs BF24hr

res_NBF.BF24 <- results(dds, contrast = c("condition", "NBF", "BF24hr"))
head(res_NBF.BF24)

## log2 fold change (MLE): condition NBF vs BF24hr 
## Wald test p-value: condition NBF vs BF24hr 
## DataFrame with 6 rows and 6 columns
##                        baseMean log2FoldChange     lfcSE      stat      pvalue
##                       <numeric>      <numeric> <numeric> <numeric>   <numeric>
## Locus10002v1rpkm9.17   155.2077      -2.138095 0.2972055 -7.193996 6.29221e-13
## Locus10003v1rpkm9.17    37.9684      -0.116472 0.4538049 -0.256657 7.97443e-01
## Locus1000v1rpkm136.27 2187.7091      -0.173208 0.0910788 -1.901733 5.72060e-02
## Locus1000v2rpkm50.24   590.7223       0.208926 0.1547960  1.349685 1.77117e-01
## Locus10016v1rpkm9.16    62.6007       0.315576 0.2876792  1.096972 2.72654e-01
## Locus10018v1rpkm9.16    67.4962      -1.539807 0.3243016 -4.748069 2.05368e-06
##                              padj
##                         <numeric>
## Locus10002v1rpkm9.17  1.81920e-11
## Locus10003v1rpkm9.17  8.74004e-01
## Locus1000v1rpkm136.27 1.37104e-01
## Locus1000v2rpkm50.24  3.15877e-01
## Locus10016v1rpkm9.16  4.29200e-01
## Locus10018v1rpkm9.16  2.15448e-05

sum(res_NBF.BF24$padj < 0.05, na.rm = TRUE)

## [1] 2613

Volcano Plot

EnhancedVolcano(
  toptable = res_NBF.BF24,
  title = "NBF vs Bloodfed 24hr",
  x = "log2FoldChange",
  y = "padj",
  lab = rownames(res_NBF.BF24),
  labSize = 4,
  pCutoff = 10e-10,
  subtitle = NULL,
  caption = NULL
)

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the EnhancedVolcano package.
##   Please report the issue to the authors.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## ℹ The deprecated feature was likely used in the EnhancedVolcano package.
##   Please report the issue to the authors.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Top Gene Plot

top_df <- data.frame(res_NBF.BF24) |> arrange(padj)
top1 <- rownames(top_df)[1]
top1

## [1] "Locus569v2rpkm130.88"

focal_gene_counts <- plotCounts(
  dds,
  gene = top1,
  intgroup = "condition",
  returnData = TRUE
)

focal_gene_counts |>
  ggplot(aes(x = condition, y = count, fill = condition)) +
  geom_boxplot(alpha = 0.5, outlier.shape = NA) +
  geom_point(size = 4, shape = 21) +
  theme_bw(base_size = 20) +
  labs(
    x = "Feeding Treatment",
    y = "Count"
  ) +
  theme(legend.position = "none")

BF3hr vs BF24hr

res_BF3.BF24 <- results(dds, contrast = c("condition", "BF3hr", "BF24hr"))
head(res_BF3.BF24)

## log2 fold change (MLE): condition BF3hr vs BF24hr 
## Wald test p-value: condition BF3hr vs BF24hr 
## DataFrame with 6 rows and 6 columns
##                        baseMean log2FoldChange     lfcSE      stat      pvalue
##                       <numeric>      <numeric> <numeric> <numeric>   <numeric>
## Locus10002v1rpkm9.17   155.2077      0.9040112  0.278473   3.24632 1.16907e-03
## Locus10003v1rpkm9.17    37.9684     -0.6731836  0.473807  -1.42080 1.55376e-01
## Locus1000v1rpkm136.27 2187.7091     -0.7430229  0.092877  -8.00007 1.24345e-15
## Locus1000v2rpkm50.24   590.7223      0.3299864  0.156404   2.10983 3.48729e-02
## Locus10016v1rpkm9.16    62.6007      0.0383842  0.300676   0.12766 8.98418e-01
## Locus10018v1rpkm9.16    67.4962     -0.7591571  0.325100  -2.33515 1.95357e-02
##                              padj
##                         <numeric>
## Locus10002v1rpkm9.17  4.04519e-03
## Locus10003v1rpkm9.17  2.56957e-01
## Locus1000v1rpkm136.27 2.67978e-14
## Locus1000v2rpkm50.24  7.56715e-02
## Locus10016v1rpkm9.16  9.31546e-01
## Locus10018v1rpkm9.16  4.70265e-02

sum(res_BF3.BF24$padj < 0.05, na.rm = TRUE)

## [1] 3456

top_df_BF <- data.frame(res_BF3.BF24) |> arrange(padj)
top1_BF <- rownames(top_df_BF)[1]
top1_BF

## [1] "Locus116v3rpkm256.78"

Heatmaps

meta <- data.frame(colData(dds))
norm_mat <- assay(dds_vst1)

top25_DE <- rownames(top_df)[1:25]
norm_mat_sel <- norm_mat[match(top25_DE, rownames(norm_mat)), ]

meta_sort <- meta |> arrange(condition) |> select(condition)

# Match column order to metadata
norm_mat_sel <- norm_mat_sel[, match(rownames(meta_sort), colnames(norm_mat_sel))]

pheatmap(
  norm_mat_sel,
  annotation_col = meta_sort,
  cluster_cols = FALSE,
  show_rownames = FALSE,
  scale = "row"
)

Digestion Gene Analysis

res <- results(dds)

res_df <- as.data.frame(res)
res_df$id <- rownames(res_df)
res_df <- res_df[!is.na(res_df$padj), ]

annotated_res <- merge(res_df, raw, by = "id")

digestion_genes <- annotated_res |>
  subset(
    padj < 0.05 &
    grepl("trypsin|protease|peptidase|chymotrypsin|carboxypeptidase|aminopeptidase",
          pred_func,
          ignore.case = TRUE)
  )

nrow(digestion_genes)

## [1] 85

sum(digestion_genes$log2FoldChange > 0)

## [1] 23

sum(digestion_genes$log2FoldChange < 0)

## [1] 62

n_total_DE <- sum(res$padj < 0.05, na.rm = TRUE)
n_digestion_DE <- nrow(digestion_genes)
percentage <- (n_digestion_DE / n_total_DE) * 100
percentage

## [1] 3.252966

Digestion Summary Plot

digestion_summary <- data.frame(
  Direction = c("Upregulated", "Downregulated"),
  Count = c(
    sum(digestion_genes$log2FoldChange > 0),
    sum(digestion_genes$log2FoldChange < 0)
  )
)

ggplot(digestion_summary, aes(x = Direction, y = Count, fill = Direction)) +
  geom_bar(stat = "identity") +
  theme_bw(base_size = 18)

Volcano Highlighting Digestion Genes

res_df$digestion <- ifelse(res_df$id %in% digestion_genes$id,
                           "Digestion gene",
                           "Other gene")

res_df$label <- ifelse(res_df$padj < 0.05 &
                       abs(res_df$log2FoldChange) > 1,
                       res_df$id,
                       "")

ggplot(res_df,
       aes(x = log2FoldChange,
           y = -log10(padj))) +

  geom_point(aes(color = digestion),
             alpha = 0.7,
             size = 2) +

  geom_text_repel(aes(label = label),
                  size = 4,
                  max.overlaps = 20) +

  geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed") +

  scale_color_manual(values = c("Digestion gene" = "red",
                                "Other gene" = "gray50")) +

  theme_bw(base_size = 16) +
  labs(
    title = "Digestion Realted DEGs",
    x = "Log2 Fold Change",
    y = "-Log10 Adjusted P-value",
    color = ""
  ) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))