################################
## Simple example of SAFE     ##
## and generating DAGs for GO ##
## 01-15-08                   ##
## Bill Barry                 ##
################################

## Source code (CHANGE DIRECTORY AS NEEDED)
#source("safe.0108.R")


##### Using the AML/ALL example that is in the SAFE vignette

library(multtest)
library(hu6800)
data(golub)
dimnames(golub)[[1]] <- golub.gnames[,3]


##### NEW: platform and annotate arguments in safe() to build categories

set.seed(12345)
results <- safe(golub, golub.cl, platform="hu6800", annotate="GO.CC",
                min.size=5,max.size=100)
results

C.mat <- getCmatrix(keyword.list= as.list(hu6800GO2ALLPROBES), 
                    present.genes = golub.gnames[,3], min.size=5,
                    max.size=100,GO.ont="CC")


##### Basic local statistics
local <- "default"
## 1) 2-sample experiment
y.vec <- golub.cl
## 2) ANOVA
#y.vec <- rep(1:3,each=13)[-39]
## 3) LM
#y.vec <- 1:38
## 4) paired t-test
#y.vec <- rep(1:19,2)*rep(c(-1,1),each=19)
#local <- "t.paired"


set.seed(12345)
results <- safe(golub, y.vec, C.mat, local=local)
#results



##### Bootstrap-based tests in SAFE

results2 <- safe(golub, y.vec, C.mat, local=local, method="bootstrap",
                args.local = list(one.sided=F,boot.test="t"))
#results2

#results2 <- safe(golub, y.vec, C.mat, local=local, method="bootstrap.q",
#                args.local = list(one.sided=F,boot.test="q"))
#results2



##### Other local statistics

## 6) z.COXPH
results <- safe(golub, 1:38, C.mat, local="z.COXPH", args.local = list(one.sided=F,
                censor=rep(T:F,19)))
results2 <- safe(golub, 1:38, C.mat, local="z.COXPH", args.local = list(one.sided=F,
                censor=rep(T:F,19)), method="bootstrap")
results2 <- safe(golub, 1:38, C.mat, local="z.COXPH", args.local = list(one.sided=F,
                censor=rep(T:F,19),covariate=matrix(rnorm(76),2,38)), method="bootstrap")


#### Confirmation of global pvalues ####
range(results@global.pval)
range(results2@global.pval)
plot(results@global.pval,results2@global.pval)
########################################                        



##### Other forms of annotation

results2 <- safe(golub, golub.cl, platform="hu6800", annotate="KEGG",
                min.size=5,max.size=100,method="bootstrap",
                args.local = list(one.sided=T,boot.test="t"))
results2

results2 <- safe(golub, golub.cl, platform="hu6800", annotate="PFAM",
                min.size=5,max.size=100,method="bootstrap",
                args.local = list(one.sided=T,boot.test="t"))
results2


##### Error terms
results <- safe(golub, golub.cl, C.mat, error="FDR.YB")
results <- safe(golub, golub.cl, C.mat, error="FWER.WY")
results <- safe(golub, golub.cl, C.mat, method="bootstrap", error="FDR.BH")
results <- safe(golub, golub.cl, C.mat, method="bootstrap", error="FWER.Bonf")
results <- safe(golub, golub.cl, C.mat, method="bootstrap", error="FWER.Holm")


##### Output and plots

results                          # Default output
safeplot(results)                # SAFEplot of most significant category
gene.results(results)            # NEW: get gene-specific results
gene.results(results, cat.name = "GO:0005941") 




## NEW: DAG plots (may look better when written to a postscript file)

safedag(results)                 # All categories
safedag(results,filter=1)        # Categores that pass the first color.cutoff
safedag(results,filter=2)
safedag(results,top="GO:0016020", filter=1)  # 


## NEW: other global statistics
results <- safe(golub, golub.cl, C.mat, global="AveDiff")
results <- safe(golub, golub.cl, C.mat, method="bootstrap", global="AveDiff")
results <- safe(golub, golub.cl, C.mat, global="Fisher",
                 args.global = list(genelist.cutoff=2))
results <- safe(golub, golub.cl, C.mat, global="Fisher",
                 args.global = list(genelist.length=200))
results <- safe(golub, golub.cl, C.mat, global="Pearson",
                 args.global = list(genelist.cutoff=2))
results <- safe(golub, golub.cl, C.mat, method="bootstrap", global="Pearson",
                args.global = list(genelist.cutoff=2),
                args.local = list(one.sided=F,boot.test="t"))