R Notebook

library("LEA")
library("cluster")
library("cowplot")
library("tidyverse") 
library("viridis")
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"

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

#mspms 41 1 -t 4500 -r 4500 1000000 -I 5 10 10 10 10 1 -m 3 4 0.2 -n 2 0.2 -n 3 0.2  -ej 0.0083 3 2 -en 0.0083  2 0.2 -ej 0.025 2 1  -en 0.025 1 1  -ej 0.1 4 1 -en 0.1 1 1 -ej 2 5 1
#mspms 41 1 -t 4500 -r 4500 1000000 -I 5 10 10 10 10 10 -m 3 4 0.2 -n 2 0.2 -n 3 0.2 -ej 0.0083 3 2 -en 0.0083 -em 0.0083 3 4 0 -ej 0.025 2 1  -en 0.025 1 1  -ej 0.1 4 1 -en 0.1 1 1 -ej 2 5 1

npops=5
ninds=20
kpops=4
for(j in 1:10){
migration_command=paste("-m 3 4 0.2 -n 2 0.2 -n 3 0.2 -ej 0.0083 3 2 -en 0.0083 2 0.2 -em 0.0083 3 4 0 -ej 0.025 2 1  -en 0.025 1 1  -ej 0.1 4 1 -en 0.1 1 1 -ej 2 5 1")
#migration_command=paste("-ej 2 2 1 -ej 2 3 1 --random-seeds 1336385744 876978214 1278603625")  
mydat<-gendata(migration_command,npops,ninds)

#cut based on silhouette scores
keep<-sort(which(silhouette(kmeans(dist(mydat),npops)[[1]],dist(mydat))[,3]>0))

left_out<-which(!(1:(npops*ninds/2) %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="")

allele_freqs=list()
for(i in 1:npops){
  allele_freqs[[i]]=colSums(mydat[1:ninds/2,])/20
}


#structure on original
obj.snmf<-snmf("~/crap.geno",K = kpops, alpha = 10, project = "new") 
qmatrix = Q(obj.snmf, K = kpops)
bobuncut=gather(data.frame(t(qmatrix)),dude,qval,1:(nrow(qmatrix))) %>%
mutate(pop=factor(rep(1:kpops,nrow(qmatrix))),dude=factor(as.numeric(gsub("X","",dude))))

#structure on cut data
obj.snmf<-snmf("~/crapcut.geno",K = kpops, alpha = 10, project = "new") 
qmatrix = Q(obj.snmf, K = kpops)

#make new matrix with 0s for left out individuals
if(length(left_out)>0){
  m <- matrix(0, ncol=kpops,nrow=(npops*ninds/2)) 
  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:kpops,nrow(m))),dude=factor(as.numeric(gsub("X","",dude))))

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_manual(guide=FALSE,values=viridis(npops)) 

    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_manual(guide=FALSE,values=viridis(npops))+
    scale_x_discrete(labels=1:(npops*ninds/2))

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

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