set.seed(42)
daysOfWeek <- c("Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday")
day.born <- sample(daysOfWeek, 1000, replace=T)
born.tuesday <- day.born == "Tuesday"

cancer <- sample(c(TRUE, FALSE), 1000, replace=T, prob=c(0.393, 0.607))	

tbl <- table(born.tuesday, cancer)
fisher.test(tbl)

simulate_study <- function() {
  daysOfWeek <- c("Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday")
  day.born <- sample(daysOfWeek, 1000, replace=T)
  born.tuesday <- day.born == "Tuesday"
  cancer <- sample(c(TRUE, FALSE), 1000, replace=T, prob=c(0.393, 0.607))
  tbl <- table(born.tuesday, cancer)
  test <- fisher.test(tbl)
  test$p.value
}


set.seed(12345)
pvals <- numeric(0)
for (i in 1:100) {
  pvals[[i]] <- simulate_study()
}

hist(pvals)

sum(pvals < 0.05)
sum(pvals < 0.2)


bonf.pvals <- p.adjust(pvals, method='bonferroni')
sum(bonf.pvals < 0.05)
sum(bonf.pvals < 0.2)


set.seed(1)
rand.data <- matrix(rnorm(12000), ncol=12)

get.p.value <- function(x) {
  test <- t.test(x[1:6], x[7:12])
  test$p.value
}

pvals <- numeric(0)
for (row in 1:nrow(rand.data)) {
  pvals[[row]] <- get.p.value(rand.data[row, ])
}

sum(pvals < 0.05)
bonf.pvals <- p.adjust(pvals, method='bonferroni')
sum(bonf.pvals < 0.05)

rand.data[951:1000, 7:12] <- rnorm(6*50, mean=2)
pvals <- numeric(0)
for (row in 1:nrow(rand.data)) {
  pvals[[row]] <- get.p.value(rand.data[row, ])
}
sum(pvals < 0.05)
sum(pvals[951:1000] < 0.05)

bonf.pvals <- p.adjust(pvals, method="bonferroni")
sum(bonf.pvals < 0.05)

BH.pvals <- p.adjust(pvals, method='BH')
sum(BH.pvals < 0.05)
sum(BH.pvals < 0.1)
sum(BH.pvals[951:1000] < 0.1)
sum(BH.pvals < 0.25)
sum(BH.pvals[951:1000] < 0.25)

library(tidyverse)
airway <- read_csv("https://denvirlab.marshall.edu/CS505/airway_cs505.csv")
gene_names <- airway %>% pull(gene_name) 

expressionMatrix <- as.matrix(airway[,-1])

pvalues <- numeric(0)
ctrlColumns <- c(1,3,5,7)
trtColumns <- c(2,4,6,8)
for (i in 1:nrow(expressionMatrix)) {
  tt <- t.test(
    expressionMatrix[i, ctrlColumns], 
    expressionMatrix[i, trtColumns], 
    paired=T
  )
  pvalues[[i]] <- tt$p.value
}

pvalues <- tibble(gene = gene_names, p.value=pvalues)
pvalues <- pvalues %>% mutate(padj = p.adjust(pval, method='BH'))

filter(pvalues, pval < 0.05)
filter(pvalues, padj < 0.05)

# Tidyverse approach:

library(broom)

airway <- read_csv("https://denvirlab.marshall.edu/CS505/airway_cs505.csv") %>% 
  pivot_longer(names_to = "Sample", values_to = "Expression", cols= -gene_name) %>%
  separate(col=Sample, into=c("Patient", "Treatment"), sep='_') %>%
  nest_by(gene_name, .key="expression_data") %>% 
  mutate(tidy(t.test(Expression ~ Treatment, paired=TRUE, data=expression_data))) %>% 
  ungroup %>% # We must apply p.adjust to the whole set of p-values, not the p-values row by row
  mutate(padj = p.adjust(p.value, method="BH"))

filter(airway, p.value < 0.05)
filter(airway, padj < 0.05)