Lets see your code. 8000 significant genesets indicate a flaw somewhere.

Hi Gordon,

Thanks for the advice, I gave camera a try, and it yielded around ~1,000 statistically significant gene sets from the same 11,000 gene sets tested. This is still much larger than what I got with preranked GSEA, so it sounds like there might be an error in my code?

The gene sets I am using are the GO Biological Process terms, the KEGG gene sets, Reactome gene sets, and two custom lists. fry/camera are run for each of these groups of gene sets separately, and then FDR-controlled separately as well. Maybe the hierarchical nature of the GO terms is inflating the percentage of significant gene sets?

Regarding my data, from DESeq2 results I see that there is on the order of 1,000 genes up/down with a log2FC cutoff of 0.585 (equivalent to a FC cutoff of 1.5). I also have another contrast that yields much fewer DEGs (on the order of 100 up/down for the same 0.585 cutoff), but that contrast is also yielding about 1,000 significant gene sets, same as the first contrast.

Here is an edited version of my script with the general workflow:

# Create sample metadata data frame
# 4 replicates in each condition (16 samples total)
conditions <- c("Sample_A_Treatment_1", "Sample_A_Treatment_2", "Sample_B_Treatment_1", "Sample_B_Treatment_2")
run_list <- paste0("DW", sprintf("%02d", 1:16))    # sample IDs
repeated_conditions <- rep(conditions, each = 4)
sample_id_list <- rep(run_list, length.out = length(run_list) * length(conditions))
samples <- data.frame(
  run = run_list,
  condition = factor(rep(repeated_conditions, length.out = length(run_list))),
  sample_id = run_list)

# Import counts
txi <- tximport(files, type = "salmon", tx2gene = tx2gene, txOut = FALSE)

# Create DGEList object
y <- DGEList(counts = txi$counts, group = samples$condition)

# Filter lowly expressed genes
design_matrix <- model.matrix(~ 0 + group, data = y$samples)
colnames(design_matrix) <- levels(y$samples$group)
keep <- filterByExpr(y, design = design_matrix)
y <- y[keep, , keep.lib.sizes=FALSE]

# Normalize for library size, estimate dispersions, fit QL model
y <- calcNormFactors(y)
y <- estimateDisp(y, design_matrix, robust=TRUE)
fit <- glmQLFit(y, design_matrix, robust=TRUE)

# Define contrasts
contrast_1 = "Sample_A_Treatment_2-Sample_A_Treatment_1"
contrast_2 = "Sample_B_Treatment_2-Sample_B_Treatment_1"
all_contrasts <- c(contrast_1, contrast_2)
contrast_matrices <- makeContrasts(contrasts = all_contrasts, levels = design_matrix)

# Load gene sets from MSigDB and custom GMT files
all_gene_sets <- list()
# code to load gene sets (skipped for brevity)
indexed_gene_sets <- lapply(all_gene_sets, function(db) {
  ids2indices(db, identifiers = rownames(y))
})

# run fry
fry_results_list <- list() # Initialize an empty list to store all results

# Loop through each defined contrast
for (i in 1:length(all_contrasts)) {

  current_contrast_name <- all_contrasts[i]
  cat(paste("\n--- Processing Contrast: ", current_contrast_name, " ---\n"))

  # Get the specific contrast matrix for the current contrast
  current_contrast_matrix <- contrast_matrices[, i, drop = FALSE] # drop=FALSE keeps it a matrix

  # Loop through each gene set database (e.g., GO:BP, KEGG, Custom)
  for (db_name in names(indexed_gene_sets)) {

    cat(paste("  Running fry for Database:", db_name, "...\n"))

    # Run the fry test for the current contrast and database
    fry_res <- fry(y,
                   index = indexed_gene_sets[[db_name]],
                   design = design_matrix,
                   contrast = current_contrast_matrix)

    # Add the database name, pathway names, and the contrast name for clarity
    fry_res <- fry_res %>%
      rownames_to_column(var = "Pathway") %>%
      mutate(
        Database = db_name,
        Contrast = current_contrast_name # Add the contrast name here
      )

    # Append the results to main list
    fry_results_list[[paste(db_name, current_contrast_name, sep="_")]] <- fry_res
  }
}

# Combine all results into one table, calculate FDR, and sort by it.
final_fry_results <- bind_rows(fry_results_list) %>%
  # Group by Contrast and Database before adjusting p-values.
  group_by(Contrast, Database) %>%
  mutate(FDR = p.adjust(PValue, method = "BH")) %>%
  ungroup() %>% # Ungroup for easier downstream manipulation
  select(Contrast, Database, Pathway, NGenes, Direction, PValue, FDR) %>%
  arrange(Contrast, PValue) # Sort by contrast and then by p-value

# save fry results
write.csv(final_fry_results, output_file, row.names = FALSE)

# run camera using same loop as for fry