#simulation for understanding the false positive
Test_n <- 2000
P_val_sim <- vector(length = Test_n)
for (i in 1:Test_n){
  P_val_sim[i] <- t.test(rnorm(20), rnorm(20))$p.value
}
sum(P_val_sim <= 0.05) / Test_n
hist(P_val_sim)


# in real data
str(GSE6536)
Genes <- levels(GSE6536$gene)
# for (i in 1: length(Genes)) {
#   wilcox.test(GSE6536$father_expr[GSE6536$gene == Genes[i]], 
#               GSE6536$mother_expr[GSE6536$gene == Genes[i]],
#               paired = T)
# }

Test_n <- 200
P_val <- vector(length = Test_n)
for (i in 1: Test_n) {
  P_val[i] <- wilcox.test(GSE6536$father_expr[GSE6536$gene == Genes[i]], 
              GSE6536$mother_expr[GSE6536$gene == Genes[i]],
              paired = T)$p.value
}
sum(P_val < 0.05) / Test_n
hist(P_val)

#Empirical Cumulative Distribution Function
plot(ecdf(P_val))
lines(ecdf(P_val_sim), col = 2)

#introudction to FDR
source("http://bioconductor.org/biocLite.R")
biocLite("qvalue")

library(qvalue) 

qobj_sim <- qvalue(P_val_sim)
sum(qobj_sim$qvalues <= 0.05)

qobj <- qvalue(P_val)
sum(qobj$qvalues <= 0.05)


#quiz
#Females have 2 X chromosomes, while males have only one.  
#For all the genes on chromosome X, how many are 
#differentially expressed between males and females?
GSE6536_X_10 <- GSE6536[GSE6536$chrom == "X" & GSE6536$avg >= 10, ]
Genes_X_10 <- levels(factor(GSE6536_X_10$gene))
P_val_X_10 <- vector(length = length(Genes_X_10))
for (i in 1: length(Genes_X_10)) {
  P_val_X_10[i] <- wilcox.test(GSE6536_X_10$father_expr[GSE6536_X_10$gene == Genes_X_10[i]], 
                          GSE6536_X_10$mother_expr[GSE6536_X_10$gene == Genes_X_10[i]],
                          paired = T)$p.value
}
sum(P_val_X_10 < 0.05)
qobj_X_10 <- qvalue(P_val_X_10)
sum(qobj_X_10$qvalues <= 0.05)

P_val_X_10_less <- vector(length = length(Genes_X_10))
for (i in 1: length(Genes_X_10)) {
  P_val_X_10_less[i] <- wilcox.test(GSE6536_X_10$father_expr[GSE6536_X_10$gene == Genes_X_10[i]], 
                          GSE6536_X_10$mother_expr[GSE6536_X_10$gene == Genes_X_10[i]],
                          paired = T, alternative = "l")$p.value
}
sum(P_val_X_10_less < 0.05)
qobj_X_10_less <- qvalue(P_val_X_10_less)
sum(qobj_X_10_less$qvalues <= 0.05)

P_val_X_10_greater <- vector(length = length(Genes_X_10))
for (i in 1: length(Genes_X_10)) {
  P_val_X_10_greater[i] <- wilcox.test(GSE6536_X_10$father_expr[GSE6536_X_10$gene == Genes_X_10[i]], 
                                    GSE6536_X_10$mother_expr[GSE6536_X_10$gene == Genes_X_10[i]],
                                    paired = T, alternative = "g")$p.value
}
sum(P_val_X_10_greater < 0.05)
qobj_X_10_greater <- qvalue(P_val_X_10_greater)
sum(qobj_X_10_greater$qvalues <= 0.05)

binom.test(6, 8)
