#R script to verify SFPs in Nipponbare and 9311

library(limma)
targets <- readTargets(file="targets.txt")
myfun <- function(x) as.numeric(x$Flags > -49.5) #excludes spots flagged "Not Found" and "Bad"
RG <- read.maimages(targets, source="genepix.median", wt.fun=myfun, other.columns="Flags") #read in file
x11()
plotMA(RG,1)
x11()
plotMA(RG,2)
x11()
plotMA(RG,3)
x11()
plotMA(RG,4)
RG.within <- normalizeWithinArrays(RG, method="loess",bc.method="none", span=0.02, weights=RG$weights)
MA<-RG.within
#MA <- normalizeBetweenArrays(RG.within, method="none",bc.method="none", weights=RG$weights)
x11()
plotMA(MA,1)
x11()
plotMA(MA,2)
x11()
plotMA(MA,3)
x11()
plotMA(MA,4)

size=dim(MA)[1]/3
design <- cbind(JvsI=c(1,1,-1,-1),DyeEffect=1)
w <- modifyWeights(MA$weights, MA$genes$Name, c("Blank","Repeat"), c(0,0))
corfit <- duplicateCorrelation(MA, design, ndups = 3, spacing=size,block=NULL, weights=w)
fit1 <- lmFit(MA, design, ndups= 3, spacing=size, block=NULL, cor = corfit$consensus, weights=w)
fit1 <- eBayes(fit1)
ordinary.t <- -abs(fit1$coef[,1]) / fit1$stdev.unscaled[,1] / fit1$sigma
ordinary.p<-pt(ordinary.t,df=fit1$df.residual)
out=cbind(fit1$genes$Name,RG$other$Flags[1:size,],RG$other$Flags[(size+1):(size*2),],RG$other$Flags[(size*2+1):(size*3),]
,fit1$stdev.unscaled,fit1$coefficients,fit1$Amean,fit1$p.value,ordinary.p)
write.table(out, file="Nipponbare_9311.txt", quote=FALSE, sep="\t", row.names=FALSE, col.names=TRUE)
top1 <- topTable(fit1, coef="GvsB", number = 10, adjust = "BH")

#R script for BSA

library(limma)
targets <- readTargets(file="targets_IsoIV.txt")
myfun <- function(x) as.numeric(x$Flags > -99.5) #excludes spots flagged "Bad"
RG <- read.maimages(targets, source="genepix.median", wt.fun=myfun)
spottypes <- readSpotTypes(file="Spottypes.txt")
RG$genes$Status <- controlStatus(spottypes, RG)
w <- modifyWeights(RG$weights, RG$genes$Status, c("Blank","Repeat","Saturated","Badprint"), c(0,0,0,0))
#x11()
#plotMA(RG,1)
#x11()
#plotMA(RG,2)
#x11()
#plotMA(RG,3)
#x11()
#plotMA(RG,4)
RG.within <- normalizeWithinArrays(RG, method="loess",bc.method="none", weights=w)
MA <- normalizeBetweenArrays(RG.within, method="none",bc.method="none", weights=w)
#x11()
#plotMA(MA,1)
#x11()
#plotMA(MA,2)
#x11()
#plotMA(MA,3)
#x11()
#plotMA(MA,4)


size=dim(MA)[1]/3
design <- cbind(TreatmentvsControl=c(-1,1,-1,1),DyeEffect=1)
w <- modifyWeights(MA$weights, MA$genes$Status, c("Blank","Repeat","Saturated","Same","Badprint","Positive","Negative","Bad"), c(0,0,0,0,0,0,0,0))
corfit <- duplicateCorrelation(MA, design, ndups = 3, spacing=size,block=NULL, weights=w)
fit1 <- lmFit(MA, design, ndups= 3, spacing=size, block=NULL, cor = corfit$consensus, weights=w)
fit1 <- eBayes(fit1)
ordinary.t <- -abs(fit1$coef[,"TreatmentvsControl"]) / fit1$stdev.unscaled[,"TreatmentvsControl"] / fit1$sigma
ordinary.p<-pt(ordinary.t,df=fit1$df.residual)

top1 <- topTable(fit1, coef="TreatmentvsControl", number = 20, adjust = "BH")
top1
write.table(top1, file="IsoIV_top.xls", quote=FALSE, sep="\t", row.names=FALSE)

fit1$p.adjust <- p.adjust(fit1$p.value[,"TreatmentvsControl"], method = "BH")

write.table(fit1, file="IsoIV_fit.xls", quote=FALSE, sep="\t", row.names=FALSE)

############################################################################################

#extract probe information
i.fit1.sfp <- fit1$genes$Status=="SFP"
fit1.sfp <- cbind(fit1$p.adjust[i.fit1.sfp],fit1$coefficients[i.fit1.sfp,"TreatmentvsControl"])
fit1.sfp.names <- fit1$genes$Name[i.fit1.sfp]
positions<-matrix(nrow=length(fit1.sfp.names),ncol=3)
for (i in 1:length(fit1.sfp.names)) 
	{
	positions[i,1]<-as.numeric(unlist(strsplit(fit1.sfp.names[i],"_"))[1])
	positions[i,2]<-as.numeric(unlist(strsplit(fit1.sfp.names[i],"_"))[2])
	if (unlist(strsplit(fit1.sfp.names[i],"_"))[3] == "J") {positions[i,3]<-fit1.sfp[i,2]/abs(fit1.sfp[i,2])} else {positions[i,3]<-(fit1.sfp[i,2]/abs(fit1.sfp[i,2])*-1)}
	}

#plot -log of adjusted p-values * direction of effect for all chromosomes
x11()
numrows<-6
par(mfcol=c(numrows,2),mar=c(.5,0,0,3),oma=c(2,2,2,2),bty="n")
for(i in 1:max(positions[,1])){
p <-positions[,1] ==i
positions.chr <- positions[p,]
positions.chr[,2] <- positions.chr[,2]/1000000
fit.chr <- log10(fit1.sfp[p,1])*positions.chr[,3]
plot(positions.chr[,2],fit.chr,main=NULL,xlim=c(0,max(positions[,2])/1000000),ylim=c(-5,10),xaxp=c(0,max(positions[,2])/1000000,10),yaxt="n",xaxt="n")
lines(c(0,max(positions.chr[,2])),c(0,0))
#chr.max<- max(positions.chr[,2])
chr.tick<-1000000
axis(2,at=c(-5:10),pos=-0.5)
lines(lowess(positions.chr[,2],fit.chr,f=1/6))
if (i==numrows) {axis(1,at=c(0:(chr.max)))}
}
#chr.max<- max(positions[,2])/1000000
chr.max<-30
axis(1,at=c(0:(chr.max)))

#Plot chromosome 12 only
x11()
par(mar=c(.5,0,0,3),oma=c(2,2,2,2),bty="n")
p <-positions[,1] == 12
positions.chr <- positions[p,]
positions.chr[,2] <- positions.chr[,2]/1000000
fit.chr <- log10(fit1.sfp[p,1])*positions.chr[,3]
plot(positions.chr[,2],fit.chr,main=NULL,xlim=c(0,max(positions.chr[,2])),ylim=c(-5,10),xaxp=c(0,max(positions.chr[,2]),10),yaxt="n",xaxt="n")
lines(c(0,max(positions.chr[,2])),c(0,0))
#chr.max<- max(positions.chr[,2])
chr.tick<-1000000
axis(2,at=c(-5:10),pos=-0.5)
lines(lowess(positions.chr[,2],fit.chr,f=1/6))
axis(1,at=c(0:(chr.max)))
#chr.max<- max(positions[,2])/1000000
chr.max<-30
axis(1,at=c(0:(chr.max)))


#plot coefficients for all chromosomes
x11()
numrows<-6
par(mfcol=c(numrows,2),mar=c(.5,0,0,3),oma=c(2,2,2,2),bty="n")
for(i in 1:max(positions[,1])){
p <-positions[,1] ==i
positions.chr <- positions[p,]
positions.chr[,2] <- positions.chr[,2]/1000000
fit.chr <-fit1.sfp[p,2]*positions.chr[,3]
plot(positions.chr[,2],fit.chr,main=NULL,xlim=c(0,max(positions[,2])/1000000),ylim=c(-0.4,0.4),xaxp=c(0,max(positions[,2])/1000000,10),yaxt="n",xaxt="n")
lines(c(0,max(positions.chr[,2])),c(0,0))
#chr.max<- max(positions.chr[,2])
chr.tick<-1000000
axis(2,at=c(-5:10),pos=-0.5)
lines(lowess(positions.chr[,2],fit.chr,f=1/3))
if (i==numrows) {axis(1,at=c(0:(chr.max)))}
}
#chr.max<- max(positions[,2])/1000000
chr.max<-30
axis(1,at=c(0:(chr.max)))


#Plot chromosome 12 only
x11()
par(mar=c(.5,0,0,3),oma=c(2,2,2,2),bty="n")
p <-positions[,1] == 12
positions.chr <- positions[p,]
positions.chr[,2] <- positions.chr[,2]/1000000
fit.chr <-fit1.sfp[p,2]*positions.chr[,3]
plot(positions.chr[,2],fit.chr,main=NULL,xlim=c(0,max(positions.chr[,2])),ylim=c(-0.4,0.4),xaxp=c(0,max(positions.chr[,2]),10),yaxt="n",xaxt="n")
lines(c(0,max(positions.chr[,2])),c(0,0))
#chr.max<- max(positions.chr[,2])
chr.tick<-1000000
axis(2,at=c(-5:10),pos=-0.5)
lines(lowess(positions.chr[,2],fit.chr,f=1/3))
axis(1,at=c(0:(chr.max)))
#chr.max<- max(positions[,2])/1000000
chr.max<-30
axis(1,at=c(0:(chr.max)))


#R script for polymorphism survey
library(limma)
accessionlist=read.delim("inputaccessions.txt", header=FALSE)



targets <- readTargets(file="targets_all.txt")



#myfun <- function(x) as.numeric(x$Flags > -49.5) #excludes spots flagged "Not Found" and "Bad"
myfun <- function(x) as.numeric(x$Flags > -50.5)  #excludes spots flagged "Bad"
#myfun <- function(x) as.numeric(x$Flags > -100.5) #allows spots flagged "Bad"
RG <- read.maimages(targets, source="genepix", wt.fun=myfun) #read in file
#show(RG)
#summary(RG$R)
spottypes <- readSpotTypes(file="SpotTypes_controls_small.txt")
#spottypes <- readSpotTypes(file="SpotTypes.txt")
RG$genes$Status <- controlStatus(spottypes, RG)
#x11() 
#plotMA(RG)
w <- modifyWeights(RG$weights, RG$genes$Status, c("Blank","Bad"), c(0,0))
poscontrols <- RG$genes$Status=="Positive"
negcontrols <- RG$genes$Status=="Negative"
repcontrols <- RG$genes$Status=="Repeat"
smallcontrols <- RG$genes$Status=="small"
blankcontrols <- RG$genes$Status=="Blank"
#normcontrols <- (poscontrols + negcontrols + repcontrols + blankcontrols)==1
#normcontrols <- (poscontrols + negcontrols)==1
normcontrols <- smallcontrols==1
#MA<-MA.RG(RG)
MA <- normalizeWithinArrays(RG, method="loess", controlspots=normcontrols,weights=w,bc.method="none")
#MA <- normalizeWithinArrays(RG, method="loess", weights=w,bc.method="none")
#x11()
#plotMA(MA)

MA.between <- normalizeBetweenArrays(RG, method="Gquantile")


#i.RG.sfp <- RG.b$genes$Status=="SFP"


rm(scorematrix)
for (accession in 1:ncol(MA.between)){

MA.accession<-MA.between[,accession]
size=dim(MA.accession)[1]/3
MA.N <- cbind(MA.accession[1:size,],MA.accession[(size+1):(size*2),],MA.accession[(size*2+1):(size*3),])
pcontrl <- MA.N$genes$Status=="Positive"
ncontrl <- MA.N$genes$Status=="Negative"
sfpcontrl <- MA.N$genes$Status=="SFP"
w <- modifyWeights(MA.N$weights, MA.N$genes$Status, c("Blank","Bad"), c(0,0))
design <- c(1, 1, 1)
fit1 <- lmFit(MA.N, design, weights=w)
fit1 <- eBayes(fit1)


i.fit1.sfp <- fit1$genes$Status=="SFP"
fit1.sfp.names <- fit1$genes$Name[i.fit1.sfp]


if (exists("scorematrix")){scorematrix<-cbind(scorematrix,fit1$coefficients[i.fit1.sfp])
} else {scorematrix<-data.frame(fit1$coefficients[i.fit1.sfp], row.names=fit1.sfp.names)}


}

names(scorematrix)<-accessionlist[,1]


#write.table(t(scorematrix), file="newsmall.txt", quote=FALSE, sep="\t", row.names=TRUE, col.names=TRUE)

write.table(scorematrix, file="genotypes.xls", quote=FALSE, sep="\t", row.names=TRUE, col.names=TRUE)