Differential Expression Analysis with additive model

Setup and Package Loading

# Set root directory
here::i_am("scripts/additive_DE.Rmd")

# Load necessary packages
suppressPackageStartupMessages({
  library(tidyverse)
  library(DESeq2)
  library(apeglm)
  library(ashr)
  library(grid)
  library(VennDiagram)
  library(here)
})

Data Import and Preprocessing

# Load abundance matrix and sample information 

abundance_matrix <- read.table(here::here("data", "DE", "abundance_matrix.txt"), header = TRUE, row.names = 1)
head(abundance_matrix)

##                         X39_Melops1040 X38_Melops1037 X37_Melops1036
## TRINITY_GG_716_c838_g1               0              2              0
## TRINITY_GG_1114_c422_g3              0              0              0
## TRINITY_GG_1463_c526_g6              0              7              0
## TRINITY_GG_1586_c745_g1              0              0              0
## TRINITY_GG_560_c228_g1               1              2              1
## TRINITY_GG_606_c418_g2               1              1              1
##                         X36_Melops1035 X35_Melops1033 X34_Melops1025
## TRINITY_GG_716_c838_g1               0              0              0
## TRINITY_GG_1114_c422_g3              0              0              0
## TRINITY_GG_1463_c526_g6              3              0              0
## TRINITY_GG_1586_c745_g1              0              0              0
## TRINITY_GG_560_c228_g1               0              0              0
## TRINITY_GG_606_c418_g2               1              3              6
##                         X33_Melops1024 X32_Melops1023 X31_Melops1022
## TRINITY_GG_716_c838_g1               0              0              0
## TRINITY_GG_1114_c422_g3              0              0              0
## TRINITY_GG_1463_c526_g6              1              0              2
## TRINITY_GG_1586_c745_g1              0              0              0
## TRINITY_GG_560_c228_g1               0              0              0
## TRINITY_GG_606_c418_g2               5              0              2
##                         X30_Melops1021 X29_Melops1211 X28_Melops1210
## TRINITY_GG_716_c838_g1               0              0              0
## TRINITY_GG_1114_c422_g3              0              0              0
## TRINITY_GG_1463_c526_g6              1              2              0
## TRINITY_GG_1586_c745_g1              0              0              3
## TRINITY_GG_560_c228_g1               0              4              0
## TRINITY_GG_606_c418_g2               2              0              1
##                         X27_Melops1209 X26_Melops1206 X25_Melops1205
## TRINITY_GG_716_c838_g1               0              0              0
## TRINITY_GG_1114_c422_g3              0              0              0
## TRINITY_GG_1463_c526_g6              0              0              0
## TRINITY_GG_1586_c745_g1              0              0              0
## TRINITY_GG_560_c228_g1               0              0              0
## TRINITY_GG_606_c418_g2               0              3              5
##                         X24_Melops1082 X23_Melops1065 X22_Melops1061
## TRINITY_GG_716_c838_g1               1              0              0
## TRINITY_GG_1114_c422_g3              0              0              0
## TRINITY_GG_1463_c526_g6              0              0              0
## TRINITY_GG_1586_c745_g1              0              0              0
## TRINITY_GG_560_c228_g1               0              0              0
## TRINITY_GG_606_c418_g2               0              4              2
##                         X21_Melops1058 X20_Melops1046 X19_Melops1131
## TRINITY_GG_716_c838_g1               0              1              0
## TRINITY_GG_1114_c422_g3              0              0              0
## TRINITY_GG_1463_c526_g6              0              0              0
## TRINITY_GG_1586_c745_g1              0              0              0
## TRINITY_GG_560_c228_g1               0              0              1
## TRINITY_GG_606_c418_g2               1              1              2
##                         X18_Melops1075 X17_Melops1072 X16_Melops1017
## TRINITY_GG_716_c838_g1               0              0              0
## TRINITY_GG_1114_c422_g3              0              0              0
## TRINITY_GG_1463_c526_g6              0              0              0
## TRINITY_GG_1586_c745_g1              0              0              0
## TRINITY_GG_560_c228_g1               1              0              1
## TRINITY_GG_606_c418_g2               1              0              0
##                         X15_Melops1014 X14_Melops1202 X13_Melops1200
## TRINITY_GG_716_c838_g1               0              0              0
## TRINITY_GG_1114_c422_g3              0              0              0
## TRINITY_GG_1463_c526_g6              0              0              0
## TRINITY_GG_1586_c745_g1              0              0              0
## TRINITY_GG_560_c228_g1               2              0              0
## TRINITY_GG_606_c418_g2               2              3              4
##                         X12_Melops1197 X11_Melops1194 X10_Melops1192
## TRINITY_GG_716_c838_g1               0              0              0
## TRINITY_GG_1114_c422_g3              1              0              0
## TRINITY_GG_1463_c526_g6              2              1              0
## TRINITY_GG_1586_c745_g1              0              0              0
## TRINITY_GG_560_c228_g1               2              0              1
## TRINITY_GG_606_c418_g2               0              3              5
##                         X09_Melops1158 X08_Melops1126 X07_Melops1015
## TRINITY_GG_716_c838_g1               0              1              0
## TRINITY_GG_1114_c422_g3              0              0              0
## TRINITY_GG_1463_c526_g6              2              1              1
## TRINITY_GG_1586_c745_g1              3              0              0
## TRINITY_GG_560_c228_g1               4              0              0
## TRINITY_GG_606_c418_g2               5              2              0
##                         X06_Melops1010 X05_Melops1259 X04_Melops1239
## TRINITY_GG_716_c838_g1               0              0              0
## TRINITY_GG_1114_c422_g3              0              0              0
## TRINITY_GG_1463_c526_g6              0              0              0
## TRINITY_GG_1586_c745_g1              0              0              0
## TRINITY_GG_560_c228_g1               0              0              0
## TRINITY_GG_606_c418_g2               4              0              4
##                         X03_Melops1218 X02_Melops1195 X01_Melops1188
## TRINITY_GG_716_c838_g1               0              0              0
## TRINITY_GG_1114_c422_g3              0              0              0
## TRINITY_GG_1463_c526_g6              0              0              0
## TRINITY_GG_1586_c745_g1              0              0              0
## TRINITY_GG_560_c228_g1               0              0              0
## TRINITY_GG_606_c418_g2               3              0              2

sample_info <- read.table(
  here::here("data", "DE", "samples.txt"),
  header = FALSE,
  col.names = c("condition", "sample")
)
sample_info

##    condition        sample
## 1        W18 39_Melops1040
## 2        W18 38_Melops1037
## 3        W18 37_Melops1036
## 4        W18 36_Melops1035
## 5        W18 35_Melops1033
## 6        W15 34_Melops1025
## 7        W15 33_Melops1024
## 8        W15 32_Melops1023
## 9        W15 31_Melops1022
## 10       W15 30_Melops1021
## 11       W12 29_Melops1211
## 12       W12 28_Melops1210
## 13       W12 27_Melops1209
## 14       W12 26_Melops1206
## 15       W12 25_Melops1205
## 16       H18 24_Melops1082
## 17       H18 23_Melops1065
## 18       H18 22_Melops1061
## 19       H18 21_Melops1058
## 20       H18 20_Melops1046
## 21       H15 19_Melops1131
## 22       H15 18_Melops1075
## 23       H15 17_Melops1072
## 24       H15 16_Melops1017
## 25       H15 15_Melops1014
## 26       H12 14_Melops1202
## 27       H12 13_Melops1200
## 28       H12 12_Melops1197
## 29       H12 11_Melops1194
## 30       H12 10_Melops1192
## 31       S18 09_Melops1158
## 32       S18 08_Melops1126
## 33       S15 07_Melops1015
## 34       S15 06_Melops1010
## 35       S12 05_Melops1259
## 36       S12 04_Melops1239
## 37       S12 03_Melops1218
## 38       S12 02_Melops1195
## 39       S12 01_Melops1188

Prepare dataset and run DESeq2

# Parse temperature and origin from condition string. Assumes condition is formatted like "W12", "S15", etc.
sample_info <- sample_info %>%
  mutate(
    origin = substr(condition, 1, 1),
    temperature = substr(condition, 2, nchar(condition))
  ) %>%
  mutate(
    origin = recode(origin, "W" = "West", "H" = "Hybrid", "S" = "South"),
    temperature = factor(temperature, levels = c("12", "15", "18")),
    origin = factor(origin)
  )


# Ensure column names of count matrix match sample names
colnames(abundance_matrix) <- sample_info$sample

# Round counts and ensure non-negative integers (required by DESeq2)
abundance_matrix <- round(pmax(abundance_matrix, 0))  

# Create DESeq2 dataset with additive model (temperature + origin) 
dds <- DESeqDataSetFromMatrix(
  countData = abundance_matrix,
  colData = sample_info,
  design = ~ temperature + origin
)

# Filter out lowly expressed transcripts (keep those with total count ≥ 10 across all sampl
dds <- dds[rowSums(counts(dds)) >= 10, ]
# Run DESeq2 pipeline (estimate size factors, dispersions, and fit model)
dds <- DESeq(dds)
# Extract coefficient names for downstream contrasts
coef_names <- resultsNames(dds)

# Plot dispersion estimates across transcripts to assess fit to the negative binomial model used by DESeq2
plotDispEsts(dds)

Define contrasts

# Temperature contrasts (across all origins)
temperature_levels <- levels(sample_info$temperature)
temperature_contrasts <- t(combn(temperature_levels, 2))
colnames(temperature_contrasts) <- c("condition1", "condition2")
temperature_contrasts <- as.data.frame(temperature_contrasts)
temperature_contrasts

##   condition1 condition2
## 1         12         15
## 2         12         18
## 3         15         18

# Origin contrasts (across all temperatures)
origin_levels <- levels(sample_info$origin)
origin_contrasts <- t(combn(origin_levels, 2))
colnames(origin_contrasts) <- c("condition1", "condition2")
origin_contrasts <- as.data.frame(origin_contrasts)
origin_contrasts

##   condition1 condition2
## 1     Hybrid      South
## 2     Hybrid       West
## 3      South       West

Helper function to check if a contrast corresponds to a direct coefficient in the DESeq2 model.

# Determines whether the comparison involves the reference level provided by DESeq2 as a coefficient. This determines whether we can use 'coef' (apeglm) or need to use 'contrast' (normal).
is_direct_coef <- function(factor_name, level1, level2, coef_names) {
  ref <- levels(sample_info[[factor_name]])[1]
  if (level2 == ref) {
    coef_pattern <- paste0(factor_name, "_", level1, "_vs_", ref)
    return(coef_pattern %in% coef_names)
  } else if (level1 == ref) {
    coef_pattern <- paste0(factor_name, "_", level2, "_vs_", ref)
    return(coef_pattern %in% coef_names)
  }
  return(FALSE)
}

Main function to run pairwise differential expression analysis for a given factor (e.g., temperature or origin).

# Automatically selects shrinkage method based on whether the comparison involves the reference level. Saves full and significant results, generates MA plots, and returns a summary table
run_contrast_analysis <- function(contrast_table, factor_name, output_dir, alpha = 0.05, lfcThreshold = 0.0, coef_names) {
  results_list <- list()
  counts_summary <- data.frame(
    Comparison = character(),
    Shrinkage = character(),
    Upregulated = integer(),
    Downregulated = integer(),
    stringsAsFactors = FALSE
  )
  # Create output directory if it doesn't exist
  dir.create(here::here(output_dir), recursive = TRUE, showWarnings = FALSE)
  
  for (i in 1:nrow(contrast_table)) {
    level1 <- contrast_table$condition1[i]
    level2 <- contrast_table$condition2[i]
    name <- paste(level2, "vs", level1, sep = "_")
    
    
    # Check if the comparison involves the reference level provided by DESeq2 as a coefficient
    # Use 'coef' with apeglm shrinkage for reference-level comparisons; use 'contrast' with normal shrinkage     otherwise
    direct_coef <- is_direct_coef(factor_name, level1, level2, coef_names)
    
    if (direct_coef) {
      coef_name <- paste0(factor_name, "_", level2, "_vs_", levels(sample_info[[factor_name]])[1])
      res <- lfcShrink(dds, coef = coef_name, type = "apeglm")
      shrinkage <- "apeglm"
    } else {
      res <- lfcShrink(dds, contrast = c(factor_name, level2, level1), type = "normal")
      shrinkage <- "normal"
    }
    
    results_list[[name]] <- res
    
    # Extract significant transcripts based on adjusted p-value and log2 fold change threshold
    sig <- res[which(res$padj < alpha & abs(res$log2FoldChange) > lfcThreshold), ]
    up <- sum(sig$log2FoldChange > 0)
    down <- sum(sig$log2FoldChange < 0)
    
    # Save full and significant results to CSV
    write.csv(as.data.frame(res), file = here::here(output_dir, paste0(name, "_DE_results.csv")), row.names = TRUE)
    write.csv(as.data.frame(sig), file = here::here(output_dir, paste0(name, "_significant_transcripts.csv")), row.names = TRUE)
    
    # Generate and save MA plot with up/downregulated transcript counts
    png(filename = here::here(output_dir, paste0(name, "_MA_plot.png")), width = 800, height = 600)
    plotMA(res, ylim = c(-10, 10), main = paste("MA Plot:", name, "-", shrinkage), alpha = alpha, colSig = "red")
    text(x = max(res$baseMean, na.rm = TRUE) * 0.7, y = 9, labels = paste("↑:", up), col = "blue", cex = 0.8)
    text(x = max(res$baseMean, na.rm = TRUE) * 0.7, y = -9, labels = paste("↓:", down), col = "blue", cex = 0.8)
    dev.off()
    
    # Also display MA plot inline in the report
    plotMA(res, ylim = c(-10, 10), main = paste("MA Plot:", name, "-", shrinkage), alpha = alpha, colSig = "red")
    text(x = max(res$baseMean, na.rm = TRUE) * 0.7, y = 9, labels = paste("↑:", up), col = "blue", cex = 0.8)
    text(x = max(res$baseMean, na.rm = TRUE) * 0.7, y = -9, labels = paste("↓:", down), col = "blue", cex = 0.8)
    
    # Add summary row to the results table
    counts_summary <- rbind(counts_summary, data.frame(
      Comparison = name,
      Shrinkage = shrinkage,
      Upregulated = up,
      Downregulated = down
    ))
  }
  # Return list of DESeq2 results and summary tab
  return(list(results = results_list, summary = counts_summary))
}

Temperature contrast

# Temperature contrasts
temperature_results <- run_contrast_analysis(
  contrast_table = temperature_contrasts,
  factor_name = "temperature",
  output_dir = "results/DE/temperature_contrast_additive_model",
  coef_names = coef_names
)

print(temperature_results$summary)

##   Comparison Shrinkage Upregulated Downregulated
## 1   15_vs_12    apeglm       10061          7824
## 2   18_vs_12    apeglm       12040          9429
## 3   18_vs_15    normal        1385          1551

## Save summary table
write.csv(
  temperature_results$summary,
  here::here("results", "DE", "temperature_contrast_additive_model", "contrast_summary.csv"),
  row.names = FALSE
)

Origin contrast

# Origin contrast
origin_results <- run_contrast_analysis(
  contrast_table = origin_contrasts,
  factor_name = "origin",
  output_dir = "results/DE/origin_contrast_additive_model",
  coef_names = coef_names
)

print(origin_results$summary)

##        Comparison Shrinkage Upregulated Downregulated
## 1 South_vs_Hybrid    apeglm          11            19
## 2  West_vs_Hybrid    apeglm          24            26
## 3   West_vs_South    normal         574           561

## Save summary table
write.csv(
  origin_results$summary,
  here::here("results", "DE", "origin_contrast_additive_model", "contrast_summary.csv"),
  row.names = FALSE
)

Venn diagrams

Define venn plotting function

# Plot Venn diagram to visualize overlap between sets of significant transcripts
plot_venn <- function(transcript_lists, contrast_name, output_dir) {
  venn_plot <- venn.diagram(
    x = transcript_lists,
    category.names = names(transcript_lists),
    filename = NULL,
    fill = c("red", "green", "blue"),
    alpha = 0.5,
    cex = 2,
    fontface = "bold"
  )
  
  png(filename = here::here(output_dir, paste0("venn_", contrast_name, ".png")), width = 800, height = 800)
  grid.newpage()
  grid.draw(venn_plot)
  dev.off()
  
  grid.newpage()
  grid.draw(venn_plot)
}

Prepare DE transcripts sets and plot venn diagrams for temperature contrast

# Prepare DE sets for temperature contrasts
temp_contrasts <- c("15_vs_12", "18_vs_15", "18_vs_12")
temp_sets <- setNames(lapply(temp_contrasts, function(name) {
  file <- here::here("results", "DE", "temperature_contrast_additive_model", paste0(name, "_significant_transcripts.csv"))
  rownames(read.csv(file, row.names = 1))
}), c("15 vs 12", "18 vs 15", "18 vs 12"))

# Plot Venn diagram
plot_venn(temp_sets, "temperature_contrasts", "results/DE/temperature_contrast_additive_model")

Prepare DE transcripts sets and plot venn diagrams for origin contrast

# Prepare DE sets for origin contrasts
origin_contrasts <- c("South_vs_Hybrid", "West_vs_Hybrid", "West_vs_South")
origin_labels <- c("South vs Hybrid", "West vs Hybrid", "West vs South")

origin_sets <- setNames(lapply(seq_along(origin_contrasts), function(i) {
  file <- here::here("results", "DE", "origin_contrast_additive_model", paste0(origin_contrasts[i], "_significant_transcripts.csv"))
  rownames(read.csv(file, row.names = 1))
}), origin_labels)

# Plot Venn diagram
plot_venn(origin_sets, "origin_contrasts", "results/DE/origin_contrast_additive_model")

Function to load transcript sets

# Load significant transcript sets used to visualize shared and unique DE transcripts across conditions
load_transcript_sets <- function(contrast_names, contrast_dir) {
  sets <- list()
  for (name in contrast_names) {
    file <- here::here("results", "DE", contrast_dir, paste0(name, "_significant_transcripts.csv"))
    if (file.exists(file)) {
      df <- read.csv(file, row.names = 1)
      sets[[name]] <- rownames(df)
    } else {
      warning(paste("File not found:", file))
      sets[[name]] <- character(0)
    }
  }
  return(sets)
}

Function to summarize Venn regions

# Summarize the number of transcripts in each Venn region (shared or unique across contrasts)
venn_summary_table <- function(sets) {
  all_genes <- unique(unlist(sets))
  summary <- sapply(all_genes, function(gene) {
    present_in <- names(sets)[sapply(sets, function(s) gene %in% s)]
    paste(sort(present_in), collapse = " & ")
  })
  as.data.frame(table(summary), stringsAsFactors = FALSE) |>
    setNames(c("Region", "Count")) |>
    dplyr::arrange(desc(Count))
}

Temperature Venn summary

temp_contrasts <- c("15_vs_12", "18_vs_15", "18_vs_12")
temp_sets <- load_transcript_sets(temp_contrasts, "temperature_contrast_additive_model")
temp_venn_summary <- venn_summary_table(temp_sets)
print(temp_venn_summary)

##                           Region Count
## 1            15_vs_12 & 18_vs_12 12549
## 2                       18_vs_12  7064
## 3                       15_vs_12  4182
## 4            18_vs_12 & 18_vs_15  1271
## 5 15_vs_12 & 18_vs_12 & 18_vs_15   585
## 6            15_vs_12 & 18_vs_15   569
## 7                       18_vs_15   511

# Save venn summary
write.csv(
  temp_venn_summary,
  here::here("results", "DE", "temperature_contrast_additive_model", "venn_summary.csv"),
  row.names = FALSE
)

Origin Venn summary

origin_contrasts <- c("South_vs_Hybrid", "West_vs_Hybrid", "West_vs_South")
origin_sets <- load_transcript_sets(origin_contrasts, "origin_contrast_additive_model")
origin_venn_summary <- venn_summary_table(origin_sets)
print(origin_venn_summary)

##                            Region Count
## 1                   West_vs_South  1110
## 2                  West_vs_Hybrid    49
## 3 South_vs_Hybrid & West_vs_South    24
## 4                 South_vs_Hybrid     6
## 5  West_vs_Hybrid & West_vs_South     1

# Save venn summary
write.csv(
  origin_venn_summary,
  here::here("results", "DE", "origin_contrast_additive_model", "venn_summary.csv"),
  row.names = FALSE
)

Transcript Extraction for functional enrichment

Tier 1: Shared temperature-responsive transcripts across origins. Candidates for general temperature response.

Biological Rationale

Tier 1 identifies transcripts that are consistently regulated by temperature across all origins. By intersecting results from all pairwise temperature contrasts, this approach isolates genes involved in a general thermal response, independent of population-specific effects.

# Extract transcripts that are consistently DE across all temperature contrasts. 
tier1_transcripts <- Reduce(intersect, temp_sets)

# Save for functional enrichment
write.csv(
  data.frame(transcript_id = tier1_transcripts),
  here::here("results", "functional_enrichment", "tier1_shared_temperature_response", "tier1_shared_temperature_transcripts.csv"),
  row.names = FALSE
)

Tier 2: Local adaptation candidates (constitutive + plasticity divergence)

Biological Rationale

Tier 2 aims to identify transcripts that may be involved in local adaptation between corkwing wrasse populations. These candidates are defined by the overlap between constitutive differences and interaction effects:

Constitutive differences are identified using the additive model (~ origin + temperature) by comparing origins (West vs South) while controlling for temperature. These represent baseline expression differences that are consistent across temperatures and may reflect fixed genetic divergence.
Interaction effects are identified using a likelihood ratio test (LRT) comparing models with and without the origin:temperature interaction term. These transcripts show origin-specific responses to temperature, indicating plasticity divergence.

By intersecting these two sets, we isolate transcripts that are: - Diverged between origins, and - Differentially responsive to temperature depending on origin

This overlap represents the strongest candidates for local adaptation, as they combine fixed and environmentally responsive divergence.

# Load constitutive origin differences (South vs West)
constitutive_file <- here::here("results", "DE", "origin_contrast_additive_model", "West_vs_South_significant_transcripts.csv")
constitutive_transcripts <- rownames(read.csv(constitutive_file, row.names = 1))

# Load interaction effect transcripts from LRT script
interaction_file <- here::here("results", "DE", "interaction_effects", "significant_interaction_transcripts.csv")
interaction_transcripts <- rownames(read.csv(interaction_file, row.names = 1))

# Overlap = strongest candidates for local adaptation
tier2_transcripts <- intersect(constitutive_transcripts, interaction_transcripts)

# Save for functional enrichment
write.csv(
  data.frame(transcript_id = tier2_transcripts),
  here::here("results", "functional_enrichment", "tier2_local_adaptation", "tier2_local_adaptation_candidates.csv"),
  row.names = FALSE
)

Tier 3: Hybrid misregulation Tier 3 is analyzed separately using the condition model, where each origin-temperature combination is treated as a distinct group. This approach identifies transcripts where hybrids differ from both parental origins at a given temperature. Misregulated transcripts can be further classified as transgressive or intermediate based on expression direction.

Summary of Transcript Extraction

data.frame(
  Tier = c("Tier 1", "Tier 2"),
  Transcripts = c(length(tier1_transcripts), length(tier2_transcripts))
)

##     Tier Transcripts
## 1 Tier 1         585
## 2 Tier 2          32

Summary

This script performs Tier 1 and Tier 2 differential expression analyses using an additive DESeq2 model to investigate temperature responses and local adaptation in Symphodus melops. Pairwise temperature contrasts revealed 585 transcripts consistently regulated by temperature across all origins (Tier 1). For Tier 2, constitutive origin differences were defined using West vs South origin fish, and their overlap with transcripts showing origin-dependent temperature responses yielded 32 strong candidates for local adaptation. These transcript sets were extracted and saved for downstream functional enrichment.

Session Info

sessionInfo()

## R version 4.5.1 (2025-06-13)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0  LAPACK version 3.10.0
## 
## locale:
##  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
##  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
##  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
## [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
## 
## time zone: Europe/Oslo
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] here_1.0.2                  VennDiagram_1.7.3          
##  [3] futile.logger_1.4.3         ashr_2.2-63                
##  [5] apeglm_1.30.0               DESeq2_1.48.2              
##  [7] SummarizedExperiment_1.38.1 Biobase_2.68.0             
##  [9] MatrixGenerics_1.20.0       matrixStats_1.5.0          
## [11] GenomicRanges_1.60.0        GenomeInfoDb_1.44.3        
## [13] IRanges_2.42.0              S4Vectors_0.46.0           
## [15] BiocGenerics_0.54.0         generics_0.1.4             
## [17] lubridate_1.9.4             forcats_1.0.1              
## [19] stringr_1.5.2               dplyr_1.1.4                
## [21] purrr_1.1.0                 readr_2.1.5                
## [23] tidyr_1.3.1                 tibble_3.3.0               
## [25] ggplot2_4.0.0               tidyverse_2.0.0            
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.1        farver_2.1.2            S7_0.2.0               
##  [4] fastmap_1.2.0           digest_0.6.37           timechange_0.3.0       
##  [7] lifecycle_1.0.4         invgamma_1.2            magrittr_2.0.4         
## [10] compiler_4.5.1          rlang_1.1.6             sass_0.4.10            
## [13] tools_4.5.1             yaml_2.3.10             lambda.r_1.2.4         
## [16] knitr_1.50              S4Arrays_1.8.1          DelayedArray_0.34.1    
## [19] plyr_1.8.9              RColorBrewer_1.1-3      abind_1.4-8            
## [22] BiocParallel_1.42.2     withr_3.0.2             numDeriv_2016.8-1.1    
## [25] scales_1.4.0            MASS_7.3-55             bbmle_1.0.25.1         
## [28] cli_3.6.5               mvtnorm_1.3-3           rmarkdown_2.30         
## [31] crayon_1.5.3            httr_1.4.7              tzdb_0.5.0             
## [34] bdsmatrix_1.3-7         cachem_1.1.0            parallel_4.5.1         
## [37] formatR_1.14            XVector_0.48.0          vctrs_0.6.5            
## [40] Matrix_1.4-0            jsonlite_2.0.0          hms_1.1.3              
## [43] mixsqp_0.3-54           irlba_2.3.5.1           locfit_1.5-9.12        
## [46] jquerylib_0.1.4         glue_1.8.0              emdbook_1.3.14         
## [49] codetools_0.2-18        stringi_1.8.7           gtable_0.3.6           
## [52] UCSC.utils_1.4.0        pillar_1.11.1           htmltools_0.5.8.1      
## [55] GenomeInfoDbData_1.2.14 truncnorm_1.0-9         R6_2.6.1               
## [58] rprojroot_2.1.1         evaluate_1.0.5          lattice_0.20-45        
## [61] futile.options_1.0.1    SQUAREM_2021.1          bslib_0.9.0            
## [64] Rcpp_1.1.0              coda_0.19-4.1           SparseArray_1.8.1      
## [67] xfun_0.53               pkgconfig_2.0.3

Differential Expression Analysis of Corwing Wrasse

Carlota Myhre de Gouveia

2025-10-03

Differential Expression Analysis with additive model

Setup and Package Loading

Data Import and Preprocessing

Prepare dataset and run DESeq2

Define contrasts

Helper function to check if a contrast corresponds to a direct coefficient in the DESeq2 model.

Main function to run pairwise differential expression analysis for a given factor (e.g., temperature or origin).

Temperature contrast

Origin contrast

Venn diagrams

Define venn plotting function

Prepare DE transcripts sets and plot venn diagrams for temperature contrast

Prepare DE transcripts sets and plot venn diagrams for origin contrast

Function to load transcript sets

Function to summarize Venn regions

Temperature Venn summary

Origin Venn summary

Transcript Extraction for functional enrichment

Tier 1: Shared temperature-responsive transcripts across origins. Candidates for general temperature response.

Biological Rationale

Tier 2: Local adaptation candidates (constitutive + plasticity divergence)

Biological Rationale

Summary of Transcript Extraction

Summary

Session Info