This markdown recreates Figures 2B and 2D from Gillard, Grønvold et al. 2020

# plot symmetry correlations


# load the expression data
combExprMat <- readRDS("data/BSNormalize/combExprMat.RDS")
# load EVE results
EVEresTbl <- readRDS("data/runEVE/EVEresTbl.RDS")
# Load dN/dS results
dndsRes <- readRDS("input-data/dNdS_aBSREL/dndsRes.RDS")
# salmonid colums in expression table
salCols = colnames(combExprMat) %in% c("Ssal","Salp","Omyk","Okis")
dupTbl <- 
  EVEresTbl %>%
  filter(!N3partial,dup_type=="duplicate") %>% 
  # xSal = mean expression of salmonids across N9 clade
  mutate( xSal = rowMeans(combExprMat[N9,salCols])) %>% 
  # add dnds value per clade
  left_join( select(dndsRes, N9, dnds=baseline_omega), by="N9") %>% 
  # place data for each duplicate in different columns
  group_by(N3) %>% 
  arrange(pval) %>% # sort by pval
  summarise(sigDir=paste0(sigDir[1], sigDir[2]), N9_A = N9[1], N9_B = N9[2], id_A = Ssal[1], id_B = Ssal[2],
             xSal_A = xSal[1], xSal_B = xSal[2],
             dnds_A = dnds[1], dnds_B = dnds[2]) %>% 
  # xPike = mean pike (outgroup) expression
  mutate( xPike = rowMeans(combExprMat[N9_A,colnames(combExprMat)=="Eluc"]))

Figure 2A: association between ohnolog regulatory categories and ohnolog expression level evolution symmetry

dupTbl %>% 
  mutate(sigDir_abs = gsub('updown|downup', 'both', sigDir)) %>%
  filter(sigDir_abs != '') %>%
  ggplot(aes(y=sigDir_abs,x=abs(xSal_A-xSal_B), fill=sigDir_abs)) + 
  geom_violin() +
  theme_bw() +
  theme(plot.margin=unit(c(1,1,1.5,1.2),"cm")) +
  stat_summary(fun=median, geom="point", shape=23, size=2) +
  scale_fill_manual(values = c(both="#993094", down="#298ABD", downdown="#1B4F8B",
                               up="#F62D31", upup="#8B1B1A")) +
  theme(axis.text.x = element_text(angle = 0, hjust = 1, size=rel(1.8),color="black")) +
  theme(axis.text.y = element_text(angle = 0, hjust = 1, size=rel(1.8),color="black")) +
  labs(x = "abs(ohnolog expression level difference)", y = "") + 
  theme(axis.title.x = element_text(size = rel(1.8), angle = 00, vjust = -1)) 

# wilcox pvalues:
upup <- dupTbl %>% 
  mutate(sigDir_abs = gsub('updown|downup', 'both', sigDir)) %>%
  filter(sigDir_abs == 'upup') %>%  mutate( diff = abs(xSal_A-xSal_B)) %>% select(diff)
  
up <- dupTbl %>% 
  mutate(sigDir_abs = gsub('updown|downup', 'both', sigDir)) %>%
  filter(sigDir_abs == 'up') %>%  mutate( diff = abs(xSal_A-xSal_B)) %>% select(diff)

downdown <- dupTbl %>% 
  mutate(sigDir_abs = gsub('updown|downup', 'both', sigDir)) %>%
  filter(sigDir_abs == 'downdown') %>%  mutate( diff = abs(xSal_A-xSal_B)) %>% select(diff)

down <- dupTbl %>% 
  mutate(sigDir_abs = gsub('updown|downup', 'both', sigDir)) %>%
  filter(sigDir_abs == 'down') %>%  mutate( diff = abs(xSal_A-xSal_B)) %>% select(diff)

upup.up <- wilcox.test(upup$diff, up$diff, alternative = 'less')
downdown.down <- wilcox.test(downdown$diff, down$diff, alternative = 'less')

stat.tests <- data.frame(test = c('both up VS one up', 'both down VS one down'),
                         p_value = c(upup.up$p.value, downdown.down$p.value))

# dev.copy(pdf, 'symmetry_expression_evolutionary_categories.pdf')
# dev.off()


library(knitr)
kable(stat.tests, caption='Wilcox test of differences in expression symmetry')
Wilcox test of differences in expression symmetry
test p_value
both up VS one up 0
both down VS one down 0




###Figure 2D: correlation between expression and sequence evolution symmetry

dupTbl <- 
  EVEresTbl %>%
  filter(!N3partial,dup_type=="duplicate") %>% 
  # xSal = mean expression of salmonids across N9 clade
  mutate( xSal = rowMeans(combExprMat[N9,salCols])) %>% 
  # add dnds value per clade
  left_join( select(dndsRes, N9, dnds=baseline_omega), by="N9") %>% 
  # place data for each duplicate in different columns
  group_by(N3, Eluc) %>% 
  arrange(pval) %>% # sort by pval
  summarise(sigDir=paste0(sigDir[1], sigDir[2]), N9_A = N9[1], N9_B = N9[2],
             xSal_A = xSal[1], xSal_B = xSal[2],
             dnds_A = dnds[1], dnds_B = dnds[2]) %>% 
  ungroup() %>% 
  # xPike = mean pike (outgroup) expression
  mutate( xPike = rowMeans(combExprMat[N9_A,colnames(combExprMat)=="Eluc"])) %>% 
  # add dN/dS for the pike branch
  left_join( by="Eluc", filter(dndsRes, grepl("Eluc",node)) %>% transmute(Eluc=sub("_Eluc","",node), dnds_Eluc = baseline_omega) )


pvalue <- dupTbl %>% drop_na(dnds_Eluc) %>% filter(dnds_Eluc < 1, sigDir != '' ) %>%
  mutate(y=dnds_Eluc, x=abs(xSal_A-xSal_B)) %>% select(y, x) %>% 
  do(tidy(cor.test(.$x, .$y))) %>% select(p.value)

dupTbl %>% 
  # filter extreme dnds values
drop_na(dnds_Eluc) %>% filter(dnds_Eluc < 1,  sigDir != '' ) %>%
  ggplot(aes(y=dnds_Eluc,x=abs(xSal_A-xSal_B))) + geom_point() + geom_smooth(method="lm",formula = y ~ x) +
  
  annotate(geom="label", x=12, y=0.35, 
           label=paste('P-value =', round(pvalue$p.value, 8)),
           color="darkblue", size=rel(6)) +
  
  theme(axis.text.x = element_text(angle = 0, hjust = 1, size=rel(1.8),color="black")) +
  theme(axis.text.y = element_text(angle = 0, hjust = 1, size=rel(1.8),color="black")) +
  labs(x = "ohnolog pair expression asymmetry", y = "Pike ortholog dN/dS") + 
  theme(axis.title.x = element_text(size = rel(1.8), angle = 00, vjust = -0.5)) +
  theme(axis.title.y = element_text(size = rel(1.8), angle = 90, vjust = -0.5)) 

# dev.copy(pdf, 'Pike_dNdS_vs_expression_evolution_asymmetry.pdf')
# dev.off()