setwd("/usr/local/data/supplements/sfp/data/barley")
library(limma)
source("../SFP_function.R")
load("barleyRNA.pm.RData")
load("pheno.RData")

################################################################
allSNP <- read.csv("SNPtable.csv",as.is=T)
temp <- which(duplicated(allSNP$Probe))
allSNP <- allSNP[-temp,]
# drop GP,MX both have changes from probe
str(allSNP <- allSNP[allSNP$SNPgenotype != "GP,MX",])



################################################################
pm.medpolishres <- get_medpolish_res(barleyRNA.pm,ps.n=11)
colnames(pm.medpolishres) <- pheno[,2]
(cl <- factor(substring(colnames(pm.medpolishres),1,1)))
design <- model.matrix(~0+cl);colnames(design) <- sub('cl','',colnames(design));contrast.matrix <- makeContrasts(M-G,levels=design)
fit <- do_limma(pm.medpolishres,design,contrast.matrix)
##save(fit,file="fit.RData")
(fixnum <- summary(fit.data <- decideTests(fit,p=0.001)));table(fit.data!=0)
##save(fit.data,file="fit.data.RData")


probeSNP <- match(as.character(allSNP$Probe),sfpName)

fixnum <- 10504
str(my.pqvalue <- p.adjust(fit$p.value[,1],'fdr'))
(thresh <-quantile(my.pqvalue, fixnum/length(my.pqvalue)))
fit.data <- matrix(rep(0,length(my.pqvalue)))
fit.data <- new("TestResults", array(0, dim(fit$t)))
dimnames(fit.data) <- dimnames(fit$t)
rownames(fit.data) <- sub('--.*$','',rownames(fit.data))
fit.data[my.pqvalue<thresh&fit$t>0] <- 1
fit.data[my.pqvalue<thresh&fit$t<0] <- -1
print(summary(fit.data))
##save(fit.data,file="fit.data.medianpolish_10504.RData",compress=T)

table(snp.geno <- sfp.check <- fit.data[rownames(fit.data)%in%allSNP$Probe,])
snp.geno[snp.geno!=0] <- 0
snp.geno[names(snp.geno)%in%allSNP$Probe[allSNP$SNPgenotype=='GP']] <- 1
snp.geno[names(snp.geno)%in%allSNP$Probe[allSNP$SNPgenotype=='MX']] <- -1
##save(snp.geno,file="snp.geno.RData",compress=T)
################

##load("snp.geno.RData")
sfp.check <-  fit.data[match(names(snp.geno),rownames(fit.data)),]
print(y <- table(snp.geno,sfp.check))
print(summary(y))
cat("Sensitivity", sum(c(y[c(1,3),-2]))/sum(c(y[c(1,3),])), "\n" )
cat("False Negative Rate", 1-sum(c(y[c(1,3),-2]))/sum(c(y[c(1,3),])), "\n") 
cat("False positive Rate", sum(c(y[2,-2]))/sum(c(y[2,])), "\n")
cat("False sequence polymorphism Rate", sum(c(y[2,-2]))/sum(c(y[,c(1,3)])), "\n")