# This script generates all figures in the paper.
# software package: maanova1.2.1, Rmaanova0.96-1
# all figures are made into black and white
# the first number is from the all probeset analysis and the second number is 
# from revision with SNP probesets removed (6-26-06 xq)
#################################################################

setwd("F:/xcui/JAX/MicroarrayAnalysis/Churchill_heterosis_ABF1/Affy/analysis06")

######################### Figure 1 in paper: clustering samples

load("../ftest.sample.yuna.RData")
ftest.sample$obsAnova$Sample.level <- c("AJ","AJ","AJ","AJxB6","AJxB6","AJxB6","B6","B6","B6","B6xAJ","B6xAJ","B6xAJ")
        #idx<-which(ftest.sample$Fs$adjPvalperm<0.05)
        # length(idx)#15861
idx<-which(ftest.sample$Fs$adjPvalperm<0.001)
length(idx)#6148

#HC of samples
cluster <- macluster(ftest.sample$obsAnova, term="Sample", idx.gene=idx, what =  "sample")
par(mfcol=c(1,1))
source("consensusNoTitle.R")
tmp = consensus(cluster, level = 0.8)
      # Note: looks good.
idx<-which(ftest.sample$Fs$adjPvalperm<0.05)
 length(idx)#15847

#HC of samples
cluster <- macluster(ftest.sample$obsAnova, "Sample", idx, what =  "sample")
par(mfcol=c(1,1))
source("consensusNoTitle.R")
consensus(cluster, level = 0.8)
 #HC of samples
cluster <- macluster(ftest.sample$obsAnova, "Sample", idx, what =  "gene")
par(mfcol=c(1,1))
source("consensusNoTitle.R")
consensus(cluster, level = 0.8)
# kmean cluster
 cluster <- macluster(ftest.sample$obsAnova, term="Sample", idx.gene=idx, what =  "gene", method="kmean", kmean.ngroups=7)
consensus(cluster, level = 0.8)


############## plot the cdf of varianceMouse (Figure 2 in paper)

# Use the shrinkage estimators to plot 
load("anova.strainSEp.RData")
library(stepfun)
# shrink the variance components
Shrink.s2.aj= JSshrinker(anova.AJ$S2,  anova.AJ$model$df)
Shrink.s2.b6= JSshrinker(anova.B6$S2,  anova.B6$model$df)
Shrink.s2.ab= JSshrinker(anova.AB$S2,  anova.AB$model$df)
Shrink.s2.ba= JSshrinker(anova.BA$S2,  anova.BA$model$df)

# replot the figure again   
  par(mfcol=c(1,1))
  x <- (1:200)/1000
  tmp <- ecdf(Shrink.s2.aj[,1])
  y <- tmp(x)
  plot(x[2:201],y,type="l", col="black",xlab="Mouse Variance", ylab="Cumulative Fraction", main=" ",xlim=c(0, 0.15), ylim=c(0.8, 1))
  tmp <- ecdf(Shrink.s2.b6[,1])
  y <- tmp(x)
  lines(x[2:201],y, lty=3)
  tmp <- ecdf(Shrink.s2.ab[,1])
  y <- tmp(x)
  lines(x[2:201],y, lty=5)
  tmp <- ecdf(Shrink.s2.ba[,1])
  y <- tmp(x)
  lines(x[2:201],y,lty=4)

  legend(0.1, 0.9, c("AJ","B6","AJxB6","B6xAJ"),  lty = c(1,3,5,4))


########## Figure 3 in paper: heterability ########################################
setwd("F:/xcui/JAX/MicroarrayAnalysis/Churchill_heterosis_ABF1/Affy/analysis06/revise")
load(file="Data.RData")
# subset the whole data set into AB and BA  F1 with the parental lines
DataAB <- subset(Data, arrays=c(1:18))
DataBA <- subset(Data, arrays=cbind(t(c(1:6)),t(c(19:24)),t(c(13:18))))
#Data <- DataAB
Data <- DataBA

# estimate the variances of grp, mouse, and meas.
model.affy2 <- makeModel(DataAB,formula=~Sample+Strain, random=~Sample+Strain)
anova.AB.affy2 <- fitmaanova(DataAB, model.affy2)
load("anova.AB.affy2.RData")
   #anova.AB.affy2$S2.level
   #[1] "Sample" "Strain"
#load("anova.AB.affy2.RData")
# raw variance
s2 <- anova.AB.affy2$S2
h23 <- s2[,2]/(s2[,1]+s2[,2]+s2[,3])
H23 <- h23[which(h23>0.01)]

# shrunken variance
JSs2 <- JSshrinker(anova.AB.affy2$S2,  anova.AB.affy2$model$df)
h24 <- JSs2[,2]/(JSs2[,1]+JSs2[,2]+JSs2[,3])
H24 <- h24[which(h24>0.01)]

median(h23)   # 0.16
median(H23)   # 0.30
median(h24)   # 0.22
median(H24)   # 0.38
length(H24)/length(h24) # 0.7041529
length(H23)/length(h23)  # 0.6989202
 
 # median heritability of heritable genes
 load("ftest.strain.yuna.RData")

# number of significant genes
idx.h = which(ftest.strain$Fs$adjPvalperm<0.05)
median(h24[idx.h]) #0.71
 
# plot them together
a <- hist(H23, breaks=100, main="Before Shrinkage", plot=F, freq=T )
b <- hist(H24, breaks=100,main="After Shrinkage", plot=F, freq=F)

plot(c(a$breaks[1],a$breaks), c(0, a$counts,0), type="s", xlab="Heritability", ylab="Frequency", lwd=1, lty =3, col="black", xlim=c(0, 1))
lines(c(b$breaks[1],b$breaks),c(0,b$counts,0), type="s", lwd=1, lty=1,col="black")
legend(0.6,600, c("naive", "shrinkage"), lty=c(3,1))

# var plot (Figure 7 in presentation)
source("varplot.lim.bw.R")

# look at the variance of mouse and genetic
plot(anova.AB.affy2$S2[,1], anova.AB.affy2$S2[,2])

idx11=which(anova.AB.affy2$S2[,1]>0.000001)
idx12=which(anova.AB.affy2$S2[,2]>0.000001)
length(idx11) # 28539
idx01= intersect(idx11,idx12) #19174
cor((anova.AB.affy2$S2[idx01,1]), (anova.AB.affy2$S2[idx01,2])) # 0.2036203
cor((anova.AB.affy2$S2[idx01,1])^0.5, (anova.AB.affy2$S2[idx01,2])^0.5)# 0.4684162
cor((anova.AB.affy2$S2[idx01,1])^0.1, (anova.AB.affy2$S2[idx01,2])^0.1) # 0.4675783

plot((anova.AB.affy2$S2[idx01,2])^0.1, (anova.AB.affy2$S2[idx01,1])^0.1, pch=".", xlab='mouse varance', ylab="Strain varance", main="Strain and mouse variance")
  
varplot.lim.bw(anova.AB.affy2, xlim=c(0,0.6))
 # plot the shrunken variances
 temp <- anova.AB.affy2
 temp$S2 <- JSs2
 varplot.lim.bw(temp, xlim=c(0,0.6))


 ####### try all data instead of the subset (used in version 8)
load(file="Data.RData")
model.strain <- makeModel(Data, formula=~Sample+Strain, random=~Sample+Strain)
anova.strain.random=fitmaanova(Data,model.strain)
save(anova.strain.random, file="anova.strain.random")
s2 <- anova.strain.random$S2
h23 <- s2[,2]/(s2[,1]+s2[,2]+s2[,3])
H23 <- h23[which(h23>0.01)]

# shrunken variance
JSs2 <- JSshrinker(anova.strain.random$S2,  anova.strain.random$model$df)
h24 <- JSs2[,2]/(JSs2[,1]+JSs2[,2]+JSs2[,3])
H24 <- h24[which(h24>0.01)]
median(h23)  #0.16  # 0.1617763
median(H23)  #0.26   # 0.2581955
median(h24)  #0.20   # 0.200348
median(H24)  #0.31   #0.3060319
length(H24)/length(h24) # 0.7412252  #0.7439736
length(H23)/length(h23)  # 0.7367242  #0.7392943
length(h23)-length(H23) # 11874  #11756
length(h24)-length(H24) #11676   #11545
length(which(h23>0.5)) #6503     #6557
length(which(h23>0.5))/length(h23) #0.1441875    #0.1454106
length(which(h24>0.5)) #8131  #8197
length(which(h24>0.5))/length(h24)#0.1802843  #0.1817799

# plot them together
a <- hist(H23, breaks=100, main="Before Shrinkage", plot=F, freq=T )
b <- hist(H24, breaks=100,main="After Shrinkage", plot=F, freq=F)

plot(c(a$breaks[1],a$breaks), c(0, a$counts,0), type="s", xlab="Heritability", ylab="Frequency",  lwd=1, lty =3, col="black", xlim=c(0, 1))
lines(c(b$breaks[1],b$breaks),c(0,b$counts,0), type="s", lwd=1, lty=1,col="black")
legend(0.6,600, c("naive", "shrinkage"), lty=c(3,1))

# var plot (Figure 7 in presentation)
source("varplot.lim.bw.R")
varplot.lim.bw(anova.strain.random, xlim=c(0,0.6))
# plot the shrunken variances
temp <- anova.strain.random
temp$S2 <- JSs2
varplot.lim.bw(temp, xlim=c(0,0.6))
     # comments: the shrinkage gain becomes much smaller

  ####### Try the BA data
Data <- DataBA

# estimate the variances of grp, mouse, and meas.
model.affy2 <- makeModel(DataBA,formula=~Sample+Strain, random=~Sample+Strain)
anova.BA.affy2 <- fitmaanova(DataBA, model.affy2)
#anova.AB.affy2$S2.level
#[1] "Sample" "Strain"

# raw variance
s2 <- anova.BA.affy2$S2
h23 <- s2[,2]/(s2[,1]+s2[,2]+s2[,3])
H23 <- h23[which(h23>0.01)]

# shrunken variance
JSs2 <- JSshrinker(anova.BA.affy2$S2,  anova.BA.affy2$model$df)
h24 <- JSs2[,2]/(JSs2[,1]+JSs2[,2]+JSs2[,3])
H24 <- h24[which(h24>0.01)]

median(h23)   # 0.16
median(H23)   # 0.29
median(h24)   # 0.22
median(H24)   # 0.37
length(H24)/length(h24) # 0.7031551
length(H23)/length(h23) # 0.697612

    # The question is: should we do separately or combined?  The two separate set show similar result with higher percent of heritability
    

################# Figure 5 in paper:  degree of dominance plot ####################################
load(file="ftest.strain.yuna.RData")
idx2 <- which(ftest.strain$Fs$adjPvalperm<0.01)#7871

load(file="anova.strain.RData")
# plot the four panels together
par(mfcol=c(2,2))
anova <- anova.strain
d <- anova$Strain[,2]-0.5*(anova$Strain[,1]+anova$Strain[,3])
a <-0.5*(anova$Strain[,1]-anova$Strain[,3])
cor(abs(a),abs(d))#0.4080499
dd <-d/a #  It could be huge because of small a
#histogram
temp <- dd
limt <- 10 # the boundary of x 
length(temp[temp< -limt])# 1509
length(temp[temp > limt])#  1641
idx <- which(temp< -limt)
temp[idx]=-limt
idx <- which(temp>limt)
temp[idx]=limt
A <- hist( temp, breaks=100,plot=F, freq=T )
plot(c(A$breaks[1],A$breaks), c(0, A$density,0), type="s", xlab="d/a",ylab="Frequency", lwd=1, lty =1, col="black", xlim=c(-limt,limt), main="AxB")
abline(v=1,col="black", lty=3)
abline(v=-1,col="black", lty=3)
abline(v=0,col="black", lty=1)

# plot the scatter plot
 idx <- which(abs(dd)<10)
plot(dd[idx],a[idx], pch=".", cex=2, xlab="d/a", ylab="a ", main="AxB", col="grey")
idx <- intersect(idx,idx2)
points(dd[idx],a[idx],pch=".", cex=2, col="black")

# the  plot for BA
anova <- anova.strain
d <- anova$Strain[,4]-0.5*(anova$Strain[,1]+anova$Strain[,3])
a <-0.5*(anova$Strain[,1]-anova$Strain[,3])
cor(abs(a),abs(d))#0.4080499
dd <-d/a #  It could be huge because of small a
#histogram
temp <- dd
limt <- 10 # the boundary of x 
length(temp[temp< -limt])#
length(temp[temp > limt])#
idx <- which(temp< -limt)
temp[idx]=-limt
idx <- which(temp>limt)
temp[idx]=limt
A <- hist( temp, breaks=100,plot=F, freq=T )
plot(c(A$breaks[1],A$breaks), c(0, A$density,0), type="s", xlab="d/a", ylab="Frequency", lwd=1, lty =1, col="black", xlim=c(-limt,limt), main="BxA")
abline(v=1,col="black", lty=3)
abline(v=-1,col="black", lty=3)
abline(v=0,col="black", lty=1)

# plot the scatter plot
 idx <- which(abs(dd)<10)
plot(dd[idx],a[idx], pch=".", cex=2, xlab="d/a", ylab="a ", main="BxA", col="grey")
idx <- intersect(idx,idx2)
points(dd[idx],a[idx],pch=".", cex=2, col="black")


# the scater plot of d/a vs a for AxB
idx <- which(abs(dd)<5)
plot(dd[idx],a[idx], pch=".", cex=3, xlab="d/a", ylab=" ", main="AxB")
idx <- intersect(idx,idx2)
points(dd[idx],a[idx],pch=".", cex=3, col="red")


################# Figure 3 in paper: ad plot

load(file="FourthRun/ttest.ab.yuna.RData")
load(file="FourthRun/ttest.f1ab.yuna.RData")
load(file="FourthRun/ttest.f1ba.yuna.RData")

anova.strain = ttest.ab$obsAnova
# additive and dominant effect
d.AB <- anova.strain$Strain[,2]-0.5*(anova.strain$Strain[,1]+anova.strain$Strain[,3])
a <-0.5*(anova.strain$Strain[,1] - anova.strain$Strain[,3])
idx.a <-which(ttest.ab$Fs$adjPvalperm<0.05)   #8363  #8846 #8950
idx.d.AB <- which(ttest.f1ab$Fs$adjPvalperm<0.05)  #612 #321 #317
idx.d.BA =   which(ttest.f1ba$Fs$adjPvalperm<0.05) #23   #20  #24
length(union(idx.d.AB, idx.d.BA)) #615   #324 #321
#overdominance
load("FourthRun/ttest.AF1ab.yuna.RData")
load("FourthRun/ttest.BF1ab.yuna.RData")
idxAF1 <- which(ttest.AF1ab$Fs$adjPvalperm<0.05)# A-AB 1578   #1121
idxBF1 <- which(ttest.BF1ab$Fs$adjPvalperm<0.05)# B-AB 3746   #3625
idxint <- intersect(idxAF1, idxBF1) #518    #406   #470

idxF1out <- which(abs(d.AB)>abs(a))#24635   #24710 #24671
idxoverd <- intersect(idxint, idxF1out)#176  #130  #158

## ad plot 
par(mfcol=c(1,1))
plot(a[idx.a],d.AB[idx.a],pch=".", cex=1, xlab="a", ylab="d", col="black")
points(a[idx.d.AB],d.AB[idx.d.AB],pch=21, cex=0.8,col="black" )
abline(0,1)
abline(0,0)
points(a[idxoverd],d.AB[idxoverd],pch=19, cex=0.5, col="grey")
abline(v=0)
abline(0,-1)
text(2,0,"pure additive")
text(0,2,"over dominance")
text(1.5,1.5,"A dominance")
text(-1.5,1.5,"B dominance")


### ############ histograms of p values of AJ vs B6, dominant effect of AB and over dominant of AB similar to A.

load(file="FourthRun/ttest.ab.yuna.RData")
pval.AvsB <- ttest.ab$Fs$Pvalperm
load("FourthRun/ttest.f1ab.yuna.RData")
pval.dAB <- ttest.f1ab$Fs$Pvalperm
 anova.strain = ttest.ab$obsAnova
idx1 <- which(anova.strain$Strain[,2]>anova.strain$Strain[,1] & anova.strain$Strain[,1]>anova.strain$Strain[,3])   # to compare with A
idx2 <- which(anova.strain$Strain[,2]<anova.strain$Strain[,1] & anova.strain$Strain[,1]<anova.strain$Strain[,3])  # compare with A

pval <- ttest.AF1ab$Fs$Pvalperm
pval.ABsimA <- pval[c(idx1,idx2)]#length = 13233

load("FourthRun/ttest.BF1ab.yuna.RData")
idx12 <- which(anova.strain$Strain[,2]>anova.strain$Strain[,3] & anova.strain$Strain[,3]>anova.strain$Strain[,1])   # to compare with B
idx22 <- which(anova.strain$Strain[,2]<anova.strain$Strain[,3] & anova.strain$Strain[,3]<anova.strain$Strain[,1]) 
pval <- ttest.AF1ab$Fs$Pvalperm
pval.ABsimB <- pval[c(idx12,idx22)]

length(union(idx12, union(idx22, union(idx1, idx2))))

# step function of histograms for p values
hist.AvsB <- hist(pval.AvsB,breaks=20, plot=F, freq=F)#A-B
hist.dAB <- hist(pval.dAB,breaks=20, plot=F, freq=F)#A-BA
hist.ABsimA <- hist(pval.ABsimA,breaks=20, plot=F, freq=F)#BA sim A

par(mfcol=c(1,1))
# plot the step lines for the three histograms
plot(c(hist.AvsB$breaks[1],hist.AvsB$breaks), c(0, hist.AvsB$density,0), type="s", xlab="p values", ylab="Density", lwd=1.5, lty =1, col="black", xlim=c(0, 1))
lines(c(hist.dAB$breaks[1],hist.dAB$breaks),c(0,hist.dAB$density,0), type="s", lwd=1, lty=3,col="black")
lines(c(hist.ABsimA$breaks[1],hist.ABsimA$breaks),c(0,hist.ABsimA$density,0), type="s", lwd=1.5, lty=1,col="grey")
legend(0.3, 15, legend=c("additive","dominant","overdominant"), col=c("black","black","grey"), lwd=c(1,1,1), lty=c(1,3,1))

# (not used in paper)
pval.overdom = c(pval.ABsimB, pval.ABsimA)
hist.overdom <- hist(pval.overdom,breaks=20, plot=F, freq=F)
# plot the step lines for the three histograms  and the AB compare to B
plot(c(hist.AvsB$breaks[1],hist.AvsB$breaks), c(0, hist.AvsB$density,0), type="s", xlab="p-values", ylab="Density", lwd=1.5, lty =1, col="black", xlim=c(0, 1))
lines(c(hist.overdom$breaks[1],hist.overdom$breaks),c(0,hist.overdom$density,0), type="s", lwd=1.5, lty=1,col="grey")
lines(c(hist.dAB$breaks[1],hist.dAB$breaks),c(0,hist.dAB$density,0), type="s", lwd=1, lty=3,col="black")
legend(0.6, 15, legend=c("additive           ","overdominant","dominant"), col=c("black","grey","black"), lwd=c(1,1,1), lty=c(1,1,3))

# try sepearate panels  (used in paper)
par(mfrow=c(1,2))
hist(pval.AvsB, breaks=50, main="additive",xlab="p-values")
hist(pval.dAB, breaks=50,  main="dominance",xlab="p-values", ylab=" ")


 ########### figure 6 RT-PCR result
 
 ################ Real-time PCR result analysis ###########################
###########################################################################

# read in the data files
RTraw.1 <- read.table(file="../RealTimePCR/04TQMA-001.gac10-01results.editxq.txt", header = T, sep = "\t", quote = "\"'", dec = ".")

RTraw.2 <- read.table(file="../RealTimePCR/04TQMA-001.gac10-04resultsv3.editxq.txt", header = T, sep = "\t", quote = "\"'", dec = ".")

mean.1 <- rowMeans(RTraw.1[,3:5]) # average the three replicate wells on the same machine
CTdiff.1 <- mean.1-rep(mean.1[25:36],times=4)   # subtract the gus reference
CTdiff.1 <- matrix(CTdiff.1, nrow=16,ncol=3, byrow=T)  # format the means into a matrix with each strain on one row.
mean.11 <- rowMeans(CTdiff.1)  # mean for each strain
var.1 <- apply(RTraw.1[,3:5],1, var)   # variance for each sample
Svar.1 <- matrix(var.1, nrow=16, ncol=3, byrow=T)   # format to each strain on one row
Svar.11 <- apply(Svar.1, 1, sum)                    # sum of the variances
Svar.12 <- sqrt((Svar.11+rep(Svar.11[9:12], times=4))/9)    # overall variance

tmp <- cbind(mean.11, Svar.12)
# make to the same order as the table (A, AxB, BxA, B)
tmp1 <- tmp[c(2,4,3,1,6,8,7,5,10,12,11,9,14,16,15,13),]
write.table(tmp1, file = "RTmean.1.txt", append = FALSE, quote = F, sep = "\t ", eol = "\n", na = "NA", dec = ".", row.names = F, col.names = F, qmethod =  "double")

mean.2 <- apply(RTraw.2[,3:5],1, mean)
CTdiff.2 <- mean.2-rep(mean.2[13:24],times=4)
CTdiff.2 <- matrix(CTdiff.2, nrow=16,ncol=3, byrow=T)
mean.22 <- apply(CTdiff.2, 1, mean)

var.1 <- apply(RTraw.2[,3:5],1, var)
Svar.1 <- matrix(var.1, nrow=16, ncol=3, byrow=T)
Svar.11 <- apply(Svar.1, 1, sum)
Svar.22 <- sqrt((Svar.11+rep(Svar.11[5:8], times=4))/9)

tmp <- cbind(mean.22, Svar.12)
# make to the same order as the table (A, AxB, BxA, B)
tmp2 <- tmp[c(1,2,4,3,5,6,8,7,9,10,12,11,13,14,16,15),]
write.table(tmp, file = "RTmean.2.txt", append = FALSE, quote = F, sep = "\t ", eol = "\n", na = "NA", dec = ".", row.names = F, col.names = F, qmethod =  "double")

### The corresponding microarray estimates
load("Data.RData")
probe.id <- c("1418148_at","1417461_at","1428004_at","1421363_at","1425614_x_at","1421179_at")

idxRT <- match(probe.id, Data$cloneid)

# load the ANOVA results for each strain
load( file="anova.strainSep.RData")
#measurement error for AJ
error.AJ <- sqrt(anova.AJ$S2[idxRT,1]/3+anova.AJ$S2[idxRT,2]/6)
#0.09493079 0.08613173 0.05695361 0.06309428 0.19247635 0.10065438
error.B6 <- sqrt(anova.B6$S2[idxRT,1]/3+anova.B6$S2[idxRT,2]/6)
#0.15811209 0.07131115 0.05701190 0.18318328 0.06410770 0.06672011
error.AB <- sqrt(anova.AB$S2[idxRT,1]/3+anova.AB$S2[idxRT,2]/6)
#0.03834914 0.08234975 0.06865848 0.08619570 0.27389717 0.02929364
error.BA <- sqrt(anova.BA$S2[idxRT,1]/3+anova.BA$S2[idxRT,2]/6)
# 0.06732608 0.03781023 0.07432088 0.05532177 0.01057519 0.04557384

# organize them into a matrix that match the RT-PCR 
RTresult= rbind(tmp1[1:8,], tmp1[13:16,], tmp2[1:4,], tmp2[9:16,])
# relative to B
RTresult[,1]= RTresult[,1]-rep(RTresult[c(4,8,12,16,20,24),1],each=4)
maresult=as.vector(rbind(error.AJ, error.AB, error.BA, error.B6))
load("anova.strain.RData")
mamean = anova.strain$Strain[idxRT,]+ matrix(rep(anova.strain$G[idxRT],each=4), nrow=6, ncol=4, byrow=T)
mamean = mamean[,c(1,2,4,3)]
maresult = cbind( as.vector(t(mamean)),maresult)
maresult[,1]= maresult[,1]-rep(maresult[c(4,8,12,16,20,24),1],each=4)

######## plot the Results
#pdf(file="RTpcrFig.pdf",height=4, width=5)
par(mfrow=c(2,3))
txt =c("Abhd1", "Cap1","3300001GO2RIK","Cyp2c39","H2-D1","Ttyh2")
for (i in 1:6){
mx = max(c((-RTresult[((4*i-3):(4*i)),1]-RTresult[((4*i-3):(4*i)),2]),  (-RTresult[((4*i-3):(4*i)),1]+RTresult[((4*i-3):(4*i)),2]) ,(maresult[((4*i-3):(4*i)),1]-RTresult[((4*i-3):(4*i)),2]), (maresult[((4*i-3):(4*i)),1]+RTresult[((4*i-3):(4*i)),2])))
mi = min(c((-RTresult[((4*i-3):(4*i)),1]-RTresult[((4*i-3):(4*i)),2]),  (-RTresult[((4*i-3):(4*i)),1]+RTresult[((4*i-3):(4*i)),2]) ,(maresult[((4*i-3):(4*i)),1]-RTresult[((4*i-3):(4*i)),2]), (maresult[((4*i-3):(4*i)),1]+RTresult[((4*i-3):(4*i)),2])))
#par(bty="n")
plot(c(1:4),-RTresult[((4*i-3):(4*i)),1], pch=".", ylim=c(mi,mx), ylab="log2(FC)", xlab="",  lty=3, main=txt[i], axes=F)
axis(1,at=c(1,2,3,4), labels=c("A","AxB","BxA","B"))
axis(2,labels=T)
lines (c(1:4),-RTresult[((4*i-3):(4*i)),1], pch=".", lty=3)
segments(c(1:4), (-RTresult[((4*i-3):(4*i)),1]-RTresult[((4*i-3):(4*i)),2]) , c(1:4), (-RTresult[((4*i-3):(4*i)),1]+RTresult[((4*i-3):(4*i)),2]) )
lines(c(1.001,2.001,3.001,4.001),maresult[((4*i-3):(4*i)),1], pch=".", lty=1)
 segments(c(1.001,2.001,3.001,4.001), (maresult[((4*i-3):(4*i)),1]-RTresult[((4*i-3):(4*i)),2]) , c(1.001,2.001,3.001,4.001), (maresult[((4*i-3):(4*i)),1]+RTresult[((4*i-3):(4*i)),2]) )
 }