###################################################
### chunk number 1: preliminaries
###################################################
library(GeneSelector)


###################################################
### chunk number 2: 
###################################################
data(toydata)
yy <- as.numeric(toydata[1,])
xx <- as.matrix(toydata[-1,])
dim(xx)
table(yy)


###################################################
### chunk number 3: 
###################################################
par(mfrow=c(2,2))
for(i in 1:4) boxplot(xx[i,]~yy, main=paste("Gene", i))


###################################################
### chunk number 4: 
###################################################

ordT <- RankingTstat(xx, yy, type="unpaired")
BaldiLongT <- RankingBaldiLong(xx, yy, type="unpaired")
FoxDimmicT <- RankingFoxDimmic(xx, yy, type="unpaired")
sam <- RankingSam(xx, yy, type="unpaired")
wilcox <- RankingWilcoxon(xx, yy, type="unpaired", pvalues=TRUE)
wilcoxefron <- RankingWilcEbam(xx, yy, type="unpaired")



###################################################
### chunk number 5: 
###################################################

class(ordT)
getSlots("GeneRanking")
str(ordT)



###################################################
### chunk number 6: 
###################################################

show(ordT)
summary(ordT)
toplist(ordT)



###################################################
### chunk number 7: 
###################################################
loo <- GenerateFoldMatrix(xx, yy, k=1)
show(loo)


###################################################
### chunk number 8: 
###################################################

loor_ordT <- GetRepeatRanking(ordT, loo)
loor_BaldiLongT <- GetRepeatRanking(BaldiLongT, loo)
loor_FoxDimmicT <- GetRepeatRanking(FoxDimmicT, loo)
loor_sam <- GetRepeatRanking(sam, loo)
loor_wilcox <- GetRepeatRanking(wilcox, loo)
loor_wilcoxefron <- GetRepeatRanking(wilcoxefron, loo)



###################################################
### chunk number 9: 
###################################################

ex1r_ordT <- GetRepeatRanking(ordT, loo, scheme = "Labelexchange")



###################################################
### chunk number 10: 
###################################################

boot <- GenerateBootMatrix(xx, yy, maxties=3, minclassize=5, repl=30)
show(boot)
boot_BaldiLongT <- GetRepeatRanking(BaldiLongT, boot)



###################################################
### chunk number 11: 
###################################################
noise_ordT <- GetRepeatRanking(ordT, varlist=list(genewise=TRUE, factor=1/10))


###################################################
### chunk number 12: 
###################################################
toplist(loor_ordT, show=FALSE)


###################################################
### chunk number 13: 
###################################################
par(mfrow=c(2,2))
plot(loor_ordT, col="blue", pch=".", cex=2.5)
plot(ex1r_ordT, col="blue", pch=".", cex=2.5)
plot(boot_BaldiLongT, col="blue", pch=".", cex=2.5)
plot(noise_ordT, frac=1/10, col="blue", pch=".", cex=2.5)


###################################################
### chunk number 14: 
###################################################
perturb_ordT <- join(ex1r_ordT, noise_ordT)
show(perturb_ordT)


###################################################
### chunk number 15: 
###################################################

stab_lm_ordT <- GetStabilityLm(loor_ordT, decay="linear")
show(stab_lm_ordT)
stab_lm_BaldiLongT <- GetStabilityLm(loor_BaldiLongT, decay="linear")
show(stab_lm_BaldiLongT)
stab_lm_FoxDimmicT <- GetStabilityLm(loor_FoxDimmicT, decay="linear")
show(stab_lm_FoxDimmicT)
stab_lm_sam <- GetStabilityLm(loor_sam, decay="linear")
show(stab_lm_sam)
stab_lm_wilcox <- GetStabilityLm(loor_wilcox, decay="linear")
show(stab_lm_wilcox)
stab_lm_wilcoxefron <- GetStabilityLm(loor_wilcoxefron, decay="linear")
show(stab_lm_wilcoxefron)



###################################################
### chunk number 16: 
###################################################

ordT_adjpval <- AdjustPvalues(ordT@pval, method="BH")
plot(ordT@pval, ordT_adjpval, xlab="raw p-values", ylab="adjusted p-values")



###################################################
### chunk number 17: 
###################################################

alphaopt <- GetAlpha(ordT@ranking, ordT_adjpval)


plot(1:length(ordT@ranking), 1-exp(-alphaopt*(1:length(ordT@ranking))),
     type="l", xlab="ranks", ylab="")

stab_overlap_ordT <- GetStabilityOverlap(loor_ordT, decay="exponential",
                                          alpha=alphaopt)

plot(stab_overlap_ordT)



###################################################
### chunk number 18: 
###################################################

rs_ordT <- RecoveryScore(loor_ordT, method="BH")
print(rs_ordT)



###################################################
### chunk number 19: 
###################################################

agg_ordT <- AggregateBayes(loor_ordT, stab_lm_ordT, tau=1)
agg_BaldiLongT <- AggregateBayes(loor_BaldiLongT, stab_lm_BaldiLongT, tau=1)
agg_FoxDimmicT <- AggregateBayes(loor_FoxDimmicT, stab_lm_FoxDimmicT, tau=1)
agg_sam <- AggregateBayes(loor_sam, stab_lm_sam, tau=1)
agg_wilcox <- AggregateBayes(loor_wilcox, stab_lm_wilcox, tau=1)
agg_wilcoxefron <- AggregateBayes(loor_wilcoxefron, stab_lm_wilcoxefron, tau=1)



###################################################
### chunk number 20: 
###################################################

agg_simple <- AggregateSimple(loor_BaldiLongT, stab_lm_BaldiLongT,
                              aggregatefun = "mean")


###################################################
### chunk number 21: 
###################################################

statlist  <- list(agg_ordT, agg_BaldiLongT, agg_FoxDimmicT, agg_sam, agg_wilcox,
                  agg_wilcoxefron)


HeatmapMethods(statlist, ind=1:100)



###################################################
### chunk number 22: 
###################################################
ordstat <-  c(4,5,2,3,1,6)


###################################################
### chunk number 23: 
###################################################

gk50 <- GeneSelector(statlist, ind=NULL, indstatistic=ordstat, threshold="user", maxrank=50)



###################################################
### chunk number 24: 
###################################################

gpval <- GeneSelector(statlist, ind=NULL, indstatistic=ordstat, threshold="BH", maxpval=0.05)



###################################################
### chunk number 25: 
###################################################

show(gk50)
show(gpval)
toplist(gpval)

SelectedGenes(gpval)



###################################################
### chunk number 26: 
###################################################

plot(gpval)



###################################################
### chunk number 27: 
###################################################

GeneInfoScreen(gpval, which=30)