library("stats")

#method 1 student's t test
#pen name of William Sealy Gosset
#soporific drug
?sleep
sleep
str(sleep)
summary(sleep)
mean(sleep$extra)
mean(sleep$extra[sleep$group == 1])
mean(sleep$extra[sleep$group == 2])
sd(sleep$extra[sleep$group == 1])
sd(sleep$extra[sleep$group == 2])
tapply(sleep$extra, sleep$group, mean)
tapply(sleep$extra, sleep$group, sd)
tapply(sleep$extra, sleep$group, summary)
table(sleep[,c(2)])

#following t distribution
x <- -100:100 / 20
y <- dnorm(x)
plot(y ~ x)
for (i in 1:10) lines(dt(x, i) ~ x, col = i)

#t test
t.test(sleep$extra[sleep$group == 1])
t.test(sleep$extra[sleep$group == 2])

#one-tailed test
t.test(sleep$extra[sleep$group == 1], alternative = "g")
t.test(sleep$extra[sleep$group == 2], alternative = "g")

#compare two drugs
t.test(sleep$extra[sleep$group == 1], sleep$extra[sleep$group == 2])

#paired t test
?t.test
t.test(sleep$extra[sleep$group == 1], sleep$extra[sleep$group == 2], paired = T)
t.test(sleep$extra[sleep$group == 1] - sleep$extra[sleep$group == 2])

shapiro.test(sleep$extra[sleep$group == 1])
shapiro.test(sleep$extra[sleep$group == 2])

#gene expression example
setwd("D:/IGDB/courses/SDC/2017/classIV")
GSE6536 <-  read.table("GSE6536_series_matrix_trio_ensembl_avg_list.txt", header = T)
str(GSE6536)
levels(GSE6536$chrom)

#male female difference on gene expression?

Gene <- "ENSG00000000003"
t.test(GSE6536$offspring_expr[GSE6536$gene == Gene & GSE6536$sex == "d"], 
       GSE6536$offspring_expr[GSE6536$gene == Gene & GSE6536$sex == "s"])

t.test(GSE6536$father_expr[GSE6536$gene == Gene], 
       GSE6536$mother_expr[GSE6536$gene == Gene],
       paired = T)



#you need to worry about if the data follows normal distribution
#to test if the data follow normal distribution

# ks.test(x = GSE6536$offspring_expr[GSE6536$gene == Gene & GSE6536$sex == "d"], 
#         y = rnorm(length(GSE6536$offspring_expr[GSE6536$gene == Gene & GSE6536$sex == "d"]), 
#                   mean(GSE6536$offspring_expr[GSE6536$gene == Gene & GSE6536$sex == "d"]),
#                   sd(GSE6536$offspring_expr[GSE6536$gene == Gene & GSE6536$sex == "d"]))
#         )
# 
# ks.test(x = GSE6536$offspring_expr[GSE6536$gene == Gene & GSE6536$sex == "d"], 
#         y = "pnorm",
#         mean(GSE6536$offspring_expr[GSE6536$gene == Gene & GSE6536$sex == "d"]),
#         sd(GSE6536$offspring_expr[GSE6536$gene == Gene & GSE6536$sex == "d"])
# )

#quiz
shapiro.test(GSE6536$mother_expr[GSE6536$gene == Gene])

#method 2 rank sum test
t.test(rnorm(20), rnorm(20))

t.test(c(rnorm(16), rnorm(4, 20)), rnorm(20))
wilcox.test(c(rnorm(16), rnorm(4, 20)), rnorm(20))
hist(c(rnorm(16), rnorm(4, 20)))
plot(density(c(rnorm(16), rnorm(4, 20))), ylim = c(0, 0.5))
lines(density(rnorm(20)), col ="red")

#quiz
wilcox.test(sleep$extra[sleep$group == 1], sleep$extra[sleep$group == 2], paired = T)

# wilcox.test(GSE6536$father_expr[GSE6536$gene == Gene], 
#        GSE6536$mother_expr[GSE6536$gene == Gene],
#        paired = T)
# wilcox.test(GSE6536$father_expr[GSE6536$gene == Gene], 
#        GSE6536$mother_expr[GSE6536$gene == Gene],
#        paired = F)
# #method 3 permutation test
# (A <- c(rnorm(16), rnorm(4, 20)))
# (B <- rnorm(20))
# 
# t.test(A, B)
# (Union <- c(A, B))
# length(Union)
# (D <- mean(A) - mean(B))
# 
# #similar to t test 
# L <- 10000
# Diff <- vector(mode = "numeric", length = L)
# for (i in 1:L) {
#   S <- sample(Union, 40, replace = F) 
#   Diff[i] <- mean(S[1:20]) - mean(S[21:40])
# }
# (D <- mean(A) - mean(B))
# (P_val <- sum(Diff >= D | Diff <= -D) / L)
# t.test (A, B)
# 
# #similar to wilcox test
# for (i in 1:L) {
#   S <- sample(Union, 40, replace = F) 
#   Diff[i] <- median(S[1:20]) > median(S[21:40])
# }
# (P_val <- sum(Diff) / L)
# 
# for (i in 1:L) {
#   S <- sample(Union, 40, replace = F) 
#   Diff[i] <- mean(S[1:20]) > mean(S[21:40])
# }
# (P_val <- sum(Diff) / L)
# 
# wilcox.test (A, B)


# ?df
# x <- 1:1000 / 100
# d <- df(x, 1, 1)
# plot(d ~ x)
# for (i in 2:10){
#   d <- df(x, 1, i)
#   lines(d ~ x)
# }
# 
# x <- 1:1000 / 100
# d <- df(x, 1, 1)
# plot(d ~ x)
# for (i in 2:10){
#   d <- df(x, i, 1)
#   lines(d ~ x)
# }
# 
# x <- 1:1000 / 100
# d <- df(x, 1, 1)
# plot(d ~ x)
# for (i in 2:10){
#   d <- df(x, i, 10)
#   lines(d ~ x)
# }

#method 3 bootstrapping
# Bootstrapping <- vector (length =10000)
# Father <- GSE6536$father_expr[GSE6536$gene == Gene]
# Mother <- GSE6536$mother_expr[GSE6536$gene == Gene]
# for (i in 1:10000) Bootstrapping[i] <- sum(sample(Father, replace = T))-sum(sample(Mother, replace = T))
# sum(Bootstrapping > 0)


#method 4 ks.test
t.test(c(rnorm(10), rnorm(10, 10)), rnorm(20, 5))
wilcox.test(c(rnorm(10), rnorm(10, 10)), rnorm(20, 5))

#but you know that it should be different
plot(density(c(rnorm(10), rnorm(10, 10))), ylim = c(0, 0.5))
lines(density(rnorm(20, 5)), col ="red")

#f test
var.test(c(rnorm(10), rnorm(10, 10)), rnorm(20, 5))

#sleep
var.test(sleep$extra[sleep$group == 1], sleep$extra[sleep$group == 2])

#quiz
var.test(GSE6536$offspring_expr[GSE6536$gene == Gene & GSE6536$sex == "d"], 
         GSE6536$offspring_expr[GSE6536$gene == Gene & GSE6536$sex == "s"])

#compare the distribution by ks. test
ks.test(c(rnorm(10), rnorm(10, 10)), rnorm(20, 5))
plot(ecdf(c(rnorm(10), rnorm(10, 10))))
lines(ecdf(rnorm(20, 5)))

ks.test(sleep$extra[sleep$group == 1], sleep$extra[sleep$group == 2])

#quiz
ks.test(GSE6536$father_expr[GSE6536$gene == Gene], 
        GSE6536$mother_expr[GSE6536$gene == Gene],
        paired = T)

