library(tidyverse)
library(admixtools)
library(igraph)
library(stringr)
snps<-read.table("~/Downloads/old_zln.snp")  #read snps data, make chr real value based on marker id
colnames(snps)=c("id","chr","cm","bp","allele1","allele2")
snps<-mutate(snps,chr=as.numeric(gsub("_[0-9]+","",snps$id,perl=TRUE)))

ogut<-read.table("~/ogut_fifthcM_map_agpv4_INCLUDE.txt") #read in genetic map
colnames(ogut)=c("id","marker","cm","chr","bp")

#fit spline and get cm values for each chr
map<-data.frame()
for(i in 1:10){
  mapi<-filter(ogut,chr==i) %>% arrange(bp)
  snpsi<-filter(snps,chr==i)
  mapi$cm=mapi$cm+3*abs(min(mapi$cm)) # ensure no negative cM (maybe make this bigger)
  splinei<-smooth.spline(x=mapi$bp,y=mapi$cm)
  map=rbind(map,data.frame(snpsi$id,snpsi$chr,predict(splinei,x=snpsi$bp)$y,snpsi$bp,snpsi$allele1,snpsi$allele2))
}
colnames(map)=colnames(snps)
write.table(map,file="~/Downloads/zln.snp",quote=F,row.names=F,col.names=F)

get pops, calculate f2

mypops<-c("Andean","Mexlow","mex","dip","balsas")
path_to_geno<-"~/Downloads/zln"
my_f2_dir = "~/Desktop/my_f2"
extract_f2(path_to_geno, my_f2_dir,overwrite=TRUE,pops=mypops,verbose=FALSE)
f2_blocks = f2_from_precomp(my_f2_dir,pops=mypops,verbose=FALSE)

two teos tree. fit to f2 data

jri_tree<-read_table("~/Desktop/jri_tree.txt")
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   from = col_character(),
##   to = col_character()
## )
jri<-qpgraph(f2_blocks,jri_tree,return_fstats=TRUE)
plot_graph(jri$edges,textsize=5,title="two teosintes tree",fix=TRUE)

filter(jri$f4,pop3=="Andean",pop4=="Mexlow") %>% select(1:4,9,10)

alternative – maize domesticated from hybrids. fit to data

mbh_tree<-read_table("~/Desktop/simple_mbh.tree")
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   from = col_character(),
##   to = col_character()
## )
mbh<-qpgraph(f2_blocks,mbh_tree,return_fstats=TRUE)
plot_graph(mbh$edges,textsize=5,title="domestication from parv-mex hybrid")

filter(mbh$f4,pop3=="Andean",pop4=="Mexlow") %>% select(1:4,9,10)

find best graph with hybrid starting point and 2 admixtures. result shows no admixture with Andes!

winner_mbh = find_graphs(f2_blocks,initgraph=graph_from_data_frame(mbh_tree, directed = TRUE),stop_gen = 100,verbose=FALSE) %>%
  slice_min(score, with_ties = FALSE)
plot_graph(winner_mbh$edges[[1]],textsize=5,title="find graphs from mbh")

mbh_found=qpgraph(f2_blocks,graph_from_data_frame(winner_mbh$edges, directed = TRUE),return_fstats=TRUE)
filter(mbh_found$f4,pop3=="Andean",pop4=="Mexlow") %>% select(1:4,9,10)

simple tree, fit to data

simple_tree<-read_table("~/Desktop/simple.tree")
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   from = col_character(),
##   to = col_character()
## )
simple<-qpgraph(f2_blocks,simple_tree,return_fstats=TRUE)
plot_graph(simple$edges,textsize=5,title="simple graph with no admixture")

filter(simple$f4,pop3=="Andean",pop4=="Mexlow") %>% select(1:4,9,10)

best graph of simple tree is same as above

winner_simple<-find_graphs(f2_blocks,initgraph=graph_from_data_frame(simple_tree, directed = TRUE),stop_gen = 100,verbose=FALSE) %>% 
  slice_min(score, with_ties = FALSE)
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
## x No new models!
plot_graph(winner_simple$edges[[1]],textsize=5)

calculate scores. in theory corrects for differences in model complexity

fits<-qpgraph_resample_multi(f2_blocks, list(mbh_tree,simple_tree,winner_mbh$edges[[1]],jri_tree), verbose=TRUE, nboot = 200)

all models fit better than the simple

print(c(compare_fits(fits[[4]]$score_test, fits[[2]]$score_test)$p_emp, 
compare_fits(fits[[3]]$score_test, fits[[2]]$score_test)$p_emp,  
compare_fits(fits[[1]]$score_test, fits[[2]]$score_test)$p_emp))  
## [1] 0.005 0.005 0.005

support for two teos fits better than either alternative model

print(c(compare_fits(fits[[4]]$score_test, fits[[3]]$score_test)$p_emp,
compare_fits(fits[[4]]$score_test, fits[[3]]$score_test)$p_emp))
## [1] 0.005 0.005

hybrid origin vs no admixture in andes finds no difference

compare_fits(fits[[1]]$score_test, fits[[3]]$score_test)$p_emp
## [1] 0.38