rm(list=ls())

####################################################
#                    1. Load data                 #
####################################################

###packages###
library(readr)
library(tidyr)
library(dplyr) 
library(broom)
library(ggplot2)
library(ggExtra)
library(maps)
library(RColorBrewer)
library(mosaic)
library(tidyverse)
library(modelr) #required
library(performance) #required, model stats
library(see) #required to visualise model plots with performance package
library(qqplotr) #required for performance package
library(MASS) #boxcox
library(olsrr) #boxcox
library(lmtest) #boxcox
library(trafo) #yeo-johnson
library(car) #qqPlot
library(coin) #mood median test
library(rcompanion) #posthoc median test
library(ggpubr) #for significance lines
library(hexbin) #for density colours
library(cowplot) #for combined plot
library(ggVennDiagram) #for venn diagram
library(data.table) #for faster reading
library(viridis) #for colour scale

###load data###
#male data
msnps<-fread("./Data/male_maf0.05.assoc", header = T)
#female data
fsnps<-fread("./Data/female_maf0.05.assoc", header = T)
#allele count data
counts<-fread("./Data/f3c.lhm.snp.frq.edited.txt", header = T)
#smaller dataset with gene names
extendedsnps<- fread("./Data/Candidate_SNPs_extended_information.txt", header = T)
#Tajima's D measure of balancing selection
balancesnps<- fread("./Data/Assoc_Data_with_TajD_v2.txt", header = T)
balancesnps<- cbind(balancesnps, counts)
filteredsnps<- subset(balancesnps, balancesnps$fdr<= 0.3) #fdr used in study
rm(balancesnps) #remove large dataset from memory
###manipulate data###
#combine into one dataset
filteredsnps <- cbind(filteredsnps, extendedsnps)
snps<- msnps
snps$M_Effect<- snps$Effect #rename effect column
snps$Effect<- NULL #remove old column
snps$F_Effect<- fsnps$Effect
#NOTE- not all are SA
rm(msnps)
rm(fsnps)
#classify SA
snps$SA <- NA
snps$SA[snps$M_Effect*snps$F_Effect < 0] <- 'Y'
snps$SA[snps$M_Effect*snps$F_Effect >= 0] <- 'N'
#fwrite(snps, file = "all_snps.csv", sep = ",")
sasnps<- subset(snps, snps$SA == "Y")
#absolute effect sizes
sasnps$M.A.effect <- gsub("\\-", "", as.character(sasnps$M_Effect))
sasnps$F.A.effect <- gsub("\\-", "", as.character(sasnps$F_Effect))
#combining effect sizes
sasnps$combined_effect <- as.numeric(sasnps$M.A.effect) + as.numeric(sasnps$F.A.effect)
#ratio of effects
sasnps$effratio<-NA
sasnps$effratio[sasnps$M_Effect < sasnps$F_Effect]<- sasnps$F_Effect[sasnps$M_Effect < sasnps$F_Effect]/sasnps$M_Effect[sasnps$M_Effect < sasnps$F_Effect]
sasnps$effratio[sasnps$F_Effect < sasnps$M_Effect]<- sasnps$M_Effect[sasnps$F_Effect < sasnps$M_Effect]/sasnps$F_Effect[sasnps$F_Effect < sasnps$M_Effect]

# To output list of SA snps
#fwrite(sasnps, file = "sa_snps_effratio.csv", sep = ",")

#Combine with smaller dataset
filteredsnps<- subset(filteredsnps, select=which(!duplicated(names(filteredsnps)))) 
mergesnps<- merge(filteredsnps, sasnps, by.x = "Predictor", by.y = "Predictor")
###Linkage Disequilibrium data###
#load data
LD_names<- c("Region_start","Region_end","D")
LD2L<- fread("./Data/LD files/RAL_D_1kb-2L.bedgraph", skip = 1)
LD2L$V1<- NULL
names(LD2L)<- LD_names
LD2R<- fread("./Data/LD files/RAL_D_1kb-2R.bedgraph", skip = 1)
LD2R$V1<- NULL
names(LD2R)<- LD_names
LD3L<- fread("./Data/LD files/RAL_D_1kb-3L.bedgraph", skip = 1)
LD3L$V1<- NULL
names(LD3L)<- LD_names
LD3R<- fread("./Data/LD files/RAL_D_1kb-3R.bedgraph", skip = 1)
LD3R$V1<- NULL
names(LD3R)<- LD_names
LDX<- fread("./Data/LD files/RAL_D_1kb-X.bedgraph", skip = 1)
LDX$V1<- NULL
names(LDX)<- LD_names
#assign LD function
assign_LD<- function(chrompos,ldtable){
  ldindex<- ceiling(chrompos/1000)
  return(ldtable$D[ldindex]) 
}
mergesnps$LD1k<- 0
#subset by chromosome
snps2R<- subset(mergesnps, mergesnps$Chrom.r5=='2R')
snps2L<- subset(mergesnps, mergesnps$Chrom.r5=='2L')
snps3R<- subset(mergesnps, mergesnps$Chrom.r5=='3R')
snps3L<- subset(mergesnps, mergesnps$Chrom.r5=='3L')
snpsX<- subset(mergesnps, mergesnps$Chrom.r5=='X')
#apply function to match LD windows
for (i in 1:nrow(snps2R)){
  snps2R$LD1k[i]<- assign_LD(snps2R$Basepair.x[i], LD2R)
}
for (i in 1:nrow(snps2L)){
  snps2L$LD1k[i]<- assign_LD(snps2L$Basepair.x[i], LD2L)
}
for (i in 1:nrow(snps3R)){
  snps3R$LD1k[i]<- assign_LD(snps3R$Basepair.x[i], LD3R)
}
for (i in 1:nrow(snps3L)){
  snps3L$LD1k[i]<- assign_LD(snps3L$Basepair.x[i], LD3L)
}
for (i in 1:nrow(snpsX)){
  snpsX$LD1k[i]<- assign_LD(snpsX$Basepair.x[i], LDX)
}
#merge tables together
mergesnps<- rbind(snps2R, snps2L, snps3R, snps3L, snpsX)
#Variance of effect sizes
M_mean<- mean(mergesnps$M_Effect)
F_mean<- mean(mergesnps$F_Effect)
effratio_mean<- mean(mergesnps$effratio)
mergesnps$effratio_variance<- NA
for(i in 1:nrow(mergesnps)){
  mergesnps$effratio_variance[i]<- (mergesnps$effratio[i] - effratio_mean)^2
}
#group functional categories (aka consequences)
mergesnps$con_class<-NA
for (i in 1:nrow(mergesnps)){
  if(mergesnps$Consequence[i] == 'missense_variant'){
#  if(mergesnps$Consequence[i] == 'missense_variant' 
#     | mergesnps$Consequence[i] == 'splice_region_variant' | mergesnps$Consequence[i] == 'splice_acceptor_variant' | mergesnps$Consequence[i] == 'splice_donor_variant'){
    mergesnps$con_class[i] <- 'missense_variant'
  }else if(mergesnps$Consequence[i] == 'synonymous_variant'){
    mergesnps$con_class[i] <- 'synonymous_variant'
#  }else if(mergesnps$Consequence[i] == 'splice_region_variant' | mergesnps$Consequence[i] == 'splice_acceptor_variant' | mergesnps$Consequence[i] == 'splice_donor_variant'){
#    mergesnps$con_class[i] <- 'splice_variant'
  }else{
    mergesnps$con_class[i] <- 'regulatory_variant'
  }
}

#Output
fwrite(mergesnps, file = "./Output/sa_snps_effratio_filtered.txt", sep = "\t")
fwrite(sasnps, file = "./Output/sa_snps_effratio.txt", sep = "\t")

####################################################
#        2. Effect size ratio analysis             #
####################################################

#Load datasets from previous section
sasnps<-fread("./Output/sa_snps_effratio.txt", header = T, sep = "\t", stringsAsFactors = T)
mergesnps<-fread("./Output/sa_snps_effratio_filtered.txt", header = T, sep = "\t", stringsAsFactors = T)
##creating locus string for database input
mergesnps$locus<- paste(mergesnps$Chrom.r5, ":", mergesnps$Pos.r5, "..", mergesnps$Pos.r5, sep = "")
###Adding in major alleles (A2 in original dataset)
#minor dataset and allele counts
mergesnps$total_alleles<- mergesnps$N_CHR/2
#record prevalence as minor alleles
mergesnps$prevalence <- "minor"
mergesnps$EAF <- mergesnps$MAF.x
mergesnps$count_A <- (as.numeric(sub('..', '', mergesnps$`{ALLELE:COUNTB}`)))/2 #minor alleles listed as B
mergesnps$count_B <- (as.numeric(sub('..', '', mergesnps$`{ALLELE:COUNTA}`)))/2

#combine into one
mergesnps <- mergesnps
mergesnps$prevalence<- as.factor(mergesnps$prevalence)
mergesnps$Consequence<- as.factor(mergesnps$Consequence)

#proportion of allele A- should match EAF
mergesnps$count_over_sample <- mergesnps$count_A/mergesnps$total_alleles
#convert drosophila effect sizes to d
#converting z-score to effect size
#where z is z-score, n is sample size
ztod <- function(z, n){
  if(z<=0){ #account for negative z-score
    z <- -1*z
    m = -1 #to convert to negative if input negative
  }else{
    m = 1
  }
  a = z*sqrt(n)
  b = 1-(sqrt((z^2)*(n^-1)))
  c = sqrt(a/b)
  d = m*c
  return(d)
}
ztod <- function(z, n){
  m <- ifelse(z<= 0, -1, 1) #account for negative z-score
  z = m*z #convert negative to positive for calculation
  a = z*sqrt(n)
  b = 1-(sqrt((z^2)*(n^-1)))
  c = sqrt(a/b)
  d = m*c #convert negative back if applicable
  return(d)
}

#Calling function converting z to d
mergesnps$male_d<- ztod(mergesnps$M_Effect, mergesnps$total_alleles)
mergesnps$female_d<- ztod(mergesnps$F_Effect, mergesnps$total_alleles)

#SA index calculation
SAindex<- function(Male,Female){
  V1<- Male*Female
  V2<- Male^2 + Female^2
  V3<- sqrt(V2/2)
  V4 <- V1/V3
  return(V4)
}
mergesnps$SAI<- SAindex(mergesnps$male_d,mergesnps$female_d)
mergesnps$SAI_w<- SAindex(mergesnps$M_Effect,mergesnps$F_Effect) #SAI of dataset effect size
#Product of effects
mergesnps$Effect_product<- mergesnps$M_Effect*mergesnps$F_Effect
#Filtering
missense_snps<- subset(mergesnps, mergesnps$Consequence == "missense_variant")

###Statistics and modelling###
#main model 
main_mod<- glm(cbind(count_A, count_B)~effratio, data = mergesnps, family = quasibinomial)
summary(main_mod)
anova(main_mod, test = 'Chi')
m_pointsize <- log(1/missense_snps$effratio_variance)
#make predictions for plotting
main_predict<-data.frame(effratio = mergesnps$effratio, EAF = predict(main_mod, newdata = mergesnps, type = 'response'))
#add in consequence
c_model <- glm(cbind(count_A, count_B)~effratio + con_class, data = mergesnps, family = quasibinomial)
summary(c_model)
anova(c_model, test = 'Chi')
#include interaction terms
c_model2 <- glm(cbind(count_A, count_B)~effratio + con_class + effratio:con_class, data = mergesnps, family = quasibinomial)
summary(c_model2)
anova(c_model2, test = 'Chi')

##Models by consequence class
glin_mod<- function(df){
  glm(cbind(count_A, count_B)~effratio, data = df, family = quasibinomial)
}
glin_mod_sum<- function(df){
  summary(glm(cbind(count_A, count_B)~effratio, data = df, family = quasibinomial))
}
anova_glin_mod<- function(df){
  anova(glm(cbind(count_A, count_B)~effratio, data = df, family = quasibinomial), test = 'Chi')
}
glin_mods<- by(mergesnps, mergesnps$con_class, glin_mod)
glin_mod_summaries<-by(mergesnps, mergesnps$con_class, glin_mod_sum)
glin_mod_summaries
glin_mod_anovas<-by(mergesnps, mergesnps$con_class, anova_glin_mod)
glin_mod_anovas
#significant for all three classes

#assign model results to table
prevalences<- unique(mergesnps$prevalence)
consequences<- unique(mergesnps$con_class)
#filter out consequences with low numbers of alleles
consequence_frame<- as.data.frame(table(mergesnps$con_class))
consequence_frame<- subset(consequence_frame, consequence_frame$Freq>=4)
consequences<- subset(consequences, consequences %in% consequence_frame$Var1)
#blank dataframe
consequence_models = NULL
#make temporary dataframe
for(j in 1:length(consequences)){ #per consequence
  t2_alleles<- subset(mergesnps, mergesnps$con_class == consequences[j])
  t_model_an<- anova(glin_mod(t2_alleles), test = 'Chi')
  t_model_sum<- summary(glin_mod(t2_alleles))
  #combine into dataframe
  consequence_models<- rbind(consequence_models, data.frame(con_class = as.character(consequences[j]), Slope = t_model_sum$coefficients[2],
                                                            Deviance = t_model_an$Deviance[2], P = t_model_an$`Pr(>Chi)`[2], SE = t_model_an$`Pr(>Chi)`[2]))
}

##Make predictions
#use only significant models
significant_models<- subset(consequence_models, consequence_models$P<= 0.05)
#blank dataframe for predictions
count_all_predicted <- data.frame(matrix(vector(),0, 4,
                                         dimnames = list(c(), c("count_over_sample_predicted", "effratio", "con_class", "prevalence"))),
                                  stringsAsFactors = T)
for(i in 1:nrow(significant_models)){
  t2_alleles<- subset(mergesnps, mergesnps$con_class == significant_models$con_class[i])
  t_model<- glin_mod(t2_alleles)
  #add predictions to dataframe
  count_all_predicted<- rbind(count_all_predicted, data.frame(count_over_sample_predicted = predict(t_model, t2_alleles, type = "response"), 
                                                              effratio = t2_alleles$effratio, 
                                                              con_class = t2_alleles$con_class, 
                                                              prevalence = t2_alleles$prevalence))
}                                        
##Preparation for plotting
#transforming effect size ratio for plotting
mergesnps$peffratio<- -log10(-1*mergesnps$effratio)
#labels for graphs
consequence_labels<-c('missense_variant' = 'Missense', 'synonymous_variant' = 'Synonymous variant', 'regulatory_variant' = 'Regulatory variant')
pmergesnps<- subset(mergesnps, (mergesnps$con_class %in% consequence_frame$Var1))

####################################################
#          3.  Recombination analysis              #
####################################################

#####Functions#####
locationwriter<- function(chromosome, loci){#function for merging co-ordinates into chr:start..end form
  chromlocation<- paste(chromosome, ":", loci, sep = "")
  return (chromlocation)
}

locationsplitstart<- function(loci){#getting start point from chr:start..end form
  trimstart<- gsub(".*:","",loci)
  trim<- gsub("\\..*", "", trimstart)
  return(as.numeric(trim))
}

locationsplitend<- function(loci){#getting end point from chr:start..end form
  trim<- gsub("^.*\\.","", loci)
  return(as.numeric(trim))
}

locationchromosomeextract<- function(loci){#gets the chromosome out of long form locus position
  trim<- gsub(":.*", "", loci)
  return(trim)
}

recomlookup<- function(chromosome, locusstart, locusend, recombination){#gets midpoint Comeron rate from table, requires for loop for datasets
  locusmidpoint <- (locusstart+locusend)/2 #calculates midpoint
  pos = which(chromosome == recombination$`Chromosome Arm` & locusmidpoint >= recombination$`Window begins` & locusmidpoint <= recombination$`Window ends`) #finds correct row in recombination dataset
  if(length(pos)>0){
    if(is.na(pos)){
      rrate<- NA
    }else{
      rrate<- recombination$`CO rate (cM/Mb/female meiosis)`[pos] #extracts recombination rate
    }
  }else{
    rrate<- NA
  }
  return(rrate)
}


#calculate mode
calcmode <- function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}

###Load data###
##import Comeron dataset
#NOTE: must use release 5 co-ordinates
comeron<- fread("./Data/Comeron.2012.10.15.txt", header = T, sep = "\t") #Release 5 co-ordinates
##import dataset for all genes
all_genes<- fread("./Data/flybase/gene_map_table_fb_2021_02.txt", sep = "\t", header = T, stringsAsFactors = T)
#subset for melanogaster genes
dmel_genes<- subset(all_genes, all_genes$organism_abbreviation == "Dmel")
#remove all genes to save memory
rm(all_genes)
dmel_genes$positiontrim<- gsub("\\(.*", "", dmel_genes$sequence_loc)

#read in converted co-ordinates
all_coordinates<- fread("./Data/flybase/Flybase_Converted_Coordinates_all_melanogaster_genes.tsv", sep = "\t", stringsAsFactors = T, header = F)
names(all_coordinates)<- c("R6", "R5", "notes")
dmel_genes<- merge(dmel_genes, all_coordinates, by.x = "positiontrim", by.y = "R6")
#get r5 coordinates
dmel_genes$chrom.r5<- locationchromosomeextract(dmel_genes$R5)
dmel_genes$start.r5<- locationsplitstart(dmel_genes$R5)
dmel_genes$end.r5<- locationsplitend(dmel_genes$R5)
#remove error values
dmel_genes_correct<- subset(dmel_genes, dmel_genes$chrom.r5 %in% comeron$`Chromosome Arm`)
dmel_genes_correct<- subset(dmel_genes_correct, !is.na(dmel_genes_correct$start.r5))
#get recombination rates- long step
#set up progress bar
pb = txtProgressBar(min = 0, max = nrow(dmel_genes_correct), initial = 0, style = 3) 
for(i in 1:nrow(dmel_genes_correct)){ #Time-consuming step!
  dmel_genes_correct$recombination_rate[i]<- recomlookup(dmel_genes_correct$chrom.r5[i], dmel_genes_correct$start.r5[i], dmel_genes_correct$end.r5[i], comeron)
  setTxtProgressBar(pb,i)
  close(pb)
}
##SA genes
#make dataset of only female beneficial alleles
mergesnps$comeron<- NA
for(i in 1:nrow(mergesnps)){
  mergesnps$comeron[i]<- recomlookup(mergesnps$Chrom.r5[i], mergesnps$Pos.r5[i], mergesnps$Pos.r5[i], comeron)
}

##Calculate means for each gene
sa_gene_Comeron_recombination<- setNames(aggregate(mergesnps$comeron~mergesnps$Gene, FUN = mean), c("Gene", "Mean.Recombination.Comeron"))
sa_gene_SAI<- setNames(aggregate(mergesnps$SAI~mergesnps$Gene, FUN = mean), c("Gene", "SAI"))
sa_gene_Comeron_recombination<- sa_gene_Comeron_recombination[order(sa_gene_Comeron_recombination$Gene),]
sa_gene_SAI<- sa_gene_SAI[order(sa_gene_SAI$Gene),]
sa_genes<- merge(sa_gene_Comeron_recombination, sa_gene_SAI)

###Human comparison###
humansnps<-read.table("./Data/human_snps_recombination.txt", header = T, sep = "\t")
humansnps$SAI<- SAindex(humansnps$Male_d, humansnps$Female_d)

#make table of only directly comparable values
humancomp<- data.frame(Gene = humansnps$Gene, 
                       SAI = humansnps$SAI, 
                       Recombination = humansnps$deCODE_sex_average,
                       con_class = humansnps$Trait_class)
humancomp$Species<- "Human"
#average SAI
human_SAI_average<- setNames(aggregate(humancomp$SAI~humancomp$Gene, FUN = mean), c("Gene", "SAI"))
humancomp_unique<- humancomp[!duplicated(humancomp$Gene),]
humancomp_unique$SAI<- human_SAI_average$SAI
#same as above for drosophila, missense genes only
missense_genes<- subset(sa_genes, sa_genes$Gene %in% subset(mergesnps, mergesnps$Consequence == "missense_variant")$Gene)

flycomp<- data.frame(Gene = sa_genes$Gene,
                     SAI = sa_genes$SAI,
                     Recombination = sa_genes$Mean.Recombination.Comeron)
flycomp$con_class = "Dros"
flycomp$Species <- "Drosophila"
#flycomp for missense variants
mis_snps<- subset(mergesnps, mergesnps$con_class == 'missense_variant')
#mean SAI by gene
flycomp<- setNames(aggregate(mergesnps$SAI~mergesnps$Gene, FUN = mean), c("Gene", "SAI"))
#get recombination rates
flycomp<- merge(flycomp, sa_genes[, c('Gene', 'Mean.Recombination.Comeron')], by = 'Gene', all.x = T)
colnames(flycomp)[which(names(flycomp) == "Mean.Recombination.Comeron")] <- "Recombination"
flycomp$con_class = "Dros"
flycomp$Species <- "Drosophila"

#combine
humanfly<- rbind(humancomp_unique, flycomp)

#modelling
speciesmod<- glm(data = humanfly, Recombination~SAI + Species + Species:SAI)
summary(speciesmod)

humanmod<- glm(data = humancomp_unique, Recombination~SAI, family = quasipoisson(link = "log"))
summary(humanmod)
anova(humanmod)

flymod<- glm(data = flycomp, Recombination~SAI)
summary(flymod)
anova(flymod)

#split SAI into low/medium/high categories
flycomp$SAI_cat<- NA
for(i in 1:nrow(flycomp)){
  if(flycomp$SAI[i]<= quantile(flycomp$SAI, .33)){
    flycomp$SAI_cat[i]<- "low" 
  }else if(flycomp$SAI[i]>=quantile(flycomp$SAI, .67)){
    flycomp$SAI_cat[i]<- "high"
  }else{
    flycomp$SAI_cat[i]<- "medium"
  }
}

humancomp$SAI_cat<- NA
for(i in 1:nrow(humancomp_unique)){
  if(humancomp_unique$SAI[i]<= quantile(humancomp_unique$SAI, .33)){
    humancomp_unique$SAI_cat[i]<- "low" 
  }else if(humancomp_unique$SAI[i]>=quantile(humancomp_unique$SAI, .67)){
    humancomp_unique$SAI_cat[i]<- "high"
  }else{
    humancomp_unique$SAI_cat[i]<- "medium"
  }
}

leveneTest(Recombination~SAI_cat, data = flycomp)
summary(aov(Recombination~SAI_cat, data = flycomp))

####Median testing####

t_flycomp<- subset(flycomp, flycomp$SAI>-1.75)
#median testing
median_test(Recombination~as.factor(SAI_cat), data = flycomp) #significant
median_test(Recombination~as.factor(SAI_cat), data = humancomp_unique) #non significant
leveneTest(Recombination~as.factor(SAI_cat), data = flycomp) #in flies, non significant
leveneTest(Recombination~as.factor(SAI_cat), data = humancomp_unique)
pairwiseMedianTest(Recombination~as.factor(SAI_cat), data = flycomp)
#prepare for plotting
rectest<- pairwiseMedianTest(Recombination~as.factor(SAI_cat), data = flycomp)
recplot<-data.frame(group1 = c('high', 'high', 'low'), group2 = c('low', 'medium', 'medium'), p = rectest$p.value,
                    p.adjust = rectest$p.adjust, y.position = c(7, 8, 6), Species = c('Drosophila', 'Drosophila', 'Drosophila'))
recplot$p.adjust.sig<- signif(recplot$p.adjust, 3)

#pairwise testing
#low v medium
median_test(Recombination~as.factor(SAI_cat), data = subset(flycomp, !(flycomp$SAI_cat == 'high')))
lvm <- median_test(Recombination~as.factor(SAI_cat), data = subset(flycomp, !(flycomp$SAI_cat == 'high')))
#low v high
median_test(Recombination~as.factor(SAI_cat), data = subset(flycomp, !(flycomp$SAI_cat == 'medium')))
lvh <- median_test(Recombination~as.factor(SAI_cat), data = subset(flycomp, !(flycomp$SAI_cat == 'medium')))
#medium v high
median_test(Recombination~as.factor(SAI_cat), data = subset(flycomp, !(flycomp$SAI_cat == 'low')))
mvh <- median_test(Recombination~as.factor(SAI_cat), data = subset(flycomp, !(flycomp$SAI_cat == 'low')))

#prepare median test for plotting
#p value must be entered manually
fmedtesting<- data.frame(group1 = c('high', 'high', 'low'), group2 = c('low', 'medium', 'medium'), p = c(0.0445, 0.597, 3.372e-06),
                        p.code = c('*', '.', '***'), y.position = c(7, 8, 6), Species = c('Drosophila', 'Drosophila', 'Drosophila'))
fmedtesting$p.plot<- signif(fmedtesting$p, 3)

#prepare human testing for plotting
hrectest<- data.frame(group1 = c('high', 'high', 'low'), group2 = c('low', 'medium', 'medium'), p = c(0.2531, 0.2531, 0.2531),
                      p.code = c('ns', 'ns', 'ns'), y.position = c(0.25, 0.4, 0.1), Species = c('Human', 'Human', 'Human'))
hvartest<- data.frame(group1 = c('high', 'high', 'low'), group2 = c('low', 'medium', 'medium'), p = c(0.9067, 0.9067, 0.9067),
                      p.code = c('ns', 'ns', 'ns'), y.position = c(-1.15, -1, -1.3), Species = c('Human', 'Human', 'Human'))

##compare recombination and SAI between humans and flies
median_test(Recombination~as.factor(Species), data = humanfly)
leveneTest(Recombination~as.factor(Species), data = humanfly)
median_test(SAI~as.factor(Species), data = humanfly)
leveneTest(SAI~as.factor(Species), data = humanfly)

###loop for all consequence types
#empty dataframe
cat_genes<- data.frame(Gene = character(), SAI = numeric(), Recombination = numeric(), 
                       SAI_cat = character(), t_consequence = character())
#dataframe will contain every SA gene, grouped by consequence, SAI averaged for that variant type only
#calculate SAI
for(i in 1:nlevels(as.factor(mergesnps$con_class))){
  t_con_class<- levels(as.factor(mergesnps$con_class))[i]
  t_snps<- subset(mergesnps, mergesnps$con_class == t_con_class)
  #mean SAI by gene
  t_genes<- setNames(aggregate(t_snps$SAI~t_snps$Gene, FUN = mean), c("Gene", "SAI"))
  #get recombination rates
  t_genes<- merge(t_genes, sa_genes[, c('Gene', 'Mean.Recombination.Comeron')], by = 'Gene', all.x = T)
  #categorise SAI
  colnames(t_genes)[which(names(t_genes) == "Mean.Recombination.Comeron")] <- "Recombination"
  for(i in 1:nrow(t_genes)){
    if(t_genes$SAI[i]<= quantile(t_genes$SAI, .33)){
      t_genes$SAI_cat[i]<- "low" 
    }else if(t_genes$SAI[i]>=quantile(t_genes$SAI, .67)){
      t_genes$SAI_cat[i]<- "high"
    }else{
      t_genes$SAI_cat[i]<- "medium"
    }
  }
  #add into dataframe
  t_genes$con_class <- t_con_class
  cat_genes<- rbind(cat_genes, t_genes)
}  
#test SAI
#empty dataframes
mediantest1<- data.frame(Test = character(), con_class = character(),
                         p.value = numeric(), chisquared = numeric())
mediantests<- data.frame(Test = character(), con_class = character(),
                         Comparison = character(), p.value = numeric(),
                         p.adjust = numeric())
levenetests<- data.frame(Test = character(), con_class = character(),
                         Comparison = character(), DF = numeric(),
                         f = numeric(), P = numeric())
#loop for each variant type
for(i in 1:nlevels(as.factor(cat_genes$con_class))){
  try({
    t_con_class<- levels(as.factor(cat_genes$con_class))[i]
    t_cat<- subset(cat_genes, cat_genes$con_class == t_con_class)
    temp_test1<- median_test(Recombination~as.factor(SAI_cat), data = t_cat)
    print(t_con_class)
    conprint<- print(temp_test1)
    #  temp_med1<- data.frame(Test = "Median", con_class = t_con_class, 
    #                         p.value = temp_test1$p-value, chisquared = temp_test1@statistic@teststatistic)
    temp_test<- pairwiseMedianTest(Recombination~as.factor(SAI_cat), data = t_cat)
    temp_med<- data.frame(Test = "Pairwise median", con_class = t_con_class, 
                          Comparison = temp_test$Comparison,
                          p.value = temp_test$p.value, p.adjust = temp_test$p.adjust)
    #  mediantest1<- rbind(mediantest1, temp_med1)
    mediantests<- rbind(mediantests, temp_med)
    #levene's test
    #low v medium
    temp_test <- leveneTest(Recombination~as.factor(SAI_cat), data = subset(t_cat, !(t_cat$SAI_cat == 'high')))
    temp_lev<- data.frame(Test = "Levene test", con_class = t_con_class, 
                          Comparison = 'low v medium', DF = temp_test$Df[2],
                          f = temp_test$`F value`[1], P = temp_test$`Pr(>F)`[1])
    levenetests<- rbind(levenetests, temp_lev)
    #low v high
    temp_test<- leveneTest(Recombination~as.factor(SAI_cat), data = subset(t_cat, !(t_cat$SAI_cat == 'medium')))
    temp_lev<- data.frame(Test = "Levene test", con_class = t_con_class, 
                          Comparison = 'low v high', DF = temp_test$Df[2],
                          f = temp_test$`F value`[1], P = temp_test$`Pr(>F)`[1])
    levenetests<- rbind(levenetests, temp_lev)
    #medium v high
    temp_test<- leveneTest(Recombination~as.factor(SAI_cat), data = subset(t_cat, !(t_cat$SAI_cat == 'low')))
    temp_lev<- data.frame(Test = "Levene test", con_class = t_con_class, 
                          Comparison = 'medium v high', DF = temp_test$Df[2],
                          f = temp_test$`F value`[1], P = temp_test$`Pr(>F)`[1])
    levenetests<- rbind(levenetests, temp_lev)
  })
}
#make dataframe of Brown-Mood p values- have to be entered manually
medianbmood<- data.frame(con_class = levels(as.factor(mediantests$con_class)),
                         P = c(0.544, #missense variant
                               0.0000291, #regulatory variant
                               0.08444 #synonymous variant
                              ))

#data for plotting
mediantests$group1<- vapply(strsplit(mediantests$Comparison," "), `[`, 1, FUN.VALUE=character(1))
mediantests$group2<- vapply(strsplit(mediantests$Comparison," "), `[`, 3, FUN.VALUE=character(1))
levenetests$group1<- vapply(strsplit(levenetests$Comparison," "), `[`, 1, FUN.VALUE=character(1))
levenetests$group2<- vapply(strsplit(levenetests$Comparison," "), `[`, 3, FUN.VALUE=character(1))
#add y values for p value labels
#mediantests$Recombination <- 6
for(i in 1:nrow(mediantests)){
  if(mediantests$group1[i] == 'high'){
    if(mediantests$group2[i] == 'medium'){
      mediantests$Recombination[i] <- 8
    }else{
      mediantests$Recombination[i] <- 7
    }
  }else{
    mediantests$Recombination[i] <- 6
  }
}
#levenetests$Recombination <- -1
for(i in 1:nrow(levenetests)){
  if(levenetests$group1[i] == 'low'){
    if(levenetests$group2[i] == 'medium'){
      levenetests$Recombination[i] <- -1
    }else{
      levenetests$Recombination[i] <- -2
    }
  }else{
    levenetests$Recombination[i] <- -3
  }
}
#significant figures
mediantests$p.adjust.sig<- signif(mediantests$p.adjust, 3)
levenetests$p.adjust.sig<- signif(levenetests$P, 3)
#plot only significant median test values
pmedianbmood<- subset(medianbmood, medianbmood$P <= 0.05)
pmediantests<- subset(mediantests, mediantests$con_class %in% pmedianbmood$con_class)
medianbmood$Recombination <- 8.5 #y value axis for plotting
medianbmood$Pplot <- paste('P =', signif(medianbmood$P, 3), sep = ' ')
#change SAI order for plotting, to reflect conflict intensity
sorder<- c('high', 'medium', 'low')

###Paralogs###

##load dataset
paralogs<- fread("./Data/paralog_recombination.txt", sep = '\t', header = T)

##All melanogaster genes with paralogs
#load data for every dmel paralog
all_paralogs<- fread("./Data/flybase/dmel_paralogs_fb_2021_02.csv", header = T, sep = ',')
#convert locus for flybase back-conversion tool
all_paralogs$Gene_location_full<- locationwriter(all_paralogs$Arm.Scaffold, all_paralogs$Location)
all_paralogs$Paralog_location_full<- locationwriter(all_paralogs$Paralog_Arm.Scaffold, all_paralogs$Paralog_Location)

parent_convert<- fread("./Data/flybase/Flybase_Converted_Coordinates_parents.tsv", sep = "\t", header = F)
paralog_convert<- fread("./Data/flybase/Flybase_Converted_Coordinates_all_paralogs.tsv", sep = "\t", header = F)
names(parent_convert)<- c("R6", "R5", "notes")
names(paralog_convert)<- c("R6", "R5", "notes")
#put back into paralog table
all_paralogs$parent_r5<- parent_convert$R5
all_paralogs$paralog_r5<- paralog_convert$R5
#get r5 coordinates
all_paralogs$parent_chromosome_r5<- locationchromosomeextract(all_paralogs$parent_r5)
all_paralogs$parent_location_start_r5<- locationsplitstart(all_paralogs$parent_r5)
all_paralogs$parent_location_end_r5<- locationsplitend(all_paralogs$parent_r5)
all_paralogs$paralog_chromosome_r5<- locationchromosomeextract(all_paralogs$paralog_r5)
all_paralogs$paralog_location_start_r5<- locationsplitstart(all_paralogs$paralog_r5)
all_paralogs$paralog_location_end_r5<- locationsplitend(all_paralogs$paralog_r5)

#filter with DIOPT score
all_paralogs<-subset(all_paralogs, all_paralogs$DIOPT_score>=2)

#remove error values
all_paralogs_correct<- subset(all_paralogs, all_paralogs$parent_chromosome_r5 %in% comeron$`Chromosome Arm`)
all_paralogs_correct<- subset(all_paralogs_correct, all_paralogs_correct$paralog_chromosome_r5 %in% comeron$`Chromosome Arm`)
all_paralogs_correct<- subset(all_paralogs_correct, !is.na(all_paralogs_correct$parent_location_start_r5))
all_paralogs_correct<- subset(all_paralogs_correct, !is.na(all_paralogs_correct$paralog_location_start_r5))

#slow steps- skip to read line to use premade file. Code used to make file below.
 # for(i in 1:nrow(all_paralogs_correct)){ #3 min
 #   all_paralogs_correct$parent_recombination[i]<- recomlookup(all_paralogs_correct$parent_chromosome_r5[i], all_paralogs_correct$parent_location_start_r5[i], all_paralogs_correct$parent_location_end_r5[i], comeron)
 #   print(i)
 # }
 # for(i in 1:nrow(all_paralogs_correct)){ #3 min
 #   all_paralogs_correct$paralog_recombination[i]<- recomlookup(all_paralogs_correct$paralog_chromosome_r5[i], all_paralogs_correct$paralog_location_start_r5[i], all_paralogs_correct$paralog_location_end_r5[i], comeron)
 #   print(i)
 # }
 # fwrite(all_paralogs_correct, './Output/all_paralogs_correct.csv')
all_paralogs_correct<- read.csv('./Output/all_paralogs_correct.csv')

#Parent versus paralog recombination model
paralog_recombination_model<- lm(data = all_paralogs_correct, paralog_recombination~parent_recombination)
summary(paralog_recombination_model)

#get counts
paralog_count<- as.data.frame(table(all_paralogs$FBgn_ID))
names(paralog_count) <- c("FBid", "paralog_count")
##melanogaster genes with and without paralogs
w_paralogs<- subset(dmel_genes_correct, (dmel_genes_correct$primary_FBid %in% all_paralogs$FBgn_ID))
w_paralogs$Paralogs <- "Y"
w_paralogs<- merge(w_paralogs, paralog_count, by.x = "primary_FBid", by.y = "FBid")
no_paralogs<- subset(dmel_genes_correct, !(dmel_genes_correct$primary_FBid %in% w_paralogs$primary_FBid))
no_paralogs$Paralogs <- "N"
no_paralogs$paralog_count<- 0
allplusnone<- rbind(no_paralogs, w_paralogs)

##mean by gene families
all_gene_families<- aggregate(all_paralogs_correct, by = list(all_paralogs_correct$FBgn_ID), FUN = mean)

#compare recombination rates of drosophila genes with and w/o paralogs
all_paralogs_correct<- subset(all_paralogs_correct, all_paralogs_correct$FBgn_ID %in% dmel_genes_correct$primary_FBid)
all_gene_families2<- subset(all_gene_families, all_gene_families$Group.1 %in% dmel_genes_correct$primary_FBid)
plot(density(all_gene_families$paralog_recombination))
lines((density(subset(no_paralogs$recombination_rate, !is.na(no_paralogs$recombination_rate)))))
all_gene_families$P<- "Y"
no_paralogs$P<- "N"
no_paralogs<-na.omit(no_paralogs)
all_none<- data.frame(all_gene_families$Group.1, all_gene_families$paralog_recombination, all_gene_families$P)
names(all_none)<- c("Gene", "Recombination", "Paralogs")
noparas<- data.frame(no_paralogs$primary_FBid, no_paralogs$recombination_rate, no_paralogs$P)
names(noparas)<- c("Gene", "Recombination", "Paralogs")
all_none2<- rbind(all_none, noparas)
#tests for paralog/no paralog comparison
wilcox.test(all_none2$Recombination~all_none2$Paralogs)
t.test(Recombination ~ Paralogs, all_none2)

#subset by consequence
missense_genes<- subset(sa_genes, sa_genes$Gene %in% subset(mergesnps, mergesnps$con_class == "missense_variant")$Gene)

#Compare SA genes with and w/o paralogs
sa_w_paralogs<- subset(w_paralogs, w_paralogs$primary_FBid %in% missense_genes$Gene)
sa_no_paralogs<- subset(no_paralogs, no_paralogs$primary_FBid %in% missense_genes$Gene)
sa_w_paralogs$Paralogs<- "Y"
sa_no_paralogs$Paralogs<- "N"
sa_no_paralogs<- subset(sa_no_paralogs, select = -c(P))
sa_all_none<-rbind(sa_w_paralogs, sa_no_paralogs)

#count number of each consequence type with and without paralogs
consequencecount<- NULL
for (i in 1:length(levels(as.factor(mergesnps$con_class)))){
  tsubset<- subset(sa_genes, sa_genes$Gene %in% subset(mergesnps, mergesnps$con_class == levels(as.factor(mergesnps$con_class))[i])$Gene)
  consequencecount$con_class[i] <- levels(as.factor(mergesnps$con_class))[i]
  consequencecount$Total[i]<- count(tsubset$Gene %in% dmel_genes$primary_FBid)
  consequencecount$W_paralogs[i]<- count(tsubset$Gene %in% w_paralogs$primary_FBid)
  consequencecount$No_paralogs[i]<- count(tsubset$Gene %in% no_paralogs$primary_FBid)
}
consequencecount<- data.frame(consequencecount)
consequencecount$proportion_paralogs<- consequencecount$W_paralogs/consequencecount$Total
#distinguish significant models
for (i in 1:nrow(consequencecount)){
  if(consequencecount$con_class[i] %in% significant_models$con_class){
    consequencecount$model_significance[i]<- "Y"
  }else{
    consequencecount$model_significance[i]<- "N"
  }
}
#chi square test
testparalogs<- data.frame(nrow(w_paralogs), nrow(no_paralogs))
names(testparalogs)<- c("W_paralogs", "No_paralogs")
for(i in 1:nrow(consequencecount)){
  testparalogs[2,]<- c(consequencecount$W_paralogs[i], consequencecount$No_paralogs[i])
  conchi<- chisq.test(testparalogs)
  consequencecount$chi[i]<- conchi$statistic
  consequencecount$chi_p[i]<- conchi$p.value
  consequencecount$chi_df[i]<- conchi$parameter
}

fwrite(consequence_models, "./Output/consequencemodels.csv", sep = ',', quote = F)
fwrite(consequencecount, "./Output/consequencecount.csv", sep = ',', quote = F)

####################################################
#                    4. Plotting                   #
####################################################


#combined minor model
p1<- ggplot(data = mergesnps, mapping = aes(x = effratio, y = EAF)) +
  #  geom_point(size = 1, alpha = 0.1)+
  ylim(0,0.55)+
  xlim(-7.5,0)+
  scale_y_continuous(breaks = c(0.1, 0.2, 0.3, 0.4, 0.5))+
  labs(x = "Effect size ratio (positive effect over negative)", y = "Effect allele frequency")+
  annotate(geom = "text", x = -1.7, y = 0.55, size = 2.9, label = "Net\n Positive")+
  annotate(geom = "text", x = -0.3, y = 0.55, size = 2.9, label = "Net\n Negative")+
  #  stat_smooth()
  stat_binhex() +
  scale_fill_gradient(low = "lightblue", high = "red", limits = c(0, 150))+
  geom_vline(xintercept = -1, linetype = "dashed", colour = "black")+ #modified to log10(1)
  geom_line(color = 'black', data = main_predict, aes(x = effratio, y = EAF))+
  theme(text = element_text(size = 12),
        axis.line = element_line(),
        panel.background = element_rect(fill = "white"),
        panel.border = element_rect(colour = 'black', fill = NA, linewidth = 1))
ggsave('Figure 1.png', p1, width = 4.7, height = 4, dpi = 300)

#Main plot- by consequence class
p2<- ggplot(data = mergesnps, mapping = aes(x = effratio, y = EAF)) +
  #  geom_point(size = 4, alpha = 0.1)+
  facet_wrap(~con_class, ncol = 3, labeller = labeller(con_class = consequence_labels))+
  stat_binhex() +
  scale_fill_gradient(low = "lightblue", high = "red", limits = c(0, 60))+
  geom_vline(xintercept = -1, linetype = "dashed", colour = "black")+ #modified to log10(1)
  labs(x = "Effect size ratio (positive effect over negative)", y = "Effect allele frequency")+
  geom_line(color = 'black', data = count_all_predicted, aes(x = effratio, y = count_over_sample_predicted))+
  theme(text = element_text(size = 24),
        axis.line = element_line(),
        panel.background = element_rect(fill = "white"),
        panel.border = element_rect(colour = 'black', fill = NA, size = 1),
        legend.position = "none",
        strip.text.x = element_text(size = 15)) #facet label text
ggsave('Figure 2.png', p2, width = 10, height = 4, dpi = 300)

#combined recombination boxplot
plot1<- ggplot(data = flycomp, aes(x = factor(SAI_cat, level = sorder), y = Recombination))+
  geom_boxplot(fill = '#F8766D', outlier.shape = NA)+
  geom_jitter(shape=16, alpha = 0.5)+
  labs(x = "Conflict intensity (SAI)", y = "Recombination rate (cM/Mb)")+
  stat_pvalue_manual(fmedtesting, size = 2.5, label = 'p.plot', y.position = 'y.position')+
  scale_x_discrete(labels = c('low', 'medium', 'high'))+
  ylim(-0.1,8.5)+
  theme(legend.position = 'none')+
  theme_classic()
plot2<- ggplot(data = humancomp_unique, aes(x = factor(SAI_cat, level = sorder), y = Recombination))+
  geom_boxplot(fill = '#00BFC4', outlier.shape = NA)+
  geom_jitter(shape=16, alpha = 0.5)+
  labs(x = "Conflict intensity (SAI)", y = "Recombination rate (cM/Mb)")+
  scale_x_discrete(labels = c('low', 'medium', 'high'))+
  ylim(-0.1,8.5)+
  theme(legend.position = 'none')+
  theme_classic()
p3<- plot_grid(plot1, plot2, labels = "AUTO")
ggsave('Figure 3.png', p3, width = 8, height = 4.5, dpi = 300)

#fly and human SAI comparison
p4<- ggplot(data = humanfly, aes(x = Species, y = SAI, fill = Species))+
  geom_boxplot(outlier.shape = NA)+
  scale_fill_discrete(c('#F8766D', '#00BFC4'), guide = 'none')+
  geom_jitter(shape=16, alpha = 0.5)+
  labs(x = "Species", y = "SAI")+
  theme(legend.position = 'none')+
  coord_cartesian(ylim = c(-2, 0), xlim = c(1, 2), clip="off")+ 
  annotate('segment', x = 0.01, xend = 0.01, y = -2, yend = 0,
           arrow=arrow(length=unit(0.3, "cm"), type = 'closed', ends = 'both'))+
  annotate('text', x = -0.08, y = -1.8, label = 'High', angle = 90)+
  annotate('text', x = -0.08, y = -0.17, label = 'Low', angle = 90)+
  annotate('text', x = -0.1, y = -1, label = 'Conflict intensity', angle = 90)+
  theme_classic()+
  theme(plot.margin = margin(l = 25, unit = "pt"))
ggsave('Figure S1.png', p4, width = 4, height = 4.5, dpi = 300)

#Recombination boxplot- fly faceted by variant type
p5<- ggplot(data = cat_genes, aes(x = factor(SAI_cat, level = sorder), y = Recombination))+
  facet_wrap(~con_class, ncol = 3, labeller = labeller(con_class = consequence_labels))+
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(shape=16, alpha = 0.5)+
  labs(x = "Conflict intensity (SAI)", y = "Recombination rate (cM/Mb)")+
  scale_x_discrete(labels = c('low', 'medium', 'high'))+
  #median comparison
  stat_pvalue_manual(pmediantests, size = 2.5, label = 'p.adjust.sig', y.position = 'Recombination')+
  geom_text(data = medianbmood, mapping = aes(x = 'medium', y = Recombination, label = Pplot, size = 20))+
  theme_classic()+
  theme(text = element_text(size = 20),
        legend.position = 'none',
        strip.text.x = element_text(size = 15))
ggsave('Figure 4.png', p5, width = 8, height = 6, dpi = 300)

p6<- ggplot(data = all_none2, mapping = aes(x = Paralogs, y = Recombination))+
  geom_violin()+
  labs(y="Recombination rate (cM/Mb)")+
  geom_jitter(shape=16, position=position_jitter(0.4), alpha = 0.01)+
  theme_bw()
ggsave('Figure 5.png', p6, width = 4, height = 4, dpi = 300)

#########################Code for revisions#####################################
#count number of paralogs
tcount<- NULL
for(i in levels(as.factor(cat_genes$Gene))){
  tsub<- subset(all_paralogs_correct, all_paralogs_correct$FBgn_ID == i)
  tcount<- rbind(tcount, data.frame(Gene = i, N_paralogs = nrow(tsub)))
}
#merge column
cat_genes2<- merge(cat_genes, tcount, by = 'Gene')
ggplot(cat_genes2, aes(SAI_cat, N_paralogs))+
  geom_violin()+
  geom_jitter(shape=16, alpha = 0.5)
ks.test(subset(cat_genes2, cat_genes2$SAI_cat == 'low')$N_paralogs, subset(cat_genes2, cat_genes2$SAI_cat == 'medium')$N_paralogs)
#reshape for plotting
cat_percent<- NULL
for(i in 1:nrow(con_count2)){
  cat_percent<- rbind(cat_percent, data.frame(SAI_cat = c(con_count2$SAI_cat[i], con_count2$SAI_cat[i]), Paralogs = c('paralogs', 'no paralogs'), 
                                              Percentage = c(con_count2$par[i]/(con_count2$par[i] + con_count2$npar[i]), con_count2$npar[i]/(con_count2$par[i] + con_count2$npar[i]))))
}
#separate by SAI
con_count2<- NULL
a <- 7085
b <- 9042
testparalogs<- data.frame(a, b)
for(j in levels(as.factor(cat_genes$SAI_cat))){
  tsub2<- subset(cat_genes, cat_genes$SAI_cat == j)
  par<- nrow(subset(tsub2, tsub2$Gene %in% all_paralogs_correct$FBgn_ID))
  npar<- nrow(subset(tsub2, !(tsub2$Gene %in% all_paralogs_correct$FBgn_ID)))
  testparalogs[2,]<- c(par, npar)
  conchi<- chisq.test(testparalogs)
  con_count2<- rbind(con_count2, data.frame(SAI_cat = j, proportion_paralogs = par/(par+npar), par = par, npar = npar, chi = conchi$statistic, chi_p = conchi$p.value, chi_df = conchi$parameter))
}

p7<- ggplot(cat_percent, aes(x = factor(SAI_cat, level = sorder), Percentage, fill = Paralogs))+
  geom_bar(stat = 'identity', position = 'stack', width = 1, colour = 'black')+
  geom_hline(yintercept = (7085/(7085+9042)), linetype = 'dashed')+
  labs(y = 'Proportion', x = 'Conflict intensity')+
  scale_x_discrete(expand = c(0,0), labels = c('weak', 'medium', 'strong'))+
  scale_y_continuous(expand = c(0,0))+
  theme_classic()
ggsave('Figure 7.png', p7, width = 4, height = 4, dpi = 300)
fwrite(con_count2, 'SA_paralog_table.csv')

sa_pars<- subset(cat_genes2, cat_genes2$N_paralogs>0)
ks.test(subset(sa_pars, sa_pars$SAI_cat == 'medium')$N_paralogs, subset(sa_pars, sa_pars$SAI_cat == 'high')$N_paralogs)
t1<- wilcox.test(subset(sa_pars, sa_pars$SAI_cat == 'medium')$N_paralogs, subset(sa_pars, sa_pars$SAI_cat == 'high')$N_paralogs)
t2<- wilcox.test(subset(sa_pars, sa_pars$SAI_cat == 'medium')$N_paralogs, subset(sa_pars, sa_pars$SAI_cat == 'low')$N_paralogs)
t3<- wilcox.test(subset(sa_pars, sa_pars$SAI_cat == 'low')$N_paralogs, subset(sa_pars, sa_pars$SAI_cat == 'high')$N_paralogs)
#compile tests
para_test<- NULL
para_test<- data.frame(comp1 = 'Medium', comp2 = 'Weak', P = t1$p.value, W = t1$statistic)
para_test<- rbind(para_test, data.frame(comp1 = 'Medium', comp2 = 'Strong', P = t2$p.value, W = t2$statistic))
para_test<- rbind(para_test, data.frame(comp1 = 'Strong', comp2 = 'Weak', P = t3$p.value, W = t3$statistic))
fwrite(para_test, 'SA_paralog_number_tests.csv')
wilcox.test(subset(sa_pars, sa_pars$SAI_cat == 'medium')$N_paralogs, subset(sa_pars, sa_pars$SAI_cat == 'low')$N_paralogs)
#calculate number of paralogs for all drosophila genes
genome_pars<- NULL
for(i in levels(as.factor(all_paralogs_correct$FBgn_ID))){
  tsub<- subset(all_paralogs_correct, all_paralogs_correct$FBgn_ID == i)
  genome_pars<- rbind(genome_pars, data.frame(Gene = i, N_paralogs = nrow(tsub), SAI_cat = 'all genes with paralogs'))
}
m1<- lm(N_paralogs~SAI_cat, sa_pars)
summary(m1)


p8<- ggplot(sa_pars, aes(x = factor(SAI_cat, level = sorder), N_paralogs))+
  labs(y = 'Number of paralogs', x = 'Conflict intensity', colour = element_blank())+
  scale_x_discrete(expand = c(0,0), labels = c('low', 'medium', 'high', 'all genes'))+
  scale_y_continuous(expand = c(0,0))+
  scale_colour_viridis(option = 'turbo')+
  geom_violin()+
  geom_jitter(shape=16, alpha = 0.5, aes(colour = N_paralogs))+
  geom_violin(data = genome_pars)+
  geom_jitter(data = genome_pars, shape=16, alpha = 0.05, aes(colour = N_paralogs))+
  theme_classic()
#ggsave('Figure 6.png', p8, width = 4, height = 4, dpi = 300)

##########Revisions#########

tmerge<- merge(flycomp, dmel_genes_correct, by.x = 'Gene', by.y = 'primary_FBid')
summary(as.factor(subset(tmerge, tmerge$chrom.r5 == 'X')$SAI_cat))
summary(as.factor(subset(tmerge, tmerge$chrom.r5 != 'X')$SAI_cat))
summary(subset(tmerge, tmerge$chrom.r5 == 'X')$recombination_rate)
summary(subset(tmerge, tmerge$chrom.r5 != 'X')$recombination_rate)
ks.test(subset(tmerge, tmerge$chrom.r5 == 'X')$recombination_rate, subset(tmerge, tmerge$chrom.r5 != 'X')$recombination_rate)
ggplot(subset(tmerge, tmerge$chrom.r5 == 'X'), aes(1, recombination_rate))+
  geom_boxplot()+
  geom_boxplot(subset(tmerge, tmerge$chrom.r5 != 'X'), mapping = aes(2, recombination_rate))


#test paralog counts of genes on x chromosome
testparalogs<- data.frame(nrow(w_paralogs), nrow(no_paralogs))
names(testparalogs)<- c("W_paralogs", "No_paralogs")
x_genes<- subset(dmel_genes_correct, dmel_genes_correct$chrom.r5 == 'X')
x_paras<- subset(x_genes, x_genes$primary_FBid %in% w_paralogs$primary_FBid)
x_noparas<- subset(x_genes, x_genes$primary_FBid %in% no_paralogs$primary_FBid)
testparalogs[2,] <- c(nrow(x_paras), nrow(x_noparas))
chisq.test(testparalogs)
chisq.test(x = c(nrow(x_paras), nrow(x_noparas)), p = c(a,b))
testparalogs[1,1]/(testparalogs[1,1]+testparalogs[1,2])
testparalogs[2,1]/(testparalogs[2,1]+testparalogs[2,2])

nrow(subset(sa_genes, (sa_genes$Gene %in% w_paralogs$primary_FBid)))

#separate by SAI
con_count2<- NULL
a <- 7085
b <- 9042
testparalogs<- data.frame(a, b)
for(j in levels(as.factor(cat_genes$SAI_cat))){
  tsub2<- subset(cat_genes, cat_genes$SAI_cat == j)
  par<- nrow(subset(tsub2, tsub2$Gene %in% all_paralogs_correct$FBgn_ID))
  npar<- nrow(subset(tsub2, !(tsub2$Gene %in% all_paralogs_correct$FBgn_ID)))
  testparalogs[2,]<- c(par, npar)
  conchi<- chisq.test(testparalogs)
  con_count2<- rbind(con_count2, data.frame(SAI_cat = j, proportion_paralogs = par/(par+npar), par = par, npar = npar, chi = conchi$statistic, chi_p = conchi$p.value, chi_df = conchi$parameter))
}

#count number of paralogs
tcount<- NULL
for(i in levels(as.factor(cat_genes$Gene))){
  tsub<- subset(all_paralogs_correct, all_paralogs_correct$FBgn_ID == i)
  tcount<- rbind(tcount, data.frame(Gene = i, N_paralogs = nrow(tsub)))
}
#merge column
cat_genes2<- merge(cat_genes, tcount, by = 'Gene')
ggplot(cat_genes2, aes(SAI_cat, N_paralogs))+
  geom_violin()+
  geom_jitter(shape=16, alpha = 0.5)
ks.test(subset(cat_genes2, cat_genes2$SAI_cat == 'low')$N_paralogs, subset(cat_genes2, cat_genes2$SAI_cat == 'medium')$N_paralogs)
#reshape for plotting
cat_percent<- NULL
for(i in 1:nrow(con_count2)){
  cat_percent<- rbind(cat_percent, data.frame(SAI_cat = c(con_count2$SAI_cat[i], con_count2$SAI_cat[i]), Paralogs = c('paralogs', 'no paralogs'), 
                                              Percentage = c(con_count2$par[i]/(con_count2$par[i] + con_count2$npar[i]), con_count2$npar[i]/(con_count2$par[i] + con_count2$npar[i]))))
}

p7<- ggplot(cat_percent, aes(x = factor(SAI_cat, level = sorder), Percentage, fill = Paralogs))+
  geom_bar(stat = 'identity', position = 'stack', width = 1, colour = 'black')+
  geom_hline(yintercept = (7085/(7085+9042)), linetype = 'dashed')+
  labs(y = 'Proportion', x = 'Conflict intensity')+
  scale_x_discrete(expand = c(0,0), labels = c('low', 'medium', 'high'))+
  scale_y_continuous(expand = c(0,0))+
  theme_classic()
#ggsave('Figure 7.png', p7, width = 4, height = 4, dpi = 300)
p9<- plot_grid(p7, p8, labels = "AUTO")
ggsave('Figure 6.png', p9, width = 8, height = 4, dpi = 300)
fwrite(con_count2, 'SA_paralog_table.csv')