## Fecundity Trait Correlations #
## created by jamilla a
## updated: 2020-02-09

library(tidyverse)
library(corrplot)


##functions -----------------------

trait.cor.test = function(data, pheno){
  idx = intersect(rownames(pheno), data$line) #find common lines
  pheno = pheno[rownames(pheno) %in% idx,] #pick only those lines
  data = data[data$line %in% idx,]
  data = data[,-dim(data)[2]] #remove the line column
  
  cor.exp = matrix(NA, dim(pheno)[2], length(as.data.frame(data)))
  pval.exp = matrix(NA, dim(pheno)[2], length(as.data.frame(data)))
  
  for(i in 1:dim(pheno)[2]){
    for(j in 1:dim(as.data.frame(data))[2]){
      cor.exp[i,j] = as.numeric(cor.test(pheno[,i], as.data.frame(data)[,j])$estimate)
      pval.exp[i,j] = cor.test(pheno[,i], as.data.frame(data)[,j])$p.value
    }
  }
  colnames(cor.exp) = colnames(data)
  colnames(pval.exp) = colnames(data)
  rownames(cor.exp) = colnames(pheno)
  rownames(pval.exp) = colnames(pheno)
  return(list('cor'= cor.exp, 'pval'= pval.exp))  
  }


##load in fecundity phenotypes ----------------
setwd("~/de Bivort 2/Fecundity/new analysis/data")
pheno = read_csv('Line Effects Data.csv')
lines = pull(pheno, Line)
pheno = as.matrix(pheno[-1])
rownames(pheno) = lines

df = read_csv('Batch Adjusted Phenotypes.csv')
df = df %>% #take the mean of vials
  group_by(Line) %>% 
  summarize(numF = mean(numF),
            numM = mean(numM),
            Fmass = mean(Fmass),
            Mmass = mean(Mmass))

plot(pheno[,1], df$numF) #check how correlated the line effects are to the mean measured vals
cor.test(pheno[,1], df$numF)


##load in phenotypes from other studies --------------------

setwd("~/de Bivort 2/Fecundity/new analysis/DGRP Data")

#load trait means
#food intake
cafe.female = read.csv('cafe.female.csv', header = F) 
cafe.female$line = str_split_fixed(cafe.female$V1, '_', 2)[,2]
#reorder line in ascending order
cafe.female = cafe.female[order(as.numeric(cafe.female$line)),] 

cafe.male = read.csv('cafe.male.csv', header = F)
cafe.male$line = str_split_fixed(cafe.male$V1, '_', 2)[,2]
cafe.male = cafe.male[order(as.numeric(cafe.male$line)),]

#chill coma
chill.female = read.csv('chillcoma.female.csv', header = F)
chill.female$line = str_split_fixed(chill.female$V1, '_', 2)[,2]
chill.female = chill.female[order(as.numeric(chill.female$line)),] 

chill.male = read.csv('chillcoma.male.csv', header = F)
chill.male$line = str_split_fixed(chill.male$V1, '_', 2)[,2]
chill.male = chill.male[order(as.numeric(chill.male$line)),]

#starvation
starve.female = read.csv('starvation.female.csv', header = F)
starve.female$line = str_split_fixed(starve.female$V1, '_', 2)[,2]
starve.female = starve.female[order(as.numeric(starve.female$line)),] 

starve.male = read.csv('starvation.male.csv', header = F)
starve.male$line = str_split_fixed(starve.male$V1, '_', 2)[,2]
starve.male = starve.male[order(as.numeric(starve.male$line)),]

#fecundity
fec.low = read.csv('Durham_DGRP_LowYeast.csv')
fec.low$line =str_split_fixed(fec.low$X, '_', 2)[,2]
fec.high = read.csv('Durham_DGRP_HighYeast.csv')
fec.high$line =str_split_fixed(fec.high$X, '_', 2)[,2]
fec.bodySize = read.csv('Durham_DGRP_Line_Means_Body_Size.csv')
fec.bodySize$line = fec.bodySize$Line

#nutrients and mean weight
nutrient = read.csv('Uncklessetal2015PLoSGenetics_Immunity_Nutrients.csv')
nutrient = nutrient[-153,] #delete empty row
nutrient$line = nutrient$Line
nutrient = nutrient[order(as.numeric(as.character(nutrient$line))),]

#growth
growth = read.csv('growth.female.male.csv')
growth$line = growth$RAL.ID
growth = growth[order(as.numeric(growth$line)),]

#horvath - dev time and viability
dt_eav = read.csv('Horvath_DevTime_Viability.csv')
dt_eav_LD = dt_eav %>% as_tibble() %>% filter(Density=='LD')
dt_eav_MD = dt_eav %>% as_tibble() %>% filter(Density=='MD')
dt_eav_HD = dt_eav %>% as_tibble() %>% filter(Density=='HD')

##make correlation matrices -----------------
#food intake
cafe.cor.test.fem = trait.cor.test(cafe.female[,-1], pheno)
cafe.cor.test.male = trait.cor.test(cafe.male[,-1], pheno)

#starvation
starve.cor.test.fem = trait.cor.test(starve.female[,-1], pheno)
starve.cor.test.fem$pval.adj = p.adjust(starve.cor.test.fem$pval, method = 'BH')
starve.cor.test.male = trait.cor.test(starve.male[,-1], pheno)

#chill coma
chill.cor.test.fem = trait.cor.test(chill.female[,-1], pheno)
chill.cor.test.male = trait.cor.test(chill.male[,-1], pheno)

#fecundity
fecLOW.cor.test = trait.cor.test(fec.low[,-1], pheno) #low yeast
fecHIGH.cor.test = trait.cor.test(fec.high[,-c(1,8)], pheno) #high yeast

fecBODY.cor.test = trait.cor.test(fec.bodySize[,-1], pheno) #body size
fecBODY.cor.test$pval.adj = p.adjust(fecBODY.cor.test$pval, method = 'BH')
fecBODY.cor.test$pval.adj = matrix(fecBODY.cor.test$pval.adj, 5, 2)

#nutrient
nutr.cor.test = trait.cor.test(nutrient[,-c(1,2)], pheno)
nutr.cor.test$pval.adj = p.adjust(nutr.cor.test$pval, method = 'BH')
nutr.cor.test$pval.adj = matrix(nutr.cor.test$pval.adj, 5, 21)
rownames(nutr.cor.test$pval.adj) = rownames(nutr.cor.test$pval)
colnames(nutr.cor.test$pval.adj) = colnames(nutr.cor.test$pval)

#growth
growth.cor.test = trait.cor.test(growth[,-1], pheno)

#dev time and viability 
dt_eav_LD.cor.test = trait.cor.test(dt_eav_LD[,-1], pheno)
dt_eav_MD.cor.test = trait.cor.test(dt_eav_MD[,-1], pheno)
dt_eav_HD.cor.test = trait.cor.test(dt_eav_HD[,-1], pheno)
dt_eav_HD.cor.test$pval.adj = p.adjust(dt_eav_HD.cor.test$pval, method = 'BH')
dt_eav_HD.cor.test$pval.adj = matrix(dt_eav_HD.cor.test$pval.adj, 5, 4)
#### combine data ---------------------------

total.trait.cor = cbind(food.fem = cafe.cor.test.fem$cor, 
                        food.male = cafe.cor.test.male$cor,
                        starve.fem = starve.cor.test.fem$cor,
                        starve.male = starve.cor.test.male$cor,
                        chill.fem = chill.cor.test.fem$cor,
                        chill.male = chill.cor.test.male$cor,
                        fecLOW = fecLOW.cor.test$cor,
                        fecHIGH = fecHIGH.cor.test$cor,
                        fecBODY = fecBODY.cor.test$cor,
                        nutrient = nutr.cor.test$cor, 
                        dt_eav_LD = dt_eav_LD.cor.test$cor,
                        dt_eav_MD = dt_eav_MD.cor.test$cor,
                        dt_eav_HD = dt_eav_HD.cor.test$cor)

colnames(total.trait.cor)[1:6] = c('food.fem', 'food.male',
                                   'starve.fem', 'starve.male',
                                   'chill.fem', 'chill.male')

total.trait.pval = cbind(food.fem = cafe.cor.test.fem$pval, 
                        food.male = cafe.cor.test.male$pval,
                        starve.fem = starve.cor.test.fem$pval,
                        starve.male = starve.cor.test.male$pval,
                        chill.fem = chill.cor.test.fem$pval,
                        chill.male = chill.cor.test.male$pval,
                        fecLOW = fecLOW.cor.test$pval,
                        fecHIGH = fecHIGH.cor.test$pval,
                        fecBODY = fecBODY.cor.test$pval,
                        nutrient = nutr.cor.test$pval,
                        dt_eav_LD = dt_eav_LD.cor.test$pval,
                        dt_eav_MD = dt_eav_MD.cor.test$pval,
                        dt_eav_HD = dt_eav_HD.cor.test$pval)

colnames(total.trait.pval)[1:6] = c('food.fem', 'food.male',
                                   'starve.fem', 'starve.male',
                                   'chill.fem', 'chill.male')

##multiple test correction
total.trait.pval.adj = p.adjust(total.trait.pval, method = 'BH')
total.trait.pval.adj = matrix(total.trait.pval.adj, dim(total.trait.pval)[1], 
                              dim(total.trait.pval)[2])
colnames(total.trait.pval.adj) = colnames(total.trait.pval)
rownames(total.trait.pval.adj) = rownames(total.trait.pval)

##plotting 
library(corrplot)
corrplot(nutr.cor.test$cor,
         method= "ellipse", 
         tl.cex =0.5, tl.srt = 45, cl.pos = 'r',
         p.mat = nutr.cor.test$pval.adj, sig.level = 0.05, insig = 'blank',
         cl.align.text = "l", cl.cex = 0.4)