R Notebook

library("LEA")
library("cluster")
library("cowplot")
library("tidyverse")
source("http://membres-timc.imag.fr/Olivier.Francois/Conversion.R")
source("http://membres-timc.imag.fr/Olivier.Francois/POPSutilities.R")

## [1] "Loading fields"
## [1] "Loading RColorBrewer"

source("http://peterhaschke.com/Code/multiplot.R")

gendata<-function(migration_command,n){
  hapdat<-matrix(unlist(lapply(strsplit(system(paste("~/src/msdir/ms ", n*10, " 1 -t 10 -r 10 1000 -I ", n, " ", paste(rep(10,3),collapse=" ") , " ", migration_command, " | tail -", n*10,sep=""),intern=T),''),as.numeric)),nrow=n*10,byrow=T)
  genodat<-hapdat[seq(1,n*10,2),]+hapdat[seq(2,n*10,2),] #make diploid
  write.table(t(genodat),file="~/crap.geno",quote=F,row.names=F,col.names=F,sep="")
  return(genodat)
}

npops=3
for(i in 1:100){
  migration_command=paste(round(runif(1),3),"-ej 0.001 2 1 -ej 0.03 3 1")
  mydat<-gendata(migration_command,npops)
  
  #cut based on silhouette scores
  keep<-sort(which(silhouette(kmeans(dist(mydat),npops)[[1]],dist(mydat))[,npops]>0))
  left_out<-which(!(1:15 %in% keep))
  if(length(left_out)<1){ next }
  
  #make new data from cut
  cutdat<-mydat[keep,]
  write.table(t(cutdat),file="~/crapcut.geno",quote=F,row.names=F,col.names=F,sep="")
  
  #structure on original
  obj.snmf<-snmf("~/crap.geno",K = npops, alpha = 10, project = "new")
  qmatrix = Q(obj.snmf, K = npops)
  bobuncut=gather(data.frame(t(qmatrix)),dude,qval,1:(nrow(qmatrix))) %>% 
    mutate(pop=factor(rep(1:npops,nrow(qmatrix))))

  #structure on cut data
  obj.snmf<-snmf("~/crapcut.geno",K = npops, alpha = 10, project = "new")
  qmatrix = Q(obj.snmf, K = npops)
  
  #make new matrix with 0s for left out individuals
  if(length(left_out)>0){
    m <- matrix(0, ncol=npops,nrow=5*npops)
    m[!(1:nrow(m) %in% left_out),] <- qmatrix
  } else{ m=qmatrix}
  bobcut=gather(data.frame(t(m)),dude,qval,1:(nrow(m))) %>% 
    mutate(pop=factor(rep(1:npops,nrow(m))))
  
  #graphs
  cut<-ggplot(bobcut,aes(x=dude,y=qval,fill=pop))+
    geom_bar(stat="identity")+
    theme(
      axis.text.y=element_blank(),
      axis.ticks.y=element_blank(),
      axis.title.x=element_blank(),
      axis.title.y=element_blank()
    )+
    scale_fill_discrete(guide=FALSE)+
    scale_x_discrete(labels=1:15)

  
    uncut<-ggplot(bobuncut,aes(x=dude,y=qval,fill=pop))+
    geom_bar(stat="identity")+
    theme(
      axis.text.y=element_blank(),
      axis.ticks.y=element_blank(),
      axis.title.x=element_blank(),
      axis.title.y=element_blank()
    )+
    scale_fill_discrete(guide=FALSE)+
    scale_x_discrete(labels=1:15)

    bob=list(cut,uncut)
        graphlist <- replicate(10, qplot(1,1), simplify = FALSE)

     title <- ggdraw() + draw_label(paste("simulation", i, "\n4Nm =",migration_command), fontface='bold')
     plots<-plot_grid(uncut,cut)
    print( plot_grid(title,plots,ncol=1, rel_heights=c(0.3, 1)))
}