R Notebook

title: “KDM_fishers_exact + correction” output: html_document date: “2025-09-22”

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

df <- read.csv("kdm_data.csv", stringsAsFactors = FALSE)
fisher_vs_lexa <- function(df, genotype_label) {
  d <- df %>% filter(Genotype == genotype_label)
  stopifnot("LexA" %in% d$RNAi)
  lex <- d %>% filter(RNAi == "LexA") %>% slice(1)
  out_list <- lapply(setdiff(unique(d$RNAi), "LexA"), function(mod) {
    mrow <- d %>% filter(RNAi == mod) %>% slice(1)
    mat <- matrix(
      c(mrow$Males, mrow$Females, lex$Males, lex$Females),
      nrow = 2, byrow = TRUE,
      dimnames = list(RNAi = c(mod, "LexA"), Sex = c("M","F"))
    )
    ft_two  <- fisher.test(mat, alternative = "two.sided")
    ft_less <- fisher.test(mat, alternative = "less")  # modifier lower than LexA
    data.frame(
      Genotype      = genotype_label,
      RNAi          = mod,
      M_mod         = mrow$Males,
      F_mod         = mrow$Females,
      M_LexA        = lex$Males,
      F_LexA        = lex$Females,
      prop_mod      = mrow$Males / (mrow$Males + mrow$Females),
      prop_LexA     = lex$Males / (lex$Males + lex$Females),
      OR_hat        = unname(ft_two$estimate),
      OR_LCL95      = ft_two$conf.int[1],
      OR_UCL95      = ft_two$conf.int[2],
      p_two_sided   = ft_two$p.value,
      p_one_less    = ft_less$p.value,
      stringsAsFactors = FALSE
    )
  })
  out_tbl <- bind_rows(out_list)
  out_tbl$p_FDR_two_sided <- p.adjust(out_tbl$p_two_sided, method = "BH")
  out_tbl$p_FDR_one_less  <- p.adjust(out_tbl$p_one_less,  method = "BH")
  out_tbl %>% arrange(p_FDR_one_less, p_one_less)
}
# Run separately for each genotype
fisher_HWT  <- fisher_vs_lexa(df, "HWT")
fisher_K16R <- fisher_vs_lexa(df, "K16R")
# Optional combined table
fisher_all <- bind_rows(fisher_HWT, fisher_K16R)
# Peek
fisher_HWT

fisher_K16R

# View(fisher_all)

```