Title: Heat mapping for gene: Nat10

Stage 1:

1. Set working directory and load data
2. For ABR change column name Test Date to Assay Date to match other datasets
3. For Plasma chemistry relabel column Na to Sodium to avoid confusion within R as NA means missing data. 
4. ABR and BW reports from the mig database need converting from long to wide data format 
5. Recode the number "Number.Of.Caudal.Vertebrae" to low (<28), normal and high (>29)
6. Calculate an AUC for the growth curve provided the mouse has data from week 4 to week 16
7. Determine the variable list for each screen.  Note currently this is set for data obtained from MGP Select
setwd("Y:\\Natasha\\Gabriel Balmus\\Nat10 heatMap")
ABR=read.csv("Nat10ABRMalesOnly.csv")
names(ABR)[names(ABR)=="Test.Date"]="Assay.Date"
BW_NKarpPull=read.csv("Nat10Bodyweightat4_9_13_16weeks.csv")
Calo=read.csv("Nat10Calorimetry.csv")
PC=read.csv("Nat10CBC_Insulin.csv")
names(PC)[names(PC)=="Na"]="Sodium"
DEXA=read.csv("Nat10DEXA.csv")
Eye=read.csv("Nat10EyeMorphology.csv")
GS=read.csv("Nat10GripStrength.csv")
Haem=read.csv("Nat10Haem.csv")
HW=read.csv("Nat10HeartWeight.csv")
ipGTT=read.csv("Nat10IPGTT.csv")
NAMPA=read.csv("Nat10NaMPADysmorphology.csv")
PBL=read.csv("Nat10PBL.csv")
Xray=read.csv("Nat10XRay.csv")

Xray$Number.Of.Caudal.Vertebrae2=Xray$Number.Of.Caudal.Vertebrae
Xray$Number.Of.Caudal.Vertebrae2[Xray$Number.Of.Caudal.Vertebrae2< 28] <- "low"
Xray$Number.Of.Caudal.Vertebrae2[Xray$Number.Of.Caudal.Vertebrae2>29] <- "high"
Xray$Number.Of.Caudal.Vertebrae2[Xray$Number.Of.Caudal.Vertebrae2>27 & Xray$Number.Of.Caudal.Vertebrae<30] <- "normal"
#convert ABR and BW  to a wide format 
library(tidyr)
ABRwide=spread(data=ABR, key=Frequency, value=Threshold, fill=NA, convert=FALSE)

#Bw has duplicate rows due to being entered twice through having different colony names an issue with BW report - mig looking into this issue
#table(duplicated(BW_NKarpPull))  #36 duplicates #2823 originals
BW_NKarpPullclean=BW_NKarpPull[!duplicated(BW_NKarpPull), ]
#length(BW_NKarpPullclean$Week) # good now 2823
BWwide=spread(data=BW_NKarpPullclean, key=Week, value=Avg.Weight, fill=NA, convert=FALSE)
names(BWwide)[names(BWwide)=="Birth.Date"]="Assay.Date"
names(BWwide)[names(BWwide)=="4"]="Week4"
names(BWwide)[names(BWwide)=="9"]="Week9"
names(BWwide)[names(BWwide)=="13"]="Week13"  #PhenStat falls over if column names are numbers!
names(BWwide)[names(BWwide)=="16"]="Week16"

GrowthCurve=read.csv("GrowthCurveData.csv")
GrowthCurveclean=GrowthCurve[!duplicated(GrowthCurve), ] #remove duplicate entries


weight_outcome<-function(x){
  min_age=min(x$Week, na.rm=TRUE)
    max_age=max(x$Week, na.rm=TRUE) 
    if(min_age>4||max_age<15){
        area="NA"
        return(area)
    }else{
        area= auc(x=x$Week, y=x$Avg.Weight) 

        return(area)            
    }                   
}

#second function splits the data up per mouse and applies the first function and outputs the information into a dataframe. 
Weight_AUC<-function(df){
    df2=df[which(df$Week>3.4), ]
    require(plyr)
    require(MESS)
    aucPerMouse=ddply(.data=df2, .variables=c("Mouse.Name","Genotype", "Gender", "Birth.Date"), .fun=weight_outcome)
    names(aucPerMouse)[names(aucPerMouse)=="V1"]="totalAUC"
    return(aucPerMouse)
}

output=Weight_AUC(GrowthCurveclean)
## Loading required package: plyr
## Loading required package: MESS
## Loading required package: geepack
## 
## Attaching package: 'MESS'
## The following object is masked from 'package:stats':
## 
##     power.t.test
names(output)[names(output)=="Birth.Date"]="Assay.Date"
output$totalAUC=as.numeric(output$totalAUC)
## Warning: NAs introduced by coercion
#sets variable list for each screen
getVariableList<-function(assayName, dataType){
  if(assayName=="DEXA" & dataType=="Continuous"){
        return(VariableList=c("Lean.Mass", "Fat.Mass", "Bone.Mineral.Content", "Nose.To.Tail.Base.Length",  "Bone.Mineral.Density", "Fat.Percentage.Estimate"))

  }else if(assayName=="HAEM" & dataType=="Continuous"){
        return(VariableList=c("Wbc","Rbc","Hgb","Hct","Plt","Mcv","Mch","Mchc","Rdw","Mpv"))        

  }else if(assayName=="PC" & dataType=="Continuous"){
        return(VariableList=c("Sodium", "K", "Cl", "Gluc", "Trigs", "Chol", "Hdl", "Ldl", "Nefac", "Glyc", "Frct", "Amy", "Tp", "Alb", "Ca", "Phos", "Iron", "Ast", "Alp", "Alt", "Creat", "Urea","Tbilc", "Ck", "Mg", "Thyrx", "Insulin")) 

  }else if(assayName=="ABR"& dataType=="Continuous"){
    return(VariableList=c("12","18","24","30","6" , "click"))  

  }else if(assayName=="BW"& dataType=="Continuous"){
   return(VariableList=c("Week4", "Week9", "Week13", "Week16"))  

  }else if(assayName=="GS"&dataType=="Continuous"){
   return(VariableList=c("Grip.Strength.Fore.Paws.1", "Grip.Strength.Fore.Paws.2" ,"Grip.Strength.Fore.Paws.3","Grip.Strength.All.Paws.1","Grip.Strength.All.Paws.2","Grip.Strength.All.Paws.3" ))  

  }else if(assayName=="ipGTT"& dataType=="Continuous"){
   return(variableList=c("Plasma.Glucose.Conc.T0","Plasma.Glucose.Conc.T15","Plasma.Glucose.Conc.T30" ,"Plasma.Glucose.Conc.T60","Plasma.Glucose.Conc.T120", "AUC" )) 

  }else if(assayName=="CALO"& dataType=="Continuous"){
     return(VariableList=c("Total.Av.Vco2", "Total.Av.Vo2", "Total.Av.Activity.Count.Total","Total.Av.Activity.Count" ,"Total.Av.Heat.Animal","Total.Av.Rer", "Cumulative.Food.Intake.Kj","Food.Intake.Rate.G","Total.Water.Intake" ,"Day.Av.Vco2"                  
,"Day.Av.Vo2","Day.Min.V02.Animal","Day.Av.Activity.Count.Total","Day.Av.Activity.Count","Day.Av.Heat.Animal" ,"Day.Av.Rer" ,"Night.Av.Vco2" ,"Night.Av.Vo2","Night.Av.Activity.Count.Total","Night.Av.Activity.Count","Night.Av.Heat.Animal","Night.Av.Rer","Total.Av.Low.Rer","Total.Av.High.Rer" , "Change.Rer"))

  }else if(assayName=="NAMPA"& dataType=="Categorical"){
     return(VariableList=c("Vibrissae.Size"            ,     "Vibrissae.Shape"               
, "Vibrissae.Colour"          ,     "Mouth.Morphology"            ,   "Incisors"                   ,    "Incisors.Morphology"          ,  "Incisor.Colour"               ,  "Pinna"                        ,  "Pinna.Size"                 ,    "Pinna.Morphology"              
, "Pinna.Position"            ,     "Pinna.Rotation"              ,   "External.Genitalia"         ,    "External.Genitalia.Morphology",  "External.Genitalia.Size"      ,  "Forelimb"                    ,   "Forelimb.Size"              ,    "Forelimb.Morphology"           
, "Forepaw"                   ,     "Forepaw.Size"                ,   "Forepaw.Morphology"         ,   "Forepaw.Footpad.Size"          , "Forepaw.Digit.Count.Normal.10" , "Forepaw.Digit.Size"           ,  "Forepaw.Digit.Morphology"     ,  "Forepaw.Digit.Fusion"          
,"Forepaw.Nail.Count"         ,    "Forepaw.Nail.Length"          ,  "Forepaw.Nail.Morphology"      ,  "Forepaw.Nail.Colour"           ,"Hindlimb"                      , "Hindlimb.Size"                 , "Hindlimb.Morphology"          ,  "Hindpaw"                       
, "Hindpaw.Size"              ,     "Hindpaw.Morphology"          ,   "Hindpaw.Footpad.Size"       ,    "Hindpaw.Digit.Count.Normal.10" , "Hindpaw.Digit.Size"          ,   "Hindpaw.Digit.Morphology"    ,   "Hindpaw.Digit.Fusion"       ,    "Hindpaw.Nail.Count"            
, "Hindpaw.Nail.Length"       ,     "Hindpaw.Nail.Morphology"     ,   "Hindpaw.Nail.Colour"        ,    "Tail"                          , "Tail.Length"                 ,   "Tail.Morphology"             ,   "Positional.Passivity"       ,    "Limb.Grasping"                 
,"Trunk.Curl"                 ,    "Transfer.Arousal"             ,  "Gait.Inc.Ataxia"             ,   "Tail.Elevation"                , "Tremor"                       ,  "Startle.Response"             ,  "Corneal.Touch.Reflex"         ,  "Pinna.Touch.Reflex"            
,"Contact.Righting.Reflex"    ,    "Headbobbingcircling"          ,  "Convulsions"                 ,   "Coat.Colouration.General"      , "Head.Coathair.Presence"        , "Head.Coathair.Colour"         ,  "Head.Coathair.Colour.Pattern" ,  "Head.Coathair.Appearance"      
,"Head.Coathair.Length"       ,    "Dorsal.Coathair.Presence"     ,  "Dorsal.Coathair.Colour"      ,   "Dorsal.Coathair.Colour.Pattern", "Dorsal.Coathair.Appearance"    , "Dorsal.Coathair.Length"       ,  "Ventral.Coathair.Presence"    ,  "Ventral.Coathair.Colour"       
,"Ventral.Coathair.Colour.Patter", "Ventral.Coathair.Appearance"   , "Ventral.Coathair.Length"     ,   "Limb.Coathair.Presence"        , "Limb.Coathair.Colour"          , "Limb.Coathair.Colour.Pattern" ,  "Limb.Coathair.Appearance"     ,  "Limb.Coathair.Length"          
,"Paw.Coathair.Presence"       ,   "Paw.Coathair.Colour"           , "Paw.Coathair.Colour.Pattern"  ,  "Paw.Coathair.Appearance"       , "Paw.Coathair.Length"          ,  "Digit.Coathair.Presence"      ,  "Digit.Coathair.Colour"       ,   "Digit.Coathair.Colour.Pattern" 
,"Digit.Coathair.Appearance" ,     "Digit.Coathair.Length"        ,  "Tail.Coathair.Presence"       ,  "Tail.Coathair.Colour"         ,  "Tail.Coathair.Colour.Pattern" ,  "Tail.Coathair.Appearance"     ,  "Tail.Coathair.Length"        ,   "Head.Skin.Appearance"          
,"Head.Skin.Colour"          ,     "Pinna.Skin.Appearance"       ,   "Pinna.Skin.Colour"            ,  "Dorsal.Skin.Appearance"        , "Dorsal.Skin.Colour"           ,  "Ventral.Skin.Appearance"      ,  "Ventral.Skin.Colour"         ,   "Limb.Skin.Appearance"          
,"Limb.Skin.Colour"          ,     "Paw.Footpad.Skin.Appearance" ,   "Paw.Footpad.Skin.Colour"     ,   "Tail.Skin.Appearance"          , "Tail.Skin.Colour"             ,  "Head.Morphology"             ,   "Head.Size"                     
,"Snout.Morphology"          ,     "Snout.Size"                   ,  "Vibrissae"    ))

    }else if(assayName=="XRAY" & dataType=="Categorical"){
   return(VariableList=c( "Caudal.Processes"   , "Sacral.Processes"  ,  "Shape.Of.Spine" ,"Shape.Of.Ribcage" , "Lumbar.Processes"  , "Thoracic.Processes" , "Shape.Of.Vertebrae" ,"Cervical.Processes", "Ulna" ,"Tibia","Teeth" ,"Skull.Shape","Scapula"  , "Shape.Of.Ribs"               
, "Radius" ,"Pelvis" , "Zygomatic.Bone" , "Maxilla", "Mandible" , "Joints" , "Humerus" ,"Fibula" , "Femur"                       
, "Clavicle" ,"Level", "Transitional.Vertebrae"    ,   "Fusion.Of.Processes"  , "Scoliosis" , "Lordosis", "Kyphosis" , "Polysyndactylism" , "Fusion.Of.Vertebrae"         
,"Syndactylism", "Rib.Fusions", "Digit.Integrity" , "Brachydactylism", "Number.Of.Ribs.Right","Number.Of.Ribs.Left" ,"Number.Of.Cervical.Vertebrae", "Number.Of.Lumbar.Vertebrae"   ,"Number.Of.Pelvic.Vertebrae"  , "Number.Of.Thoracic.Vertebra" , "Number.Of.Caudal.Vertebrae2" ))
  }else if (assayName=="HW"& dataType=="Continuous"){
   return(VariableList=c("Heart.Weight"))         
  }else if (assayName=="GrowthCurve"& dataType=="Continuous"){
   return(VariableList=c("totalAUC"))   
  }else if (assayName=="EYE"& dataType=="Categorical"){
   return(VariableList=c("Bulging.Eye", "Eye.Hemorrhage" ,"Eyelid.Morphology" , "Eyelid.Closure" ,"Narrow.Eye.Opening","Cornea", "Corneal.Opacity","Corneal.Vascularisation","Fusion.Between.Cornea.And.Lens","Iris.Pupil","Synechia","Pupil.Position", "Pupil.Shape" , "Pupil.Dilation"  ,"Pupil.Light.Response" , "Iris.Pigmentation" , "Lens" , "Lens.Opacity" , "Retina", "Retinal.Pigmentation", "Retinal.Structure" , "Optic.Disc" , "Retinal.Blood.Vessels" ,"Retinal.Blood.Vessels.Structur" ,"Retinal.Blood.Vessels.Pattern",  "Persistence.Of.Hyaloid.Vascula"))     
  }else if (assayName=="PBL"& dataType=="Continuous"){
   return(variableList=c("Total.T.Cell.Percentage"               , "Alpha.Beta.T.Cell.Percentage","Cd4.Alpha.Beta.T.Cell.Percentage"    ,   "Cd8.Alpha.Beta.T.Cell.Percentage"    ,   "Gamma.Delta.T.Cell.Percentage" ,"Nkt.Cell.Percentage"  , "Nk.Cell.Percentage" ,"Cd4.Cd25.Alpha.Beta.Reg.Percent", "Cd4.Cd44.Cd62l.Alpha.Beta.Percent"   ,   "Cd8.Cd44.Cd62l.Alpha.Beta.Percent"   ,   "Klrg1.Mature.Nk.Cell.Percentage"    ,   "Cd4.Klrg1.Alpha.Beta.T.Cell.Percentage" ,"Cd8.Klrg1.Alpha.Beta.T.Cell.Percentage", "B.Cell.Percentage" ,"Igd.Mature.B.Cell.Percentage"  ,"Monocyte.Percentage" , "Ly6c.Pos.Monocyte.Percentage","Ly6c.Neg.Monocyte.Percentage"  ,"Neutrophil.Percentage"  , "Eosinophil.Percentage"   , "Total.T.Cell.Number" , "Alpha.Beta.T.Cell.Number"  ,"Cd4.Alpha.Beta.T.Cell.Number" ,"Cd8.Alpha.Beta.T.Cell.Number","Gamma.Delta.T.Cell.Number", "Nkt.Cell.Number"  , "Nk.Cell.Number" ,"Cd4.Cd25.Alpha.Beta.Reg.Number"       ,  "Cd4.Cd44.Cd62l.Alpha.Beta.Number"      ,"Cd8.Cd44.Cd62l.Alpha.Beta.Number" ,"Klrg1.Mature.Nk.Cell.Number"        ,    "Cd4.Klrg1.Alpha.Beta.T.Cell.Number","Cd8.Klrg1.Alpha.Beta.T.Cell.Number"   ,  "B.Cell.Number"  , "Igd.Mature.B.Cell.Number" , "Monocyte.Number","Ly6c.Pos.Monocyte.Number" ,"Ly6c.Neg.Monocyte.Number","Neutrophil.Number" , "Eosinophil.Number"   )) 
    }
}   

Stage 2:

1. Run PhenStat for each screen and data type
2. Collect results by analysis methods
3. Adjust for multiple testing and generate a column p_adjusted
4. Assess significance - generate a collumn titled HeatMap and callSummary
5. Write files out

Note: This has been not been wrapped into a single wrapper as which screens you want might depend on the scientist and data availability.

library(PhenStat)

#This wrapper runs PhenStat for a dataset and a variable
#Example usage:  PhenStatWrapper(HW, GenotypeInterest="Nat10/+", VariableInterest="Heart.Weight", dataType="Continuous", analysisMethod="MM")


PhenStatWrapper<-function(df, GenotypeInterest, VariableInterest, dataType, analysisMethod){
  test=PhenList(df, testGenotype=GenotypeInterest, refGenotype="+/+", dataset.colname.sex="Gender", dataset.values.male="Male",
            dataset.values.female="Female", dataset.clean=TRUE, dataset.colname.batch="Assay.Date", dataset.colname.genotype="Genotype", outputMessages=FALSE)
    result <- testDataset(test,depVariable=VariableInterest, method=analysisMethod, dataPointsThreshold=2, transformValues=FALSE, outputMessages=FALSE, equation="withoutWeight")
    output=vectorOutput(result)
    output=c(output, GenotypeInterest)
    return(output)
}

#This wrapper runs PhenStat for multiple varibles for a screen
#Example usage 1: multipleVariablesForDataset(df=HW, GenotypeInterest="Nat10/+", assayName="HW",  dataType="Continuous", analysisMethod="MM")
#Example usage 2: multipleVariablesForDataset(df=ABRwide, GenotypeInterest="Nat10/+", assayName="ABR",  dataType="Continuous", analysisMethod="RR")
#this function is written to return a printed statement if a variable for a dataset fails to fit the model eg "Thrx_Failure to fit PhenStat""

multipleVariablesForDataset<-function(df, GenotypeInterest, assayName, dataType, analysisMethod){

  finalOutput=c()

  variablesToTest=getVariableList(assayName, dataType)

  for (bob in variablesToTest){
    tryCatch(
                {PS_output=PhenStatWrapper(df, GenotypeInterest, VariableInterest=bob, dataType, analysisMethod)
             PS_output=c(PS_output, as.character(bob), as.character(dataType))
           finalOutput=rbind(PS_output, finalOutput)
                    },
                    error=function(e){
            print(paste(bob, "Failure to fit PhenStat", sep="_"))                       
                    }
        )        
  } 
  return(finalOutput) 
}

HWresults=multipleVariablesForDataset(df=HW, GenotypeInterest="Nat10/+", assayName="HW",  dataType="Continuous", analysisMethod="MM")
PCresults=multipleVariablesForDataset(df=PC, GenotypeInterest="Nat10/+", assayName="PC",  dataType="Continuous", analysisMethod="MM")
CALOresults=multipleVariablesForDataset(df=Calo, GenotypeInterest="Nat10/+", assayName="CALO",  dataType="Continuous", analysisMethod="MM")
DEXAresults=multipleVariablesForDataset(df=DEXA, GenotypeInterest="Nat10/+", assayName="DEXA",  dataType="Continuous", analysisMethod="MM")
EYEresults=multipleVariablesForDataset(df=Eye, GenotypeInterest="Nat10/+", assayName="EYE",  dataType="Categorical", analysisMethod="FE")
GSresults=multipleVariablesForDataset(df=GS, GenotypeInterest="Nat10/+", assayName="GS",  dataType="Continuous", analysisMethod="MM")
HAEMresults=multipleVariablesForDataset(df=Haem, GenotypeInterest="Nat10/+", assayName="HAEM",  dataType="Continuous", analysisMethod="MM")
ipGTTresults=multipleVariablesForDataset(df=ipGTT, GenotypeInterest="Nat10/+", assayName="ipGTT",  dataType="Continuous", analysisMethod="MM")
PBLresults=multipleVariablesForDataset(df=PBL, GenotypeInterest="Nat10/+", assayName="PBL",  dataType="Continuous", analysisMethod="MM")
NAMPAresults=multipleVariablesForDataset(df=NAMPA, GenotypeInterest="Nat10/+", assayName="NAMPA",  dataType="Categorical", analysisMethod="FE")
XRAYresults=multipleVariablesForDataset(df=Xray, GenotypeInterest="Nat10/+", assayName="XRAY",  dataType="Categorical", analysisMethod="FE")
BWresults=multipleVariablesForDataset(df=BWwide, GenotypeInterest="Nat10/+", assayName="BW",  dataType="Continuous", analysisMethod="MM")
GrowthCurveresults=multipleVariablesForDataset(df=output, GenotypeInterest="Nat10/+", assayName="GrowthCurve",  dataType="Continuous", analysisMethod="MM")

results=rbind(GrowthCurveresults, BWresults,HWresults,PCresults,  CALOresults, DEXAresults, EYEresults, GSresults, HAEMresults, ipGTTresults,PBLresults, NAMPAresults, XRAYresults )
write.csv(results, "temp.csv")
results=read.csv("temp.csv")
ContData=subset(results, results$X.3 == "Continuous")
ContData$p_adjusted=p.adjust(ContData$Genotype.contribution, method="fdr")
print(c("Number rows significant in Continuous dataset:", nrow(ContData[ContData$p_adjusted<=0.05,])))
## [1] "Number rows significant in Continuous dataset:"
## [2] "3"
print(c("Number rows in Continuous dataset:", nrow(ContData))) 
## [1] "Number rows in Continuous dataset:"
## [2] "126"
ContData$callSummary=ifelse(ContData$p_adjusted<0.05, "significant", "not significant")
ContData$HeatMap=ifelse(ContData$p_adjusted<0.05, 1, 0)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:plyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
DataTovisualise= select(ContData, Dependent.variable, HeatMap, p_adjusted,Genotype.percentage.change )
filter(DataTovisualise, HeatMap==1)
##            Dependent.variable HeatMap  p_adjusted
## 1           Neutrophil.Number       1 0.032577763
## 2 Klrg1.Mature.Nk.Cell.Number       1 0.002271607
## 3              Nk.Cell.Number       1 0.002271607
##       Genotype.percentage.change
## 1 Female: -31.85%, Male: 162.21%
## 2     Female: 9.2%, Male: 64.01%
## 3     Female: 7.99%, Male: 57.6%
CatData=subset(results, results$X.3 == "Categorical")
CatData$MinPvalue <- with(CatData, pmin(Sex.MvKO.p.val, Sex.FvKO.p.val, na.rm = TRUE))
CatData$p_adjusted=p.adjust(CatData$MinPvalue, method="fdr")
print(c("Number rows significant in Categorical dataset:", nrow(CatData[CatData$p_adjusted<=0.05,]))) 
## [1] "Number rows significant in Categorical dataset:"
## [2] "0"
print(c("Number rows  in Categorical dataset:", nrow(CatData))) 
## [1] "Number rows  in Categorical dataset:"
## [2] "185"
CatData$callSummary=ifelse(CatData$p_adjusted<0.05, "significant", "not significant")
CatData$HeatMap=ifelse(CatData$p_adjusted<0.05, 1, 0)
#problem  ABR - low n 4 - males only and lack of variation for most frequency means the MM won't fit
ABRresults=multipleVariablesForDataset(df=ABRwide, GenotypeInterest="Nat10/+", assayName="ABR",  dataType="Continuous", analysisMethod="RR")
write.csv(ABRresults, "temp.csv")
ABRresults=read.csv("temp.csv")
library(tidyr)
ABRresults=separate(ABRresults, col=Genotype.p.Val, into=c("Genotype_pValue_lowTest", "Genotype_pValue_highTest"), sep="\\,")
ABRresults$MinPvalue <- with(ABRresults, pmin(Genotype_pValue_lowTest, Genotype_pValue_highTest, na.rm = TRUE))
ABRresults$p_adjusted=p.adjust(ABRresults$MinPvalue, method="fdr")
print(c("Number rows significant in ABR dataset:", nrow(ABRresults[ABRresults$p_adjusted<=0.05,])))
## [1] "Number rows significant in ABR dataset:"
## [2] "0"
print(c("Number rows  in ABR dataset:", nrow(ABRresults)))
## [1] "Number rows  in ABR dataset:" "6"
ABRresults$callSummary=ifelse(ABRresults$p_adjusted<0.05, "significant", "not significant")
ABRresults$HeatMap=ifelse(ABRresults$p_adjusted<0.05, 1, 0)
write.csv(ABRresults, "ABRresults_Nat10Het.csv")
write.csv(ContData, "contResults_Nat10Het.csv")
write.csv(CatData, "catResults_Nat10Het.csv")
SummaryData=rbind(select(CatData,Dependent.variable, HeatMap), select(ContData,Dependent.variable, HeatMap), select(ABRresults,Dependent.variable, HeatMap))