########################################################################
# This script is for reanalyzing the F1 data considering the SNPs
# between AJ and B6 in the revision of the paper.
# 6/21/2006 8:52:11 AM
# software package: rmaanova0.98-7 (maybe 1.00 because installed 1.00, but # version function says 0.98, and matestHyunaF.R
###########################################################################


########### test for the SNP removal #################################### 

# read in the data
raw = ReadAffy()
result = rma(raw)
 write.exprs(result, file='rma.expr.affy.txt') 
 
# modify the CDF file to remove the SNP containing probes 
library(ProbeFilter)
xyList<-scan("F:/xcui/MicroarrayDataSets/SNPremoveSoftware/AJB6/probe_snp_AJ_B6_noGenotype.txt", what=list('numeric', 'numeric', 'char', 'char', 'numeric'), sep='\t', skip=1)
#data@cdfName<-'Mm430_Mm_ENTREZG_7'
ProbeFilter(raw, xyList=xyList)
rslt<-rma(raw)
setwd("F:/xcui/JAX/MicroarrayAnalysis/Churchill_heterosis_ABF1/Affy/Analysis06/revise")
write.exprs(rslt, file='rma.expr.customcdf.txt')  

################## read in preprocessed data 

raw.data.custom <- read.madata(datafile="rma.expr.customcdf.txt",designfile="AffyDesignFinal.txt",pmt=2,spotflag=F, cloneid=1, header=T)
     # already logarithm.
raw.data.affy <- read.madata(datafile="rma.expr.affy.txt",designfile="AffyDesignFinal.txt",pmt=2,spotflag=F, cloneid=1, header=T)

# find the probe sets that are removed from the custom CDF
comp = raw.data.affy$cloneid %in% raw.data.custom$cloneid
 idx = which(comp==FALSE)
 raw.data.affy$cloneid[idx] #[1] "1416809_at"   "1435628_x_at" "1435697_a_at" "1436822_x_at" "1438050_x_at" "1438936_s_at" "1452452_at"   "1452540_a_at"
 
 # look at the numerical differences between the two CDF
 plot(raw.data.affy$data[-idx,1], raw.data.custom$data[,1]) # quite bit of difference.  all probe sets change numerically.
 
# read in annotation (description, probeID=cloneid, unigeneID). 
anno <- read.delim("../description.txt")
    # Note: anno order differs from cloneid in raw.data


###############  redo analysis on data with SNP probes removed ####################

Data <- createData(raw.data.custom)
id = Data$cloneid

# replace the data fieled with raw data to avoid log again
Data$data <- raw.data.custom$data
save(Data,file="Data.RData")
# recalculate the column means
Data$colmeans =  colMeans(Data$data)
save(Data, file="Data.RData")

############### compute the variance components using fitmaanova ################ 

DataAJ <- subset(Data,arrays=c(1:6))
DataB6 <- subset(Data,arrays=c(13:18))
DataAB <- subset(Data,arrays=c(7:12))
DataBA <- subset(Data,arrays=c(19:24))

model.var <- makeModel(DataAJ,formula=~Sample, random=~Sample)
anova.AJ <- fitmaanova(DataAJ,model.var)
anova.AB <- fitmaanova(DataAB,model.var)
anova.B6 <- fitmaanova(DataB6,model.var)
anova.BA <- fitmaanova(DataBA,model.var)
save(anova.AJ, anova.AB, anova.B6, anova.BA, file="anova.strainSep.RData")

###################################### fitmaanova ####################

model.affy <- makeModel(DataAB,formula=~Sample+Strain, random=~Sample)
save(DataAB,DataBA, model.affy, file="DataModelABBA.affy.RData")

anova.AB.affy <- fitmaanova(DataAB, model.affy)
anova.BA.affy <- fitmaanova(DataBA, model.affy)
save(anova.AB.affy,anova.BA.affy, file="anova.ABBA.affy.RData")

load(file="anova.ABBA.affy.RData")

# 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 in the list of probesets affected by SNP
    snp.list = read.delim("probeSetsABsnp.txt",header=F, sep="\t")
    snp.list = as.character(snp.list[[1]])
  # load in the original Data for comparison
  load(file="../Data.RData")
  id.org = Data$cloneid  

######################### F test of the 12 samples for clustering the samples  ##########

source("matestHyunaF.R")
Data <- createData(raw.data.custom)
Data$data <- raw.data.custom$data
# recalculate the column means
Data$colmeans =  colMeans(Data$data)

model <- makeModel(Data,formula=~Sample)
ftest.sample <- matest(Data,model,term="Sample", n.perm=100,nnodes=1, shuffle="sample", test.method=c(1,0,1,1))
 ftest.sample = adjPval(ftest.sample, "adaptive")
  save(ftest.sample, file="ftest.sample.yuna.RData")
  
  idx  <- which(ftest.sample$Fs$adjPvalperm<0.01)
  length(idx)  #10464  #10372
  idx  <- which(ftest.sample$Fs$adjPvalperm<0.005)
  length(idx) #8994   #8919
 idx  <- which(ftest.sample$Fs$adjPvalperm<0.05)
  length(idx)  #15861 #15649
  idx  <- which(ftest.sample$Fs$adjPvalperm<0.1)
  length(idx) #19888  #19693
  
 idx<-which(ftest.sample$Fs$Pvalperm<0.00001)
 length(idx)#3024  #2618
idx<-which(ftest.sample$Fs$Pvalperm<0.0001)
 length(idx)#5402 #5352
idx<-which(ftest.sample$Fs$Pvalperm<0.001)
 length(idx)#8399 #8362

ftest.sample$obsAnova$Sample.level <- c("AJ","AJ","AJ","AJxB6","AJxB6","AJxB6","B6","B6","B6","B6xAJ","B6xAJ","B6xAJ")

  volcano(ftest.sample, threshold=c(0.001,0.0001,0.001,0.001),
             method=c("fdr","fdr","fdr","fdr"), title="Volcano Plot",
             highlight.flag=TRUE, onScreen=TRUE)
  
## compare the clustering
      idx<-which(ftest.sample$Fs$adjPvalperm<0.001)      
      length(idx)#6148 #6084
      load(file="../ftest.sample.yuna.RData")
      idx.org = which(ftest.sample$Fs$adjPvalperm<0.001)      
      length(intersect(setdiff( id.org[idx.org],id[ idx]), snp.list))  #83
      length(setdiff( id.org[idx.org],id[ idx]))  #111
      length(setdiff( id[ idx],id.org[idx.org]))  #47
      length(intersect(setdiff( id[ idx],id.org[idx.org]), snp.list)) #45
      
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: The same structure as before with one bootstrapping number increase.


################### F test for mouse ########################3 
model <- makeModel(Data, formula=~Sample+Strain)
ftest.mouse <- matest(Data,model,term="Sample", n.perm=100,nnodes=1, shuffle="sample")
 ftest.mouse = adjPval(ftest.mouse, "adaptive")
  save(ftest.mouse, file="ftest.mouse.yuna.RData")
  idx<-which(ftest.mouse$Fs$Pvalperm<0.001)
 length(idx)#1356   #1239
 idx  <- which(ftest.mouse$Fs$adjPvalperm<0.05)
  length(idx)  #1666 $1452

###############  F test for overall group or genotype effect Using Hyuna's new perm ################################ Note: the permutations were run on cheaha.

library(maanova)
load(file="Data.RData")
model.strain <- makeModel(Data, formula=~Sample+Strain, random=~Sample)
ftest.strain <- matest(Data,model.strain,term="Strain", n.perm=100, nnodes=1, shuffle="sample")
ftest.strain = adjPval(ftest.strain, "adaptive")
save(ftest.strain, file="ftest.strain.yuna.RData")

load("ftest.strain.yuna.RData")

# number of significant genes
length(which(ftest.strain$Fs$adjPvalperm<0.05))# 7871 #7819
length(which(ftest.strain$Fs$Pvalperm<0.001))  #3131  #3110

# overlap with other tests at FDR0.05
idx.1 = which(ftest.strain$Fs$adjPvalperm<0.05)  #7871
idx.2 = which(ftest.mouse$Fs$adjPvalper<0.05) #1666
length(intersect(idx.1, idx.2))   #268

########## matest for a and d effects  ################################

##### compare AB and BA vs mean of the parental strains #
   # Note: strain.level order: A, AxB, B, BxA
# contrast matrix
#      A     AxB  B    BxA
#[1,]  1.0    0 -1.0    0
#[2,]  0.5   -1  0.5    0
#[3,]  0.5    0  0.5   -1
#[4,]  0.0    1  0.0   -1

#load(file="Data.RData")
model.strain <- makeModel(Data, formula=~Sample+Strain, random=~Sample)
model.strain <- makeModel(Data, formula=~Sample+Strain, random=~Sample)
ttest.ab <- matest(Data,model.strain,term="Strain", Contrast=matrix(c(1,0,-1,0),1,byrow=T),test.method=c(1,0,0,1),n.perm=100, nnodes=1, shuffle="sample")

ttest.ab = adjPval(ttest.ab, "adaptive")
save (ttest.ab, file="ttest.ab.yuna.RData")
load(file="ttest.ab.yuna.RData")
length(which(ttest.ab$Fs$adjPvalperm<0.05))# 8363   #  8846
length(which(ttest.ab$Fs$Pvalperm<0.001))  #3769   #3959
idx = which(ttest.ab$Fs$adjPvalperm<0.05)
fd.wSNP = ttest.ab$obsAnova$Strain[,1]-ttest.ab$obsAnova$Strain[,3]
load(file="../ttest.ab.yuna.RData")
fd.woSNP = ttest.ab$obsAnova$Strain[,1]-ttest.ab$obsAnova$Strain[,3]
idx.org = which(ttest.ab$Fs$adjPvalperm<0.05)
length(intersect(setdiff( id.org[idx.org],id[ idx]), snp.list))  #95
length(setdiff( id.org[idx.org],id[ idx]))  #96
# look at the level of A vs B to see whether these are all A<B
idx.tmp = match( setdiff( id.org[idx.org],id[ idx]), id.org)
com.tmp =ttest.ab$obsAnova$Strain[idx.tmp, 1]- ttest.ab$obsAnova$Strain[idx.tmp, 3]
length(which(com.tmp<0))/length(com.tmp)# 0.4895833

#### compare the fole change of A to B before and after removing SNP 
diff.snp = fd.wSNP[match(snp.list,id.org)]- fd.woSNP[match(snp.list,id)]
length(which(diff.snp>0))/length(snp.list) # 0.5065161
#### end

length(setdiff( id[ idx],id.org[idx.org]))  #579
length(intersect(setdiff( id[ idx],id.org[idx.org]), snp.list)) #143

length(intersect(snp.list, id.org[idx]))  #768
# compare the two F1s
model.strain <- makeModel(Data, formula=~Sample+Strain, random=~Sample)
ttest.f1s <- matest(Data,model.strain,term="Strain", Contrast=matrix(c(0,1,0,-1),1,byrow=T),test.method=c(1,0,0,1),n.perm=100, nnodes=1, shuffle="sample")
ttest.f1s = adjPval(ttest.f1s, "adaptive")
save (ttest.f1s, file="ttest.f1s.yuna.RData")
 load("ttest.f1s.yuna.RData")
 length(which(ttest.f1s$Fs$adjPvalperm<0.05))# 755   #137
length(which(ttest.f1s$Fs$Pvalperm<0.001)) #820      #539
       idx = which(ttest.f1s$Fs$adjPvalperm<0.05)
      load(file="../ttest.f1s.yuna.RData")
      idx.org = which(ttest.f1s$Fs$adjPvalperm<0.05)
      length(intersect(setdiff( id.org[idx.org],id[ idx]), snp.list))  #57
      length(setdiff( id.org[idx.org],id[ idx]))  #620
      length(setdiff( id[ idx],id.org[idx.org]))  #2
      length(intersect(setdiff( id[ idx],id.org[idx.org]), snp.list))  #2

# test to compare the two F1s and the two parents
  model.strain <- makeModel(Data, formula=~Sample+Strain, random=~Sample)
ttest.ab.f1s <- matest(Data,model.strain,term="Strain", Contrast=matrix(c(1,-1,1,-1),1,byrow=T),test.method=c(1,0,0,1),n.perm=100, nnodes=1, shuffle="sample")
ttest.ab.f1s = adjPval(ttest.ab.f1s, "adaptive")
save (ttest.ab.f1s, file="ttest.ab.f1.yuna.RData")
 load("ttest.ab.f1.yuna.RData")
  length(which(ttest.ab.f1$Fs$adjPvalperm<0.05)) #  # 41
        idx =  which(ttest.ab.f1$Fs$adjPvalperm<0.05)
        load("../ttest.ab.f1.yuna.RData")
        idx.org =  which(ttest.ab.f1$Fs$adjPvalperm<0.05) #68
        length(intersect(setdiff( id.org[idx.org],id[ idx]), snp.list))  #6
        length(setdiff( id.org[idx.org],id[ idx]))  #64
        length(setdiff( id[ idx],id.org[idx.org]))  #64
        length(intersect(setdiff( id[ idx],id.org[idx.org]), snp.list))  #2

# test for dominance of AB 
 ttest.f1ab <- matest(Data,model.strain,term="Strain", Contrast=matrix(c(0.5,-1,0.5,0),1,byrow=T),test.method=c(1,0,0,1),n.perm=100, nnodes=1, shuffle="sample")
   # Note: strain.level order: A, AxB, B, BxA
ttest.f1ab = adjPval(ttest.f1ab, "adaptive")
save (ttest.f1ab, file="ttest.f1ab.yuna.RData")
load ( file="ttest.f1ab.yuna.RData")
length(which(ttest.f1ab$Fs$adjPvalperm<0.05)) #  #321
idx =  which(ttest.f1ab$Fs$adjPvalperm<0.05)
load ( file="../ttest.f1ab.yuna.RData")
idx.org =  which(ttest.f1ab$Fs$adjPvalperm<0.05)  #612
      length(intersect(setdiff( id.org[idx.org],id[ idx]), snp.list))  #29
      length(setdiff( id.org[idx.org],id[ idx]))  #294
      length(setdiff( id[ idx],id.org[idx.org]))  #3
      length(intersect(setdiff( id[ idx],id.org[idx.org]), snp.list))  #3

 # test for dominance of BA 
   ttest.f1ba <- matest(Data,model.strain,term="Strain", Contrast=matrix(c(0.5,0.5,-1),1,byrow=T),test.method=c(1,0,0,1),n.perm=100, nnodes=1, shuffle="sample")
  # Note: strain.level order: A, AxB, B, BxA
  ttest.f1ba = adjPval(ttest.f1ba, "adaptive")
  save (ttest.f1ba, file="ttest.f1ba.yuna.RData")
  load ( file="ttest.f1ba.yuna.RData")
  length(which(ttest.f1ba$Fs$adjPvalperm<0.05))# 23 #20
  idx = which(ttest.f1ba$Fs$adjPvalperm<0.05)
  load ( file="../ttest.f1ba.yuna.RData")
  idx.org = which(ttest.f1ba$Fs$adjPvalperm<0.05) 
      length(intersect(setdiff( id.org[idx.org],id[ idx]), snp.list))  #2
      length(setdiff( id.org[idx.org],id[ idx]))  #3
      length(setdiff( id[ idx],id.org[idx.org]))  #0
      length(intersect(setdiff( id[ idx],id.org[idx.org]), snp.list))  #0

# overlap between sig genes between F1s and dominant genes
  load("ttest.f1ab.yuna.RData")
  load("ttest.f1ba.yuna.RData")
  load("ttest.f1s.yuna.RData")
  idx.tmp1=   which(ttest.f1ab$Fs$adjPvalperm<0.05)
  idx.tmp12=   which(ttest.f1ba$Fs$adjPvalperm<0.05)
    
  idx.tmp2 = which(ttest.f1s$Fs$adjPvalperm<0.05)
  length(intersect(union(idx.tmp1,idx.tmp12), idx.tmp2))#265  #71
  length(union(union(idx.tmp1,idx.tmp12), idx.tmp2))       #   #390
    
####################### ttest for F1 vs each parental strain ########

library(maanova)
source(file="matestHyunaF.R")
load(file="Data.RData")
model.strain <- makeModel(Data, formula=~Sample+Strain, random=~Sample)

# A vs 2 AB
ttest.AF1ab <- matest(Data,model.strain,term="Strain", Contrast=matrix(c(1,-1,0,0),1,byrow=T),n.perm=100, nnodes=1, shuffle="sample")
ttest.AF1ab = adjPval(ttest.AF1ab)
save (ttest.AF1ab, file="ttest.AF1ab.yuna.RData")

# A vs 2 BA
ttest.AF1ba <- matest(Data,model.strain,term="Strain", Contrast=matrix(c( 1,0,0,-1),1,byrow=T),n.perm=100, nnodes=1, shuffle="sample")
ttest.AF1ba = adjPval(ttest.AF1ba)
save (ttest.AF1ba, file="ttest.AF1ba.yuna.RData")
# B vs AB
ttest.BF1ab <- matest(Data,model.strain,term="Strain", Contrast=matrix(c(0,-1,1,0),1,byrow=T),n.perm=100, nnodes=1, shuffle="sample")
 ttest.BF1ab = adjPval(ttest.BF1ab, "adaptive")
save (ttest.BF1ab, file="ttest.BF1ab.yuna.RData")
# B vs BA
ttest.BF1ba <- matest(Data,model.strain,term="Strain", Contrast=matrix(c(0,0,1,-1),1,byrow=T),n.perm=100, nnodes=1, shuffle="sample")
 ttest.BF1ba = adjPval(ttest.BF1ba, "adaptive")
save (ttest.BF1ba, file="ttest.BF1ba.yuna.RData")

################# identify overdominant genes

 load("ttest.AF1ba.yuna.RData")
 load("ttest.AF1ab.yuna.RData")
 load("ttest.BF1ba.yuna.RData")
 load("ttest.BF1ab.yuna.RData")
 
 #load(file="anova.strain.RData")
 anova.strain = ttest.AF1ab$obsAnova
 # for ab
  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]) 
  idx.ab.a = which(ttest.AF1ab$Fs$adjPvalperm[c(idx1,idx2)]<0.05)#100  #67
  
  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])
  idx.ab.b = which(ttest.BF1ab$Fs$adjPvalperm[c(idx12,idx22)]<0.05) #88  #83
  idx.ab = union(idx.ab.a, idx.ab.b)#188  #150
  
        #  get in the original
        load(file="../anova.strain.RData")
        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])
        load("../ttest.AF1ab.yuna.RData") 
        idx.ab.a.org = which(ttest.AF1ab$Fs$adjPvalperm[c(idx1,idx2)]<0.05)#100 
        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]) 
        load("../ttest.BF1ab.yuna.RData")
        idx.ab.b.org = which(ttest.BF1ab$Fs$adjPvalperm[c(idx12,idx22)]<0.05) 
        idx.ab.org =  union(idx.ab.a.org, idx.ab.b.org) #188
        length(intersect(setdiff( id.org[idx.ab.org],id[ idx.ab]), snp.list))  #15
        length(setdiff( id.org[idx.ab.org],id[ idx.ab]))  #184
        length(setdiff( id[ idx.ab],id.org[idx.ab.org]))  #146
        length(intersect(setdiff( id[ idx.ab],id.org[idx.ab.org]), snp.list))  #10
      
  # for ba
   anova.strain = ttest.AF1ba$obsAnova
  idx1 <- which(anova.strain$Strain[,4]>anova.strain$Strain[,1] & anova.strain$Strain[,1]>anova.strain$Strain[,3])   # to compare with A
  idx2 <- which(anova.strain$Strain[,4]<anova.strain$Strain[,1] & anova.strain$Strain[,1]<anova.strain$Strain[,3]) 
  idx.ba.a = which(ttest.AF1ba$Fs$adjPvalperm[c(idx1,idx2)]<0.05)#6   #2
  
  idx12 <- which(anova.strain$Strain[,4]>anova.strain$Strain[,3] & anova.strain$Strain[,3]>anova.strain$Strain[,1])   # to compare with B
  idx22 <- which(anova.strain$Strain[,4]<anova.strain$Strain[,3] & anova.strain$Strain[,3]<anova.strain$Strain[,1])
  idx.ba.b = which(ttest.BF1ba$Fs$adjPvalperm[c(idx12,idx22)]<0.05) #27
  idx.ba = union(idx.ba.a, idx.ba.b)#33    #25
  length(union(idx.ab, idx.ba)) #221      #175
  
  # get the original
    load("../ttest.AF1ba.yuna.RData")
    anova.strain = ttest.AF1ba$obsAnova
    idx1 <- which(anova.strain$Strain[,4]>anova.strain$Strain[,1] & anova.strain$Strain[,1]>anova.strain$Strain[,3])   # to compare with A
    idx2 <- which(anova.strain$Strain[,4]<anova.strain$Strain[,1] & anova.strain$Strain[,1]<anova.strain$Strain[,3]) 
    idx.ba.a = which(ttest.AF1ba$Fs$adjPvalperm[c(idx1,idx2)]<0.05)#6   #2
    load("../ttest.BF1ba.yuna.RData")
    idx12 <- which(anova.strain$Strain[,4]>anova.strain$Strain[,3] & anova.strain$Strain[,3]>anova.strain$Strain[,1])   # to compare with B
    idx22 <- which(anova.strain$Strain[,4]<anova.strain$Strain[,3] & anova.strain$Strain[,3]<anova.strain$Strain[,1])
    idx.ba.b = which(ttest.BF1ba$Fs$adjPvalperm[c(idx12,idx22)]<0.05) #27
    idx.ba.org = union(idx.ba.a, idx.ba.b)#33   
    length(intersect(setdiff( id.org[idx.ba.org],id[ idx.ba]), snp.list))  #1
    length(setdiff( id.org[idx.ba.org],id[ idx.ba]))  #33
    length(setdiff( id[ idx.ba],id.org[idx.ba.org]))  #25
    length(intersect(setdiff( id[ idx.ba],id.org[idx.ba.org]), snp.list))  #1


############### estimating the proportion of differentially expressed genes

# try q-value package
library(qvalue)
  q.sample = qvalue(ftest.sample$Fs$Pvalperm)
  1-q.sample$pi0# 0.735526   # 0.7356779
  q.a = qvalue(ftest.mouse$Fs$Pvalperm)
  1-q.a$pi0# 0     #0
  q.a = qvalue(ftest.strain$Fs$Pvalperm)
  1-q.a$pi0# 0.67  #0.6732143
  q.a = qvalue(ttest.ab$Fs$Pvalperm)
  1-q.a$pi0# 0.55    # 0.5513194
  q.a = qvalue(ttest.f1ab$Fs$Pvalperm)
  1-q.a$pi0# 0.30     # 0.2907216
  q.a = qvalue(ttest.f1ba$Fs$Pvalperm)
  1-q.a$pi0# 0.19   #0.1654382
  q.a = qvalue(ttest.f1s$Fs$Pvalperm)
  1-q.a$pi0# 0.31   #0.3026897

 # write out the pvalues for mix-o-matic
    out =  ftest.sample$Fs$Pvalperm
    write.table(out, file = "p.f.sample.txt", append = FALSE, quote = F, sep = "\t ", eol = "\n", na = "NA", dec = ".", row.names = F, col.names = F, qmethod =  "double") 
    out =  ftest.mouse$Fs$Pvalperm
    write.table(out, file = "p.f.mouse.txt", append = FALSE, quote = F, sep = "\t ", eol = "\n", na = "NA", dec = ".", row.names = F, col.names = F, qmethod =  "double")
    out =  ftest.strain$Fs$Pvalperm
    write.table(out, file = "p.f.strain.txt", append = FALSE, quote = F, sep = "\t ", eol = "\n", na = "NA", dec = ".", row.names = F, col.names = F, qmethod =  "double")     
    out =  ttest.ab$Fs$Pvalperm
    write.table(out, file = "p.t.ab.txt", append = FALSE, quote = F, sep = "\t ", eol = "\n", na = "NA", dec = ".", row.names = F, col.names = F, qmethod =  "double") 
    out =  ttest.f1ab$Fs$Pvalperm
    write.table(out, file = "p.t.f1ab.txt", append = FALSE, quote = F, sep = "\t ", eol = "\n", na = "NA", dec = ".", row.names = F, col.names = F, qmethod =  "double")
    out =  ttest.f1ba$Fs$Pvalperm
    write.table(out, file = "p.t.f1ba.txt", append = FALSE, quote = F, sep = "\t ", eol = "\n", na = "NA", dec = ".", row.names = F, col.names = F, qmethod =  "double") 
    out =  ttest.f1s$Fs$Pvalperm
    write.table(out, file = "p.t.f1s.txt", append = FALSE, quote = F, sep = "\t ", eol = "\n", na = "NA", dec = ".", row.names = F, col.names = F, qmethod =  "double")   

################ 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]) )
   }

 ################# check wheterh there are SNP in the validated genes
 
 snp.list = read.delim("probeSetsABsnp.txt",header=F, sep="\t")
 # match(probe.id,as.character(snp.list[[1]]))
# NA NA NA NA NA NA

##### Compare the expression summary before and after customizing the CDF
setwd("F:/xcui/JAX/MicroarrayAnalysis/Churchill_heterosis_ABF1/Affy/Analysis06/revise")
load(file="Data.RData")
Data.snp = Data
load(file="../Data.RData")
snp.list = read.delim("probeSetsABsnp.txt",header=F, sep="\t")
snp.list = as.character(snp.list[[1]])

idx.snp = match(snp.list, Data$cloneid)
idx.nosnp = setdiff((1:length(Data$cloneid)), idx.snp)

idx.2 = match(Data$cloneid[idx.nosnp], Data.snp$cloneid)

plot(Data$data[idx.nosnp,1], Data.snp$data[idx.2,1]) # 0.9999996

id.snp2 = intersect(snp.list, Data.snp$cloneid)
idx.snp2 = match(id.snp2, Data.snp$cloneid)
idx.1 = match(Data.snp$cloneid[idx.snp2], Data$cloneid)
plot(Data$data[idx.1,1], Data.snp$data[idx.snp2,1])   #0.97
  # conclusion: the snp-containing probe sets are unstable but the other probe sets are stable from RMA preprocessing.