#### FIGURE 9

# Requires the user to have downloaded seed.gl.Rdata to the working directory 
# and to set the working directory

#set working directory so seed.gl.Rdata loads
setwd("D:/workspace/R_analysis/R_ms_distance") #CHANGE ME
p2 <- readRDS("seed_gl.Rdata")

# Requires installation and loading of the following

# install.packages("dartR")
library(dartR)
# install.packages("gganimate")
library(gganimate)
# install.packages("gridExtra")
library(gridExtra)
# install.packages("magick")
library(magick)

# function to create a pcoa animation (based on gl.pcoa.plot)

gl.pcoa.plot2 <- 
  function(glPca,
           x,
           scale = FALSE,
           ellipse = FALSE,
           plevel = 0.95,
           pop.labels = "pop",
           interactive = FALSE,
           as.pop = NULL,
           hadjust = 1.5,
           vadjust = 1,
           xaxis = 1,
           yaxis = 2,
           zaxis = NULL,
           pt.size = 2,
           pt.colors = NULL,
           pt.shapes = NULL,
           label.size = 1,
           axis.label.size = 1.5,
           save2tmp = FALSE,
           verbose = NULL) {
    
    hold_x <- x
    hold_glPca <- glPca
    
    # SET VERBOSITY
    verbose <- gl.check.verbosity(verbose)
    
    # FLAG SCRIPT START
    funname <- match.call()[[1]]
    
    # CHECK DATATYPE
    datatype1 <- "glPca"
    datatype2 <- "SNP"

    # SCRIPT SPECIFIC ERROR CHECKING
    
    if (interactive | !is.null(zaxis)) {
      pkg <- "plotly"
      if (!(requireNamespace(pkg, quietly = TRUE))) {
        stop(
          error(
            "Package ",
            pkg,
            " needed for this function to work. Please install it."
          )
        )
      }
    }
    
    if (datatype1=="list") {
      pkg <- "gganimate"
      if (!(requireNamespace(pkg, quietly = TRUE))) {
        stop(
          error(
            "Package ",
            pkg,
            " needed for this function to work. Please install it."
          )
        )
      }
      pkg <- "tibble"
      if (!(requireNamespace(pkg, quietly = TRUE))) {
        stop(
          error(
            "Package ",
            pkg,
            " needed for this function to work. Please install it."
          )
        )
      }
      x <- x[[1]]
      glPca <- glPca[[1]]
    }
    
    if (pop.labels != "none" &&
        pop.labels != "ind" &&
        pop.labels != "pop" && pop.labels != "legend") {
      cat(
        warn(
          "  Warning: Parameter 'pop.labels' must be one of none|ind|pop|legend, set to 'pop'\n"
        )
      )
      pop.labels <- "pop"
    }
    if (plevel < 0 | plevel > 1) {
      cat(warn(
        "  Warning: Parameter 'plevel' must fall between 0 and 1, set to 0.95\n"
      ))
      plevel <- 0.95
    }
    if (hadjust < 0 | hadjust > 3) {
      cat(warn(
        "  Warning: Parameter 'hadjust' must fall between 0 and 3, set to 1.5\n"
      ))
      hadjust <- 1.5
    }
    if (vadjust < 0 | hadjust > 3) {
      cat(warn(
        "  Warning: Parameter 'vadjust' must fall between 0 and 3, set to 1.5\n"
      ))
      vadjust <- 1.5
    }
    if (xaxis < 1 | xaxis > ncol(glPca$scores)) {
      cat(
        warn(
          "  Warning: X-axis must be specified to lie between 1 and the number of retained dimensions of the ordination",
          ncol(glPca$scores),
          "; set to 1\n"
        )
      )
      xaxis <- 1
    }
    if (yaxis < 1 | yaxis > ncol(glPca$scores)) {
      cat(
        warn(
          "  Warning: Y-axis must be specified to lie between 1 and the number of retained dimensions of the ordination",
          ncol(glPca$scores),
          "; set to 2\n"
        )
      )
      yaxis <- 2
    }
    if (!is.null(zaxis)) {
      if (zaxis < 1 | zaxis > ncol(glPca$scores)) {
        cat(
          warn(
            "  Warning: Z-axis must be specified to lie between 1 and the number of retained dimensions of the ordination",
            ncol(glPca$scores),
            "; set to 3\n"
          )
        )
        zaxis <- 3
      }
    }
    
    # Assign the new population list if as.pop is specified
    if(datatype2 %in% c("SNP","SilicoDArT")){
      pop.hold <- pop(x)
      if (!is.null(as.pop)) {
        if (as.pop %in% names(x@other$ind.metrics)) {
          pop(x) <- as.matrix(x@other$ind.metrics[as.pop])
          if (verbose >= 2) {
            cat(
              report(
                "  Temporarily setting population assignments to",
                as.pop,
                "as specified by the as.pop parameter\n"
              )
            )
          }
        } else {
          stop(
            error(
              "Fatal Error: individual metric assigned to 'pop' does not exist. Check names(gl@other$loc.metrics) and select again\n"
            )
          )
        }
      }
    }
    
    # If an fd object, pull out the genlight object
    if (datatype2 == "fd") {
      x <- x$fd
      datatype2 <- utils.check.datatype(x, verbose = 0)
    }
    axis.label.size <- axis.label.size * 10
    
    # DO THE JOB
    # Set NULL to variables to pass CRAN checks
    gen <- NULL  
    
    if(datatype1=="list"){
      
      gen_number <- length(hold_x)
      df_sim <- as.data.frame(matrix(ncol = 5))
      colnames(df_sim) <- c("PCoAx","PCoAy","ind","pop","gen")
      
      test_pos_neg <- as.data.frame(matrix(nrow = gen_number,ncol = 3 ))
      colnames(test_pos_neg) <- c("gen","test_x","test_y")
      
      # the direction of the PCA axes are chosen at random 
      # this is to set the same direction in every generation
      # first get the individual with more variance for axis x and y 
      # for the first generation of the simulations
      ind_x_axis <- which.max(abs(hold_glPca[[1]]$scores[,xaxis]))
      ind_y_axis <- which.max(abs(hold_glPca[[1]]$scores[,yaxis]))
      
      # check whether is positive or negative
      test_pos_neg[1, "test_x"] <- 
        if(hold_glPca[[1]]$scores[ind_x_axis,xaxis]>=0)"positive"else"negative"
      test_pos_neg[1, "test_y"]  <- 
        if(hold_glPca[[1]]$scores[ind_y_axis,yaxis]>=0)"positive"else"negative"
      for(sim_i in 1:gen_number){
        glPca <- hold_glPca[[sim_i]]
        x <- hold_x[[sim_i]]
        m <- cbind(glPca$scores[, xaxis], glPca$scores[, yaxis])
        df <- data.frame(m)
        # Convert the eigenvalues to percentages
        # s <- sum(glPca$eig[glPca$eig >= 0])
        # e <- round(glPca$eig * 100 / s, 1)
        # Labels for the axes and points
        xlab <- paste("PCA Axis", xaxis)
        ylab <- paste("PCA Axis", yaxis)
        ind <- indNames(x)
        pop <- factor(pop(x))
        gen <- unique(x$other$sim.vars$generation)
        df <- cbind(df, ind, pop,unique(x$other$sim.vars$generation))
        colnames(df) <- c("PCoAx", "PCoAy", "ind", "pop","gen")
        
        test_pos_neg[ sim_i, "test_x"] <- 
          if(hold_glPca[[sim_i]]$scores[ind_x_axis,xaxis]>=0)"positive"else"negative"
        test_pos_neg[ sim_i, "test_y"]  <- 
          if(hold_glPca[[sim_i]]$scores[ind_y_axis,yaxis]>=0)"positive"else"negative"
        
        if(test_pos_neg[1, "test_x"] != test_pos_neg[ sim_i, "test_x"]){
          df$PCoAx <- df$PCoAx * -1
          # test_pos_neg[ sim_i, "test_x"] <- test_pos_neg[ axis_ind-1, "test_x"] 
        }
        
        if(test_pos_neg[ 1, "test_y"] != test_pos_neg[ sim_i, "test_y"]){
          df$PCoAy <- df$PCoAy * -1
          
          # test_pos_neg[ sim_i, "test_y"] <- test_pos_neg[ axis_ind-1, "test_y"] 
        }
        
        df_sim <- rbind(df_sim,df)
      }
      df_sim <- tibble::as_tibble(df_sim)
      df_sim <- df_sim[-1,]
      
      p  <- ggplot(df_sim, aes(PCoAx, PCoAy, colour = pop)) +
        geom_point(size=3) +
        labs(title = 'Generation: {frame_time}', x = xlab, y = ylab) +
        gganimate::transition_time(gen) +
        gganimate::ease_aes('linear')
      return(p)
    }
    
    PCoAx <- PCoAy <- NULL
    
    # Create a dataframe to hold the required scores
    if (is.null(zaxis)) {
      m <- cbind(glPca$scores[, xaxis], glPca$scores[, yaxis])
    } else {
      m <-
        cbind(glPca$scores[, xaxis], glPca$scores[, yaxis], glPca$scores[, zaxis])
    }
    df <- data.frame(m)
    
    # Convert the eigenvalues to percentages
    s <- sum(glPca$eig[glPca$eig >= 0])
    e <- round(glPca$eig * 100 / s, 1)
    
    # Labels for the axes and points
    
    if (datatype2 == "SNP" | datatype2 == "SilicoDArT") {
      xlab <- paste("PCA Axis", xaxis, "(", e[xaxis], "%)")
      ylab <- paste("PCA Axis", yaxis, "(", e[yaxis], "%)")
      if (!is.null(zaxis)) {
        zlab <- paste("PCA Axis", zaxis, "(", e[zaxis], "%)")
      }
      
      ind <- indNames(x)
      pop <- factor(pop(x))
      df <- cbind(df, ind, pop)
      if (is.null(zaxis)) {
        colnames(df) <- c("PCoAx", "PCoAy", "ind", "pop")
      } else {
        colnames(df) <- c("PCoAx", "PCoAy", "PCoAz", "ind", "pop")
      }
      
    } else {
      # datatype2 == 'dist'
      xlab <- paste("PCoA Axis", xaxis, "(", e[xaxis], "%)")
      ylab <- paste("PCoA Axis", yaxis, "(", e[yaxis], "%)")
      if (!is.null(zaxis)) {
        zlab <- paste("PCA Axis", zaxis, "(", e[zaxis], "%)")
      }
      
      ind <- rownames(as.matrix(x))
      pop <- ind
      df <- cbind(df, ind, pop)
      if (is.null(zaxis)) {
        colnames(df) <- c("PCoAx", "PCoAy", "ind", "pop")
      } else {
        colnames(df) <- c("PCoAx", "PCoAy", "PCoAz", "ind", "pop")
      }
      if (interactive) {
        cat(
          warn(
            "  Sorry, interactive labels are not available for an ordination generated from a Distance Matrix\n"
          )
        )
        cat(warn(
          "  Labelling the plot with names taken from the Distance Matrix\n"
        ))
      }
      pop.labels <- "pop"
    }
    
    ####### 2D PLOT
    if (is.null(zaxis)) {
      # If population labels
      
      if (pop.labels == "pop") {
        if (datatype2 == "SNP") {
          if (verbose >= 2)
            cat(report(
              "  Plotting populations in a space defined by the SNPs\n"
            ))
        } else if (datatype2 == "SilicoDArT") {
          if (verbose >= 2)
            cat(
              report(
                "  Plotting populations in a space defined by the presence/absence data\n"
              )
            )
        } else {
          if (verbose >= 2)
            cat(report("  Plotting entities from the Distance Matrix\n"))
        }
        
        # Plot
        if (is.null(pt.shapes)) {
          plott <-
            ggplot(df,
                   aes(
                     x = PCoAx,
                     y = PCoAy,
                     group = pop,
                     color = pop
                   ))
        } else {
          plott <-
            ggplot(df,
                   aes(
                     x = PCoAx,
                     y = PCoAy,
                     group = pop,
                     color = pop,
                     shape = pop
                   ))
        }
        plott <- plott + geom_point(size = pt.size, aes(color = pop)) + 
          directlabels::geom_dl(aes(label = pop),
                                method = list("smart.grid", 
                                              cex = label.size)) + 
          theme(axis.title = element_text(face = "bold.italic", 
                                          size = axis.label.size,
                                          color = "black"),
                axis.text.x = element_text(face = "bold", 
                                           angle = 0,
                                           vjust = 0.5,
                                           size = axis.label.size),
                axis.text.y = element_text(face = "bold", 
                                           angle = 0,
                                           vjust = 0.5,
                                           size = axis.label.size)) +
          labs(x = xlab, y = ylab)
        
        if (!is.null(pt.shapes)) {
          plott <- plott + scale_shape_manual(values = pt.shapes)
        }
        if (!is.null(pt.colors)) {
          plott <- plott + scale_color_manual(values = pt.colors)
        }
        plott <-
          plott + geom_hline(yintercept = 0) + 
          geom_vline(xintercept = 0) + 
          theme(legend.position = "none")
        # Scale the axes in proportion to % explained, if requested if(scale==TRUE) { plott <- plott +
        # coord_fixed(ratio=e[yaxis]/e[xaxis]) }
        if (scale == TRUE) {
          plott <- plott + coord_fixed(ratio = e[yaxis]/ e[xaxis])
        }
        # Add ellipses if requested
        if (ellipse == TRUE) {
          plott <- plott + stat_ellipse(type = "norm", level = plevel)
        }
      }
      
      # If interactive labels
      
      if (interactive) {
        cat(report("  Displaying an interactive plot\n"))
        cat(
          warn(
            "  NOTE: Returning the ordination scores, not a ggplot2 compatable object\n"
          )
        )
        
        # Plot
        plott <-
          ggplot(df, aes(
            x = PCoAx,
            y = PCoAy,
            label = ind
          )) + geom_point(size = pt.size, aes(color = pop)) + theme(
            axis.title = element_text(
              face = "bold.italic",
              size = axis.label.size,
              color = "black"
            ),
            axis.text.x = element_text(
              face = "bold",
              angle = 0,
              vjust = 0.5,
              size = axis.label.size
            ),
            axis.text.y = element_text(
              face = "bold",
              angle = 0,
              vjust = 0.5,
              size = axis.label.size
            ),
            legend.title = element_text(
              color = "black",
              size = axis.label.size,
              face = "bold"
            ),
            legend.text = element_text(
              color = "black",
              size = axis.label.size,
              face = "bold"
            )
          ) +
          labs(x = xlab, y = ylab) + geom_hline(yintercept = 0) + geom_vline(xintercept = 0) + theme(legend.position = "none")
        # Scale the axes in proportion to % explained, if requested if(scale==TRUE) { plott <- plott +
        # coord_fixed(ratio=e[yaxis]/e[xaxis]) }
        if (scale == TRUE) {
          plott <- plott + coord_fixed(ratio = 1)
        }
        # Add ellipses if requested
        if (ellipse == TRUE) {
          plott <-
            plott + stat_ellipse(aes(color = pop),
                                 type = "norm",
                                 level = plevel)
        }
        cat(warn(
          "  Ignore any warning on the number of shape categories\n"
        ))
      }
      
      # If labels = legend
      
      if (pop.labels == "legend") {
        if (verbose >= 2)
          cat(report("  Plotting populations identified by a legend\n"))
        
        # Plot
        Population <- pop
        if (is.null(pt.shapes)) {
          plott <-
            ggplot(df,
                   aes(
                     x = PCoAx,
                     y = PCoAy,
                     group = Population,
                     color = Population
                   ))
        } else {
          plott <-
            ggplot(
              df,
              aes(
                x = PCoAx,
                y = PCoAy,
                group = pop,
                color = Population,
                shape = Population
              )
            )
        }
        plott <-
          plott + geom_point(size = pt.size, aes(color = pop)) + theme(
            axis.title = element_text(
              face = "bold.italic",
              size = axis.label.size,
              color = "black"
            ),
            axis.text.x = element_text(
              face = "bold",
              angle = 0,
              vjust = 0.5,
              size = axis.label.size
            ),
            axis.text.y = element_text(
              face = "bold",
              angle = 0,
              vjust = 0.5,
              size = axis.label.size
            ),
            legend.title = element_text(
              color = "black",
              size = axis.label.size,
              face = "bold"
            ),
            legend.text = element_text(
              color = "black",
              size = axis.label.size,
              face = "bold"
            )
          ) + labs(x = xlab, y = ylab)
        if (!is.null(pt.shapes)) {
          plott <- plott + scale_shape_manual(values = pt.shapes)
        }
        if (!is.null(pt.colors)) {
          plott <- plott + scale_color_manual(values = pt.colors)
        }
        plott <-
          plott + geom_hline(yintercept = 0) + geom_vline(xintercept = 0)
        # Scale the axes in proportion to % explained, if requested if(scale==TRUE) { plott <- plott +
        # coord_fixed(ratio=e[yaxis]/e[xaxis]) }
        if (scale == TRUE) {
          plott <- plott + coord_fixed(ratio = 1)
        }
        # Add ellipses if requested
        if (ellipse == TRUE) {
          plott <- plott + stat_ellipse(type = "norm", level = plevel)
        }
      }
      
      # If labels = none
      
      if (pop.labels == "none" | pop.labels == FALSE) {
        if (verbose >= 0)
          cat(report("  Plotting points with no labels\n"))
        
        # Plot
        if (is.null(pt.shapes)) {
          plott <- ggplot(df, aes(
            x = PCoAx,
            y = PCoAy,
            color = pop
          ))
        } else {
          plott <-
            ggplot(df,
                   aes(
                     x = PCoAx,
                     y = PCoAy,
                     color = pop,
                     shape = pop
                   ))
        }
        plott <-
          plott + geom_point(size = pt.size, aes(color = pop)) + theme(
            axis.title = element_text(
              face = "bold.italic",
              size = axis.label.size,
              color = "black"
            ),
            axis.text.x = element_text(
              face = "bold",
              angle = 0,
              vjust = 0.5,
              size = axis.label.size
            ),
            axis.text.y = element_text(
              face = "bold",
              angle = 0,
              vjust = 0.5,
              size = axis.label.size
            )
          ) + labs(x = xlab, y = ylab)
        if (!is.null(pt.shapes)) {
          plott <- plott + scale_shape_manual(values = pt.shapes)
        }
        if (!is.null(pt.colors)) {
          plott <- plott + scale_color_manual(values = pt.colors)
        }
        plott <-
          plott + geom_hline(yintercept = 0) + geom_vline(xintercept = 0) + theme(legend.position = "none")
        # Scale the axes in proportion to % explained, if requested if(scale==TRUE) { plott <- plott +
        # coord_fixed(ratio=e[yaxis]/e[xaxis]) }
        if (scale == TRUE) {
          plott <- plott + coord_fixed(ratio = 1)
        }
        # Add ellipses if requested
        if (ellipse == TRUE) {
          plott <- plott + stat_ellipse(type = "norm", level = plevel)
        }
      }
      
      if (verbose >= 2) {
        cat(report("  Preparing plot .... please wait\n"))
      }
      if (interactive) {
        plott <- plotly::ggplotly(plott)
        show(plott)
      } else {
        show(plott)
      }
    }  # End 2D plot
    
    ##### IF 3D PLOT
    if (!is.null(zaxis)) {
      if (verbose >= 2) {
        cat(
          report(
            "  Displaying a three dimensional plot, mouse over for details for each point\n"
          )
        )
      }
      plott <-
        plotly::plot_ly(
          df,
          x = ~ PCoAx,
          y = ~ PCoAy,
          z = ~ PCoAz,
          marker = list(size = pt.size * 2),
          colors = pt.colors,
          text = ind
        ) %>%
        plotly::add_markers(color = ~ pop) %>%
        plotly::layout(
          legend = list(title = list(text = "Populations")),
          scene = list(
            xaxis = list(
              title = xlab,
              titlefont = list(size = axis.label.size / 2)
            ),
            yaxis = list(
              title = ylab,
              titlefont = list(size = axis.label.size / 2)
            ),
            zaxis = list(
              title = zlab,
              titlefont = list(size = axis.label.size / 2)
            )
          )
        )
      show(plott)
      if (verbose >= 2) {
        cat(warn("  May need to zoom out to place 3D plot within bounds\n"))
      }
    }
    
    # creating temp file names
    if (save2tmp) {
      temp_plot <- tempfile(pattern = "Plot_")
      match_call <-
        paste0(names(match.call()),
               "_",
               as.character(match.call()),
               collapse = "_")
      # saving to tempdir
      saveRDS(list(match_call, plott), file = temp_plot)
      if (verbose >= 2) {
        cat(report("  Saving the ggplot to the session tempfile\n"))
      }
      temp_table <- tempfile(pattern = "Table_")
      saveRDS(list(match_call, df), file = temp_table)
      if (verbose >= 2) {
        cat(report("  Saving tabulation to the session tempfile\n"))
      }
    }
    # FLAG SCRIPT END
    
    if (verbose >= 1) {
      cat(report("Completed:", funname, "\n"))
    }
    
    invisible(plott)
  }


############################### start of script

### Create offsprings between randomly selected parents from one population

#create a sample data set of one populations
pc <- gl.pcoa(p2)

set.seed(333)  #comment for testing other pairs (this creates the figures in the paper)

#sample two parents from "Emmac_MDB_nth"
p12 <- sample(1:nInd(p2), 2, replace = FALSE)


#rename parents so they appear in the plot
pnew <- as.character(pop(p2))
pnew[p12[1]] <- "father"
pnew[p12[2]] <- "mother"
pop(p2)<- factor(pnew)



pca <- gl.pcoa.plot2(pc,p2, verbose = 0)


#run a loop and create 2-8 offsprings (takes some minutes)
pca_off<- list()
pca_gl <- list()
poff <- list()
for (noff in 2:8) { 
  ppar<- p2[p12,]
  offspring <- gl.sim.offspring(ppar[1,], ppar[2,], noffpermother = noff)
  pop(offspring)<- rep("offspring", noff)
  p2_off <- rbind(p2, offspring)
  
  p3 <- new("dartR", as.matrix(p2_off))
  pop(p3)<- pop(p2_off)
  p3 <- gl.compliance.check(p3)
  pca_off[[noff]] <-gl.pcoa(p3)
  pca_gl[[noff]] <- p3
  poff[[noff]] <- p3
}


#name the generations for sim.vars (necessary for the animation)
gls <- poff
gls[[1]]<- p2
for ( i in 1:length(gls)) {
  gls[[i]]@other$sim.vars<- data.frame(generation=i)
  
}

#create the animation
pca_off[[1]]<- pc
anim <- gl.pcoa.plot(glPca=pca_off[1:8],x=gls)
# playing the animation
animate(anim,duration=10, renderer=magick_renderer())


# uncomment to save the animation
#anim_save("pca_offspring_simulation.gif")

#plot single images #1
gl.pcoa.plot(glPca=pca_off[[1]], gls[[1]])

#plot single images #8
gl.pcoa.plot(glPca=pca_off[[8]], gls[[8]])