```{r include=FALSE}
library(tidyverse)
library(ggforce)
library(shiny)
library(gridExtra)
library(scales)

generate_parent_matrix <- function(pop_size, ploidy, generations, force_coal=FALSE) {
  
  # Creates a matrix with values describing the parent of each individual
  # The parent is chosen at random from the previous generation (equal probability)
  #
  # Inputs:
  #   pop_size: int for population size
  #   ploidy: int for number of paired individuals
  #   generations: int for total number of generations
  #
  # Returns the matrix described above
  #
  # **Should this be long or wide? Wide fits better with the current coalescent diagram, but other matrix in generate_tree() needs to be long**
  # **This is where we can implement different mating models, currently simplest version with equal probability of any parent**
  
  parent_matrix <- matrix(NA, nrow=(pop_size*ploidy), ncol=generations)
  for (gen in ncol(parent_matrix):2) {
    if (gen==2 & force_coal==TRUE) {
      for (indiv in 1:nrow(parent_matrix)) {
        parent_matrix[indiv, gen] <- nrow(parent_matrix) %/% 2
      }
    } else {
      for (indiv in 1:nrow(parent_matrix)) {
        parent_matrix[indiv, gen] <- sample(1:nrow(parent_matrix), 1)+(gen-2)*nrow(parent_matrix)
      }
    }
  }
  return(parent_matrix)
}

scatter_mutations <- function(parent_matrix, mutation_rate, starting_pop=c()) {
  
  # Mutates individuals (changes allele number), following Infinite Sites and Infinite Alleles Assumptions
  # All individuals start with the same conditions
  #
  # Inputs:
  #   parent_matrix: matrix describing parentage of each individual
  #   mutation_rate: float between 0 and 1 describing the probability of mutation during reproduction
  #
  # Returns data.frame with the mutation value for each individual in the matrix
  mutation_counter <- 1
  rows <- nrow(parent_matrix)
  cols <- ncol(parent_matrix)
  if (length(starting_pop) == 0) {
    starting_pop <- rep(0, rows)
  }
  poly_matrix <- matrix(starting_pop, nrow=rows, ncol=cols)
  for (gen in 1:cols) {
    for (i in 1:rows) {
      if (gen != 1) {
        if (runif(1) <= mutation_rate) {
          poly_matrix[i, gen] <- mutation_counter
          mutation_counter <- mutation_counter + 1
        } else {
          poly_matrix[i, gen] <- poly_matrix[parent_matrix[i, gen]]
        }
      }
    }
  }
  poly_dataframe <- data.frame("indiv"=1:(cols*rows), "allele"=as.factor(poly_matrix[1:(cols*rows)]))
  return(poly_dataframe)
}

find_lineage <- function(indiv, parent_matrix) {
  
  # Recursively determines the lineages of an individual based on the parent_matrix provided
  #
  # Inputs:
  #   indiv: int for individual id
  #   parent_matrix: matrix of parents
  #
  # Returns a vector containing the individual ids in the associated lineage
  
  if (is.na(parent_matrix[indiv])) {
    return(indiv)
  } else {
    return(c(indiv, find_lineage(indiv=parent_matrix[indiv], parent_matrix=parent_matrix)))
  }
}

flatten_lineage_list <- function(indiv, lineage_list) {
  
  # Converts list element for the lineage of an individual to a single data.frame that can be more easily handled
  #
  # Inputs:
  #   indiv: int for individual id
  #   lineage_list: list of all lineages, output from lapply of find_lineage()
  #
  # Returns a data.frame with two columns:
  #   starting_indiv: individual id that the lineage is associated with
  #   ancestor: ancestor's id
  # Useful to group by starting_indiv for plotting and other analyses
  
  
  lineage_dataframe <- data.frame(as.factor(indiv), lineage_list[indiv])
  names(lineage_dataframe) <- c("starting_indiv", "ancestor")
  return(lineage_dataframe)
}

locate_coalescent_events <- function(lineages, tracked, pop_size, ploidy, generations) {
  
  # Locates the coalescent events by identifying the most recent point of shared ancestry between tracked lineages
  #
  # Inputs:
  #   lineages: data.frame of flattened lineages
  #   tracked: vector of tracked lineages ids
  #   pop_size: int for population size
  #   ploidy: int for number of paired chromosomes
  #   generations: int for total number of generations
  #
  # Returns a data.frame with one row per coalescent event and columns:
  #   ancestor: ancestor's id
  #   number_of_tracked_lineages: number of tracked lineages that shared this ancestor
  #   tracked_lineages: ids of tracked lineages which share this ancestor, "_" delimited
  #   generations_in_past: number of generations in past that coalescent event occurred, this assumes:
  #     - all tracked lineages are from the most recent (rightmost) generation
  #     - all tracked lineages are from the same generation (somewhat redundant to the first assumption, but also important)
  #
  # **May want to remove some of these parameters (pop_size, ploidy, generations) in the future, if there is an elegant way around them**
  
  shared_ancestry <- lineages %>%
    filter(starting_indiv %in% tracked) %>%
    group_by(ancestor) %>%
    summarize(
      number_of_tracked_lineages=n(),
      tracked_lineages=paste(starting_indiv, collapse=" ")
    ) %>%
    filter(number_of_tracked_lineages > 1)
  
  if (nrow(shared_ancestry) > 0) {
    return(
      as.data.frame(
        shared_ancestry %>%
          arrange(number_of_tracked_lineages) %>%
          mutate(generations_in_past = (generations-1) - ((ancestor-1) %/% (pop_size*ploidy))) %>%
          group_by(tracked_lineages) %>%
          filter(generations_in_past == min(generations_in_past))
      )
    )
  } else {
    return(
      data.frame(
        ancestor=character(),
        number_of_tracked_lineages=numeric(),
        tracked_lineages=character(),
        generations_in_past=numeric()
      )
    )
  }
}

identify_deepest_ancestor <- function(indiv, flat_lineages) {
  
  # Locates the oldest ancestor of a provided lineage
  # Useful for knowing how many trees there are in the tracked lineages and which lineages belong to each tree
  #
  # Inputs:
  #   indiv: int for individual id
  #   flat_lineages: data.frame returned from flatten_lineage_list
  #
  # Returns ancestor's id
  
  return(min(flat_lineages[which(flat_lineages$starting_indiv==indiv),]$ancestor))
}

order_lineages <- function(coalesced_lineages, tracked) {
  
  # Calculates an order of lineages that reduces crossovers in tree diagram
  #
  # Inputs:
  #   coalesced_lineages: vector of coalesced lineages from coalescent events data.frame
  #   tracked: vector of all tracked lineages, even if they don't coalesce
  #
  # Returns a vector of the ordered lineages
  
  if (length(coalesced_lineages) > 0) {
    lineages_order <- as.numeric(unlist(strsplit(coalesced_lineages[1], " ")))
    if (length(coalesced_lineages) > 1) {
      for (event in 2:length(coalesced_lineages)) {
        inserting <- as.numeric(unlist(strsplit(coalesced_lineages[event], " ")))
        reordered_inserting <- inserting[order(match(inserting, lineages_order))]
        if (!all(is.na(match(inserting, lineages_order)))) {
          lineages_order <- unique(c(lineages_order[1:min(match(reordered_inserting, lineages_order), na.rm=TRUE)], reordered_inserting, lineages_order[min(match(reordered_inserting, lineages_order), na.rm=TRUE):length(lineages_order)])) 
        } else {
          lineages_order <- c(lineages_order, reordered_inserting)
        }
      }
    }
    lineages_order <- unique(c(lineages_order, tracked)) 
  } else {
    lineages_order <- tracked
  }
  return(lineages_order)
}

generate_tree <- function(coalescent_events, tracked, generations) {
  
  # Calculates values for a coalescent tree based on previously located coalescent events
  # These values are completely relative and are not based on any biologically relevant scale
  #
  # Inputs:
  #   coalescent_events: data.frame containing data on the coalescent events
  #   tracked: vector of all tracked lineages, even if they don't coalesce
  #   generations: int for total number of generations
  #
  # Returns a data.frame that can be plotted to show the coalescent tree
  
  if (length(tracked) > 0) {
    lineages_order <- order_lineages(coalesced_lineages=coalescent_events$tracked_lineages, tracked=tracked)
    coalescent_tree_matrix <- matrix(c(1:length(lineages_order)), nrow=generations, ncol=length(lineages_order), byrow=TRUE)
    colnames(coalescent_tree_matrix) <- lineages_order
    if (nrow(coalescent_events) > 0) {
      for (event in 1:nrow(coalescent_events)) {
        coalescent_tree_matrix[1:(generations-coalescent_events$generations_in_past[event]),c(unlist(strsplit(coalescent_events$tracked_lineages[event], " ")))] <- mean(coalescent_tree_matrix[(generations-coalescent_events$generations_in_past[event]),c(unlist(strsplit(coalescent_events$tracked_lineages[event], " ")))])
      } 
    }
    coalescent_tree_dataframe <- data.frame("generation"=1:generations, as.data.frame(coalescent_tree_matrix)) %>%
      gather(lineage, value, -c(generation))
    coalescent_tree_dataframe$lineage <- as.numeric(gsub("X", "", coalescent_tree_dataframe$lineage))
    return(coalescent_tree_dataframe) 
  } else {
    return(data.frame("generation"=numeric(), "lineage"=numeric(), "value"=numeric()))
  }
}

locate_mutation_events <- function(allele, polymorphisms) {
  
  # Identifies the earliest instance of each mutation
  #
  # Inputs:
  #   allele: int for the allele of interest
  #   polymorphisms: data.frame containing polymorphism data for each individual
  #
  # Returns the individual id for the earliest individual with the allele of interest
  
  indiv <- min(polymorphisms[which(polymorphisms$allele==allele),]$indiv)
  return(indiv)
}

identify_allele_seq <- function(event, tracked, tree) {
  
  # Determines which allele of a mutated SNP each individual has based on location of mutation on tree
  # Unknown alleles, marked with "?", occur when an individual does not coalesce with the tree containing the mutation
  #
  # Inputs:
  #   event: mutation event from tracked_mutation_events row
  #   tracked: vector of all tracked lineages, even if they don't coalesce
  #   tree: data.frame with individuals and their associated tree
  #
  # Returns vector with each individual's genotype for the mutation event of interest
  
  allele <- rep("?", length(tracked))
  names(allele) <- tracked
  lineages_on_mutated_tree <- tree$lineage #[which(tree$tree==event["tree"]),]
  allele[which(tracked %in% lineages_on_mutated_tree)] <- event["ancestral_allele"]
  allele[which(tracked %in% as.numeric(unlist(strsplit(event["lineages"], split=" "))))] <- event["derived_allele"]
  return(allele)
}

select_alternate_ancestor <- function(indiv, orig_matrix) {
  if (!is.na(orig_matrix[indiv])) {
    if (sample(0:1, 1) == 1) {
      if ((orig_matrix[indiv]+1) %% 2 == 0) {
        return(orig_matrix[indiv]+1)
      } else {
        return(orig_matrix[indiv]-1)
      }   
    } else {
      return(orig_matrix[indiv])
    }
  } else {
    return(NA)
  }
}

show_pedigree_vs_tree <- function(plot_pedigree=FALSE) {
  flat_lineages_alt <- reactive({
    alt_matrix <- matrix(sapply(1:length(values$parent_matrix), select_alternate_ancestor, orig_matrix=values$parent_matrix), nrow=(pop_size_app1*ploidy_app1), ncol=generations_app1)
    full_lineages_alt <- lapply(1:(pop_size_app1*ploidy_app1*generations_app1), find_lineage, parent_matrix=alt_matrix)
    flat_lineages_alt <- do.call(rbind, lapply(1:length(full_lineages_alt), flatten_lineage_list, lineage_list=full_lineages_alt))
    flat_lineages_alt$pedigree <- (flat_lineages_alt$ancestor+(ploidy_app1-1)) %/% ploidy_app1
    flat_lineages_alt
  })
  
  
  coalescent_tree_alt <- reactive({
    coalescent_events_alt <- locate_coalescent_events(lineages=flat_lineages_alt(), tracked=tracked_app1, pop_size=pop_size_app1, ploidy=ploidy_app1, generations=generations_app1)
    generate_tree(coalescent_events=coalescent_events_alt, tracked=tracked_app1, generations=generations_app1)
  })
  
  labels_alt <- reactive({
    coalescent_tree_alt()[which(coalescent_tree_alt()$generation==generations_app1),]
  })
  
  if (plot_pedigree==TRUE) {
    output$pedigree <- renderPlot({
    ggplot() +
      theme_void() +
      theme(legend.position="none") +
      coord_fixed() +
      geom_line(
        data=flat_lineages_alt(),
        aes(
          (pedigree-1) %/% (pop_size_app1),
          (pop_size_app1) - ((pedigree-1) %% (pop_size_app1)),
          group=starting_indiv
        ),
        size=1
      ) +
      geom_circle(
        data=data.frame("point"=1:(pop_size_app1*generations_app1)),
        aes(
          x0=(point-1) %/% (pop_size_app1),
          y0=(pop_size_app1) - ((point-1) %% (pop_size_app1)),
          r=0.25
        ),
        size=1,
        color="black",
        fill="white"
      )
  }) 
  }
  
  output$coalescent <- renderPlot({
    ggplot() +
    theme_void() +
    theme(legend.position="none") +
    coord_fixed() +
    geom_circle(
      data=data.frame("point"=1:(pop_size_app1*generations_app1)),
      aes(
        x0=(point-1) %/% (pop_size_app1),
        y0=(pop_size_app1) - ((point-1) %% (pop_size_app1)),
        r=0.25
      ),
      size=1,
      color="black",
      fill="white"
    ) +
    geom_point(
      data=data.frame("point"=1:(pop_size_app1*ploidy_app1*generations_app1)),
      aes(
        (point-1) %/% (pop_size_app1*ploidy_app1),
        (pop_size_app1) - (((point-1) %/% ploidy_app1) %% (pop_size_app1)) - ((point-1) %% ploidy_app1 - (ploidy_app1-1)/2)/5
      ),
      size=2,
      color="grey"
    ) +
    geom_text(
      data=data.frame("point"=((pop_size_app1*ploidy_app1*generations_app1-pop_size_app1*ploidy_app1)+1):(pop_size_app1*ploidy_app1*generations_app1)),
      aes(
        (point-1) %/% (pop_size_app1*ploidy_app1) + generations_app1/20,
        (pop_size_app1) - (((point-1) %/% ploidy_app1) %% (pop_size_app1)) - ((point-1) %% ploidy_app1 - (ploidy_app1-1)/2)/5,
        label=point-(pop_size_app1*ploidy_app1*generations_app1-pop_size_app1*ploidy_app1)
      )
    ) +
    geom_line(
      data=flat_lineages_alt(),
      aes(
        (ancestor-1) %/% (pop_size_app1*ploidy_app1),
        (pop_size_app1) - (((ancestor-1) %/% ploidy_app1) %% (pop_size_app1)) - ((ancestor-1) %% ploidy_app1 - (ploidy_app1-1)/2)/5,
        group=starting_indiv
      ),
      size=1,
      color="grey"
    ) +
    geom_line(
      data=flat_lineages_alt()[which(flat_lineages_alt()$starting_indiv %in% tracked_app1),],
      aes(
        (ancestor-1) %/% (pop_size_app1*ploidy_app1),
        (pop_size_app1) - (((ancestor-1) %/% ploidy_app1) %% (pop_size_app1)) - ((ancestor-1) %% ploidy_app1 - (ploidy_app1-1)/2)/5,
        group=starting_indiv
      ),
      size=1,
      color="#4285f4ff"
    )
  })
  
  output$tree <- renderPlot({
    ggplot() +
      theme_void() +
      theme(legend.position="none") +
      geom_line(data=coalescent_tree_alt(), aes(generation, value, group=lineage), size=1) +
      geom_text(data=labels_alt(), aes(generation+generations_app1/20, value, label=lineage-(pop_size_app1*ploidy_app1*generations_app1-pop_size_app1*ploidy_app1)))
  })
}
```

<br>

In evolutionary biology, we often think about the phylogeny of a group, the underlying relationships between the samples. When the comparison is of species (particularly those that are distantly related), it is common for these relationships to be relatively well-defined, allowing researchers to represent the shared history as a single phylogenetic tree of life. As we compare groups with greater relatedness, such as looking at a random sample of humans across the globe, the genetic relationships are no longer clear-cut. In this post, I take a step back and highlight the origins and consequences of Mendelian segregation in the context of modern population genetic models.

In an early general biology class, I was introduced to Gregor Mendel’s work from the 1880’s, where he performed crosses between varieties of pea plants and meticulously studied the inheritance of traits over generations. I was hooked by the idea of relatively simple patterns that would allow us to predict the frequency of traits in subsequent generations and remember coming home and working through my own “calculations” of his results (lots of Punnett squares). In support of a child with a newfound interest in genetics, my parents bought the family 23andMe DNA testing kits. We sent out three samples (my mom, my sister, and I) and read through our results together. Despite how closely related we are, there are clear differences across our genomes.

<br>

<center>![<br><br>*A comparison between my sister’s genome and my own generated by 23andMe.*](www/23andMe_SelenaJamesComparison.png)</center>

<br>

Between my sister and me, there are regions that we have in common (completely or half identically) and other regions we don’t share. Sequences that are marked as “Completely Identical” represent regions that we both inherited from my mother and my father; we would say that we share two common genetic ancestors (one for each chromosome) in the previous generation. Here, it is important to separate general ancestry from genetic ancestry; whereas ancestry is often used to refer to anyone above you in your family tree, your genetic ancestors are the subset of those individuals from whom you’ve inherited regions of your genome. Continuing further, sequences that are marked as “Half Identical” represent regions where we share a single common genetic ancestor in the previous generation; from this diagram, we cannot discern whether that ancestor is from our mother or father. Lastly, sequences that are marked as “Not Identical” represent regions that require us to look much deeper into the past to identify a common genetic ancestor. Each locus, or position, in our genome reflects a pattern of inheritance that extends back, from generation to generation, to the beginning of life.

By unraveling patterns of sharing between genomes in a population back over multiple generations, researchers can learn about important factors that affect genetic inheritance, such as migration, population size, and natural selection. A lot has been discovered with regards to patterns of inheritance in the over one hundred years since Mendel. Here, I highlight the principles behind coalescent theory and explain the consequences of the random process of genetic inheritance that we observe in our genomes as it plays out over multiple generations.

In 1982, Kingman published “The Coalescent”, describing coalescent theory, a model of inheritance where ancestral lineages in a population are tracked backwards through time, a process independently described by a number of researchers at that time; see [Kingman (2000)](https://watermark.silverchair.com/genetics1461.pdf?token=AQECAHi208BE49Ooan9kkhW_Ercy7Dm3ZL_9Cf3qfKAc485ysgAAAtwwggLYBgkqhkiG9w0BBwagggLJMIICxQIBADCCAr4GCSqGSIb3DQEHATAeBglghkgBZQMEAS4wEQQMQa5opcvr3Ua50KeSAgEQgIICj9O3zN1GSl_jlemxMejM0TaHx1MJsGOR6_keQEKDuZpEp4Ii8Wtx4y7TVrfEx_cPW8zICEqgBNVm8Lf9wHqPdLaK0d6LHjvCwLjDHchw0PxEZnF1LcFy-VzilRvGFIudbTsDvCMj-XYeBKzUeu7ee5FpzBmQ0gAe37uLWpwfQH2xrWj-GutAEOKb4c-spui8x-D6mh88fGxh3Y2WFBsJLnUN8S5i96EegFU0e9xjESDeRwjfMd3FwMuSJ0lLal74jEWGD1AParr8MT3_9gJ_nnso9GG1O91-WrLaaupbk2CqOAraMJstHCMivxAwSNoMZWimgh5hEcAT3DryJg_cwebemQdFu7vMJYjHe_2uriW35Uqe7IZxGWw-T7UxyH36ywFpCghOwBTGTSAY9pDf73YBG77ATlT3VuLvdNBOA9eU9uwAEMD8Pe2WpkLXUmxvVAn539WJ4VLgxT7YKF8XcZuTLvBp9c0Tu0OSlHw02ZqLNs9S6FCfYroOE-r6wguPlY0MgExkqR2Tf9vCxbyDIhIectNG8FubLew9V_Ipg8nChMdF6sNA_OCwJZiA6RboxCHdYziT3FDUFbcGeWSH98RvAWI2-kWiBIyfgkVlQsQOFWahy7398S-cfLwSvgQYV01RjWME1T6UQYZ_W-x4zumQXKpqGSrxQwQKb5py26xMbPz1PiHg0UVi-8zyXBedw9AobhEKOC_Yjq3C_evcUaMdVDNE7aqRv2IAaDYIMZtVCxXzniqgoN5yT5cvJixHUJ20gNBKIXvuj-0FF8xp7Nad4mi677Si3Rfzzc9bVq7P5mCZJ12Z-fSidz_V7AkMdvbPPZiwtbJbbcVvdOsd4BQO8oSEIz4gTLjoRmqE_Hk) for a short overview of this history. Let’s start by thinking of a diploid population forward in time. The term “diploid” refers to organisms that have two copies of each locus in their genome; humans, for example, have 23 pairs of chromosomes. To keep things simple, assume that we randomly select two individuals from the populations to parent each offspring in the next generation. At a particular locus, each individual then randomly passes one of their two copies (alleles) on to their offspring. We repeat this random process until the second generation has the same number of individuals as the starting population, thus maintaining a constant population size over time. Below is an example simulation which shows this played out in a population over five generations that we can use to highlight important features of this model.

<br>

```{r fig.align="center", fig.cap="*This model is called the neutral Wright-Fisher model. Each allele in the starting population is given a different color to help distinguish them throughout the simulation. All alleles are neutral, meaning that they do not affect the likelihood of producing offspring. Note that this model allows for self-fertilization (selfing), a mechanism where a single individual parents both copies of their offspring’s genome. This mechanism is both grounded in biology (more commonly observed in plants than animals) and helps to simplify calculations within the model.*"}
set.seed(10003)

pop_size <- 3
ploidy <- 2
generations <- 5
tracked <- c(sample(((pop_size*ploidy*generations-pop_size*ploidy)+1):(pop_size*ploidy*generations), pop_size*ploidy))

parent_matrix <- generate_parent_matrix(pop_size=pop_size, ploidy=ploidy, generations=generations)
full_lineages <- lapply(1:(pop_size*ploidy*generations), find_lineage, parent_matrix=parent_matrix)
flat_lineages <- do.call(rbind, lapply(1:length(full_lineages), flatten_lineage_list, lineage_list=full_lineages))
flat_lineages$pedigree <- (flat_lineages$ancestor+(ploidy-1)) %/% ploidy
coalescent_events <- locate_coalescent_events(lineages=flat_lineages, tracked=tracked, pop_size=pop_size, ploidy=ploidy, generations=generations)
coalescent_tree <- generate_tree(coalescent_events=coalescent_events, tracked=tracked, generations=generations)
labels <- coalescent_tree[which(coalescent_tree$generation==generations),]

mutation_rate <- 0
polymorphisms <- scatter_mutations(
  parent_matrix=parent_matrix, 
  mutation_rate=mutation_rate, 
  starting_pop=c(0,1,2,3,4,5)
)

ggplot() +
  theme_void() +
  theme(legend.position="none") +
  coord_fixed() +
  geom_circle(
    data=data.frame("point"=1:(pop_size*generations)),
    aes(
      x0=(point-1) %/% (pop_size),
      y0=(pop_size) - ((point-1) %% (pop_size)),
      r=0.25
    ),
    size=1,
    color="black",
    fill="white"
  ) +
  geom_point(
    data=data.frame("point"=1:(pop_size*ploidy*generations)),
    aes(
      (point-1) %/% (pop_size*ploidy),
      (pop_size) - (((point-1) %/% ploidy) %% (pop_size)) - ((point-1) %% ploidy - (ploidy-1)/2)/5
    ),
    size=2
  ) +
  geom_text(
    data=data.frame("point"=((pop_size*ploidy*generations-pop_size*ploidy)+1):(pop_size*ploidy*generations)),
    aes(
      (point-1) %/% (pop_size*ploidy) + 0.5,
      (pop_size) - (((point-1) %/% ploidy) %% (pop_size)) - ((point-1) %% ploidy - (ploidy-1)/2)/5,
      label=point-(pop_size*ploidy*generations-pop_size*ploidy)
    )
  ) +
  geom_line(
    data=flat_lineages,
    aes(
      (ancestor-1) %/% (pop_size*ploidy),
      (pop_size) - (((ancestor-1) %/% ploidy) %% (pop_size)) - ((ancestor-1) %% ploidy - (ploidy-1)/2)/5,
      group=starting_indiv
    ),
    size=1
  ) +
  geom_point(
    data=polymorphisms,
    aes(
      (indiv-1) %/% (pop_size*ploidy),
      (pop_size) - (((indiv-1) %/% ploidy) %% (pop_size)) - ((indiv-1) %% ploidy - (ploidy-1)/2)/5,
      size=2,
      color=allele
    )
  )
```

<br>

This is a population of three diploid individuals tracked over five generations (left to right : past to present). Each individual is represented as a circle with two dots, one for each copy of an allele at a locus. Lines connecting copies show how alleles are passed down over time, so that we can track lineages through the generations (by following a line). Due to the randomness of parentage, it is common that not all of the individuals in the parental generation produce offspring in the next generation. Because of this, overtime, certain lineages are lost and others may grow in frequency; this process is known as genetic drift. In the above simulation, the green allele grows in frequency until it is found in every individual in the population. Without mutation, the overall diversity (number of different alleles) of a population with fixed finite size will decrease over time.

<br>

```{r fig.align="center"}
ggplot() +
  theme_void() +
  theme(legend.position="none") +
  coord_fixed() +
  geom_circle(
    data=data.frame("point"=1:(pop_size*generations)),
    aes(
      x0=(point-1) %/% (pop_size),
      y0=(pop_size) - ((point-1) %% (pop_size)),
      r=0.25
    ),
    size=1,
    color="black",
    fill="white"
  ) +
  geom_point(
    data=data.frame("point"=1:(pop_size*ploidy*generations)),
    aes(
      (point-1) %/% (pop_size*ploidy),
      (pop_size) - (((point-1) %/% ploidy) %% (pop_size)) - ((point-1) %% ploidy - (ploidy-1)/2)/5
    ),
    color="grey",
    size=2
  ) +
  geom_text(
    data=data.frame("point"=((pop_size*ploidy*generations-pop_size*ploidy)+1):(pop_size*ploidy*generations)),
    aes(
      (point-1) %/% (pop_size*ploidy) + 0.5,
      (pop_size) - (((point-1) %/% ploidy) %% (pop_size)) - ((point-1) %% ploidy - (ploidy-1)/2)/5,
      label=point-(pop_size*ploidy*generations-pop_size*ploidy)
    )
  ) +
  geom_line(
    data=flat_lineages,
    aes(
      (ancestor-1) %/% (pop_size*ploidy),
      (pop_size) - (((ancestor-1) %/% ploidy) %% (pop_size)) - ((ancestor-1) %% ploidy - (ploidy-1)/2)/5,
      group=starting_indiv
    ),
    color="grey",
    size=1
  ) +
  geom_line(
    data=flat_lineages[which(flat_lineages$starting_indiv %in% tracked),],
    aes(
      (ancestor-1) %/% (pop_size*ploidy),
      (pop_size) - (((ancestor-1) %/% ploidy) %% (pop_size)) - ((ancestor-1) %% ploidy - (ploidy-1)/2)/5,
      group=starting_indiv
    ),
    color="#4285f4ff",
    size=1
  ) +
  geom_point(
    data=data.frame("point"=c(3,12,21,24)),
    aes(
      (point-1) %/% (pop_size*ploidy),
      (pop_size) - (((point-1) %/% ploidy) %% (pop_size)) - ((point-1) %% ploidy - (ploidy-1)/2)/5
    ),
    color="red",
    size=3
  )
```

```{r fig.align="center", fig.cap="<br><br>*The same simulation as above, now tracing the lineages of the alleles in the present day. All coalescence events are marked with red points. By untangling these lines, we have a genealogical tree that relates all of the alleles.*"}
coalescent_events <- locate_coalescent_events(lineages=flat_lineages, tracked=tracked, pop_size=pop_size, ploidy=ploidy, generations=generations)
coalescent_tree <- generate_tree(coalescent_events=coalescent_events, tracked=tracked, generations=generations)
labels <- coalescent_tree[which(coalescent_tree$generation==generations),]

ggplot() +
  theme_void() +
  theme(legend.position="none") +
  coord_fixed() +
  geom_line(data=coalescent_tree, aes(generation, value, group=lineage), size=1) +
  geom_text(data=labels, aes(generation+0.5, value, label=lineage-(pop_size*ploidy*generations-pop_size*ploidy)))
```

<br>

To understand the coalescent process, we can track present day lineages backwards in time. For example, in the diagram above, you can choose an allele in the present day (on the right) and follow its lineage back through time (from right to left), tracing its blue line all the way back to its ancestors. By tracking multiple lineages at once, we can begin to make sense of the relationships between alleles sampled in the current day. A coalescence event, the namesake of this model, occurs when two copies share a common ancestor. Lineages 4 and 5 from the present generation coalesce one generation in the past; their lines merge into one as we scan to the left. In this simulation, all of the present day lineages share a common ancestor at some point within the past five generations, so we can describe their relationship using a single genealogical tree (shortened to ‘gene tree’ from here on). The root of the tree on the left is the most recent common ancestor of all the lineages, and closely related lineages split more recently on the right side of the tree.

The properties of this gene tree are related to the size of the population; we expect that larger populations will have a gene tree that extends deeper into the past. This is because lineages are less likely to coalesce in a given generation when there are more potential parents due to the larger population size. As this is a random process, it is possible for all lineages to coalesce very quickly, though this is increasingly unlikely at larger population sizes.

Historically when looking at a family tree, we would not have known which copies each parent passed on to their offspring. A pedigree is a graphical representation of parentage, with no genetic information. Pedigrees are very useful in genetics as they allow us to approximate genetic relatedness without having to directly observe which lineages were passed down from parent to offspring, but importantly, they are not a one-to-one replacement for the original simulation. For a given pedigree, there are many possible simulations, and therefore gene trees, that satisfy the defined relationships between individuals. Below is a simple pedigree of two full siblings, who share both parents, with the corresponding four possible inheritance patterns.

<br>

```{r fig.align="center", fig.cap="*Top: A pedigree of two parents (left) and two offspring (right). Using this diagram alone, it is unknown which allele is passed from parent to offspring. Bottom: The four possible genetic inheritance patterns that correspond to the above pedigree.*"}
pop_size <- 2
ploidy <- 2
generations <- 2

parent_matrix <- matrix(NA, nrow=(pop_size*ploidy), ncol=generations)
parent_matrix[1, 2] <- 1
parent_matrix[2, 2] <- 3
parent_matrix[3, 2] <- 2
parent_matrix[4, 2] <- 4

full_lineages <- lapply(1:(pop_size*ploidy*generations), find_lineage, parent_matrix=parent_matrix)
flat_lineages <- do.call(rbind, lapply(1:length(full_lineages), flatten_lineage_list, lineage_list=full_lineages))
flat_lineages$pedigree <- (flat_lineages$ancestor+(ploidy-1)) %/% ploidy
coalescent_events <- locate_coalescent_events(lineages=flat_lineages, tracked=tracked, pop_size=pop_size, ploidy=ploidy, generations=generations)
coalescent_tree <- generate_tree(coalescent_events=coalescent_events, tracked=tracked, generations=generations)
labels <- coalescent_tree[which(coalescent_tree$generation==generations),]

ped <- ggplot() +
  theme_void() +
  theme(legend.position="none") +
  coord_fixed() +
  geom_line(
    data=flat_lineages,
    aes(
      (pedigree-1) %/% (pop_size),
      (pop_size) - ((pedigree-1) %% (pop_size)),
      group=starting_indiv
    ),
    size=1
  ) +
  geom_circle(
    data=data.frame("point"=1:(pop_size*generations)),
    aes(
      x0=(point-1) %/% (pop_size),
      y0=(pop_size) - ((point-1) %% (pop_size)),
      r=0.25
    ),
    size=1,
    color="black",
    fill="white"
  ) +
  geom_text(
    data=data.frame("point"=((pop_size*ploidy*generations-pop_size*ploidy)+1):(pop_size*ploidy*generations)),
    aes(
      (point-1) %/% (pop_size*ploidy),
      (pop_size) - (((point-1) %/% ploidy) %% (pop_size)) - ((point-1) %% ploidy - (ploidy-1)/2)/5,
      label=point
    ),
    color="white"
  )

coal_1 <- ggplot() +
  theme_void() +
  theme(legend.position="none") +
  coord_fixed() +
  geom_circle(
    data=data.frame("point"=1:(pop_size*generations)),
    aes(
      x0=(point-1) %/% (pop_size),
      y0=(pop_size) - ((point-1) %% (pop_size)),
      r=0.25
    ),
    size=1,
    color="black",
    fill="white"
  ) +
  geom_point(
    data=data.frame("point"=1:(pop_size*ploidy*generations)),
    aes(
      (point-1) %/% (pop_size*ploidy),
      (pop_size) - (((point-1) %/% ploidy) %% (pop_size)) - ((point-1) %% ploidy - (ploidy-1)/2)/5
    ),
    color="grey",
    size=2
  ) +
  geom_line(
    data=flat_lineages,
    aes(
      (ancestor-1) %/% (pop_size*ploidy),
      (pop_size) - (((ancestor-1) %/% ploidy) %% (pop_size)) - ((ancestor-1) %% ploidy - (ploidy-1)/2)/5,
      group=starting_indiv
    ),
    color="grey",
    size=1
  ) +
  geom_line(
    data=flat_lineages[which(flat_lineages$starting_indiv %in% c()),],
    aes(
      (ancestor-1) %/% (pop_size*ploidy),
      (pop_size) - (((ancestor-1) %/% ploidy) %% (pop_size)) - ((ancestor-1) %% ploidy - (ploidy-1)/2)/5,
      group=starting_indiv
    ),
    color="#4285f4ff",
    size=1
  )

parent_matrix <- matrix(NA, nrow=(pop_size*ploidy), ncol=generations)
parent_matrix[1, 2] <- 1
parent_matrix[2, 2] <- 3
parent_matrix[3, 2] <- 1
parent_matrix[4, 2] <- 4

full_lineages <- lapply(1:(pop_size*ploidy*generations), find_lineage, parent_matrix=parent_matrix)
flat_lineages <- do.call(rbind, lapply(1:length(full_lineages), flatten_lineage_list, lineage_list=full_lineages))
flat_lineages$pedigree <- (flat_lineages$ancestor+(ploidy-1)) %/% ploidy
coalescent_events <- locate_coalescent_events(lineages=flat_lineages, tracked=tracked, pop_size=pop_size, ploidy=ploidy, generations=generations)
coalescent_tree <- generate_tree(coalescent_events=coalescent_events, tracked=tracked, generations=generations)
labels <- coalescent_tree[which(coalescent_tree$generation==generations),]

coal_2 <- ggplot() +
  theme_void() +
  theme(legend.position="none") +
  coord_fixed() +
  geom_circle(
    data=data.frame("point"=1:(pop_size*generations)),
    aes(
      x0=(point-1) %/% (pop_size),
      y0=(pop_size) - ((point-1) %% (pop_size)),
      r=0.25
    ),
    size=1,
    color="black",
    fill="white"
  ) +
  geom_point(
    data=data.frame("point"=1:(pop_size*ploidy*generations)),
    aes(
      (point-1) %/% (pop_size*ploidy),
      (pop_size) - (((point-1) %/% ploidy) %% (pop_size)) - ((point-1) %% ploidy - (ploidy-1)/2)/5
    ),
    color="grey",
    size=2
  ) +
  geom_line(
    data=flat_lineages,
    aes(
      (ancestor-1) %/% (pop_size*ploidy),
      (pop_size) - (((ancestor-1) %/% ploidy) %% (pop_size)) - ((ancestor-1) %% ploidy - (ploidy-1)/2)/5,
      group=starting_indiv
    ),
    color="grey",
    size=1
  ) +
  geom_line(
    data=flat_lineages[which(flat_lineages$starting_indiv %in% c(5,7)),],
    aes(
      (ancestor-1) %/% (pop_size*ploidy),
      (pop_size) - (((ancestor-1) %/% ploidy) %% (pop_size)) - ((ancestor-1) %% ploidy - (ploidy-1)/2)/5,
      group=starting_indiv
    ),
    color="#4285f4ff",
    size=1
  ) +
  geom_point(
    data=data.frame("point"=c(1,5,7)),
    aes(
      (point-1) %/% (pop_size*ploidy),
      (pop_size) - (((point-1) %/% ploidy) %% (pop_size)) - ((point-1) %% ploidy - (ploidy-1)/2)/5
    ),
    color="#4285f4ff",
    size=2
  )

parent_matrix <- matrix(NA, nrow=(pop_size*ploidy), ncol=generations)
parent_matrix[1, 2] <- 1
parent_matrix[2, 2] <- 4
parent_matrix[3, 2] <- 2
parent_matrix[4, 2] <- 4

full_lineages <- lapply(1:(pop_size*ploidy*generations), find_lineage, parent_matrix=parent_matrix)
flat_lineages <- do.call(rbind, lapply(1:length(full_lineages), flatten_lineage_list, lineage_list=full_lineages))
flat_lineages$pedigree <- (flat_lineages$ancestor+(ploidy-1)) %/% ploidy
coalescent_events <- locate_coalescent_events(lineages=flat_lineages, tracked=tracked, pop_size=pop_size, ploidy=ploidy, generations=generations)
coalescent_tree <- generate_tree(coalescent_events=coalescent_events, tracked=tracked, generations=generations)
labels <- coalescent_tree[which(coalescent_tree$generation==generations),]

coal_3 <- ggplot() +
  theme_void() +
  theme(legend.position="none") +
  coord_fixed() +
  geom_circle(
    data=data.frame("point"=1:(pop_size*generations)),
    aes(
      x0=(point-1) %/% (pop_size),
      y0=(pop_size) - ((point-1) %% (pop_size)),
      r=0.25
    ),
    size=1,
    color="black",
    fill="white"
  ) +
  geom_point(
    data=data.frame("point"=1:(pop_size*ploidy*generations)),
    aes(
      (point-1) %/% (pop_size*ploidy),
      (pop_size) - (((point-1) %/% ploidy) %% (pop_size)) - ((point-1) %% ploidy - (ploidy-1)/2)/5
    ),
    color="grey",
    size=2
  ) +
  geom_line(
    data=flat_lineages,
    aes(
      (ancestor-1) %/% (pop_size*ploidy),
      (pop_size) - (((ancestor-1) %/% ploidy) %% (pop_size)) - ((ancestor-1) %% ploidy - (ploidy-1)/2)/5,
      group=starting_indiv
    ),
    color="grey",
    size=1
  ) +
  geom_line(
    data=flat_lineages[which(flat_lineages$starting_indiv %in% c(6,8)),],
    aes(
      (ancestor-1) %/% (pop_size*ploidy),
      (pop_size) - (((ancestor-1) %/% ploidy) %% (pop_size)) - ((ancestor-1) %% ploidy - (ploidy-1)/2)/5,
      group=starting_indiv
    ),
    color="#4285f4ff",
    size=1
  ) +
  geom_point(
    data=data.frame("point"=c(4,6,8)),
    aes(
      (point-1) %/% (pop_size*ploidy),
      (pop_size) - (((point-1) %/% ploidy) %% (pop_size)) - ((point-1) %% ploidy - (ploidy-1)/2)/5
    ),
    color="#4285f4ff",
    size=2
  )

parent_matrix <- matrix(NA, nrow=(pop_size*ploidy), ncol=generations)
parent_matrix[1, 2] <- 1
parent_matrix[2, 2] <- 4
parent_matrix[3, 2] <- 1
parent_matrix[4, 2] <- 4

full_lineages <- lapply(1:(pop_size*ploidy*generations), find_lineage, parent_matrix=parent_matrix)
flat_lineages <- do.call(rbind, lapply(1:length(full_lineages), flatten_lineage_list, lineage_list=full_lineages))
flat_lineages$pedigree <- (flat_lineages$ancestor+(ploidy-1)) %/% ploidy
coalescent_events <- locate_coalescent_events(lineages=flat_lineages, tracked=tracked, pop_size=pop_size, ploidy=ploidy, generations=generations)
coalescent_tree <- generate_tree(coalescent_events=coalescent_events, tracked=tracked, generations=generations)
labels <- coalescent_tree[which(coalescent_tree$generation==generations),]

coal_4 <- ggplot() +
  theme_void() +
  theme(legend.position="none") +
  coord_fixed() +
  geom_circle(
    data=data.frame("point"=1:(pop_size*generations)),
    aes(
      x0=(point-1) %/% (pop_size),
      y0=(pop_size) - ((point-1) %% (pop_size)),
      r=0.25
    ),
    size=1,
    color="black",
    fill="white"
  ) +
  geom_point(
    data=data.frame("point"=1:(pop_size*ploidy*generations)),
    aes(
      (point-1) %/% (pop_size*ploidy),
      (pop_size) - (((point-1) %/% ploidy) %% (pop_size)) - ((point-1) %% ploidy - (ploidy-1)/2)/5
    ),
    color="grey",
    size=2
  ) +
  geom_line(
    data=flat_lineages,
    aes(
      (ancestor-1) %/% (pop_size*ploidy),
      (pop_size) - (((ancestor-1) %/% ploidy) %% (pop_size)) - ((ancestor-1) %% ploidy - (ploidy-1)/2)/5,
      group=starting_indiv
    ),
    color="grey",
    size=1
  ) +
  geom_line(
    data=flat_lineages[which(flat_lineages$starting_indiv %in% c(5,6,7,8)),],
    aes(
      (ancestor-1) %/% (pop_size*ploidy),
      (pop_size) - (((ancestor-1) %/% ploidy) %% (pop_size)) - ((ancestor-1) %% ploidy - (ploidy-1)/2)/5,
      group=starting_indiv
    ),
    color="#4285f4ff",
    size=1
  ) +
  geom_point(
    data=data.frame("point"=c(1,4,5,6,7,8)),
    aes(
      (point-1) %/% (pop_size*ploidy),
      (pop_size) - (((point-1) %/% ploidy) %% (pop_size)) - ((point-1) %% ploidy - (ploidy-1)/2)/5
    ),
    color="#4285f4ff",
    size=2
  )

lay <- rbind(c(1,1,1,1),
             c(2,3,4,5))
grid.arrange(ped, coal_1, coal_2, coal_3, coal_4, layout_matrix=lay)
```

<br>

Whether you and your (biological) sibling share 0, 1, or 2 alleles at a given locus of your genome is purely chance, reflecting the randomness of Mendelian segregation. Full siblings are expected to share 50% of their autosomal DNA. At any given locus the transmission is random; there is a 25% chance that they share both copies (right diagram), 50% chance that they share one copy (middle two diagrams), and 25% chance that they do not share any copies (left diagram). While in the context of the simulation the randomness of which allele is transmitted by one parent seems trivial it has cascading effects when played out across generations. If we hold the pedigree constant by rerun the randomness of transmission, we generate drastically different gene trees, all of which share the same underlying pedigree. Here is an interactive tool that lets you rerun simulations over a given pedigree to see the changes in possible gene trees:

<br>

```{r echo=FALSE, warning=FALSE, fig.align='center'}
pop_size_app1 <- 5
ploidy_app1 <- 2
generations_app1 <- 10
tracked_app1 <- c(sample(((pop_size_app1*ploidy_app1*generations_app1-pop_size_app1*ploidy_app1)+1):(pop_size_app1*ploidy_app1*generations_app1), pop_size_app1*ploidy_app1))
values <- reactiveValues(
  parent_matrix=generate_parent_matrix(pop_size=pop_size_app1, ploidy=ploidy_app1, generations=generations_app1)
)

wellPanel(style="background-color: white; border: 3px solid #000000; max-height: 100vh;", align='center',
  actionButton("rerun_pedigree", "Try New Pedigree"),
  actionButton("rerun_coalescent", "Generate Another Tree"),
  hr(style = "border-top: 1px solid #000000;"),
  plotOutput("pedigree", height="30vh"),
  fluidRow(style="width: 100%; margin-left: -25px; padding: 20px 5px 0px 5px;",
    splitLayout(cellWidths=c("75%", "1%", "24%"), plotOutput("coalescent", height="30vh"), br(),plotOutput("tree", height="30vh"))
  )
)

show_pedigree_vs_tree(plot_pedigree=TRUE)

observeEvent(input$rerun_pedigree, {
  values$parent_matrix <- generate_parent_matrix(pop_size=pop_size_app1, ploidy=ploidy_app1, generations=generations_app1)
  show_pedigree_vs_tree(plot_pedigree=TRUE)
})

observeEvent(input$rerun_coalescent, {
  show_pedigree_vs_tree()
})
```

<br>

But why should we care about this many-to-one relationship between gene trees and pedigrees? Though he didn’t know the mechanisms, Mendel discovered that the alleles underlying many traits assort independently of one another. Today we understand that alleles on different chromosomes will share a pedigree but are independently transmitted to the next generation from each other. Thus loci on different chromosomes are very unlikely to share a gene tree due to this independent assortment (vanishingly unlikely in large populations). Completely independent assortment is equivalent to clicking the "Another Tree From Pedigree" button in the above simulation. Even loci on the same chromosome can be separated through the process of recombination which increases the independence between loci and so reduces the correlation between their gene trees. Because the loci on a chromosome can sort independently from generation to generation, a population’s genetic history cannot be illustrated with a single gene tree. To truly sketch the history, you would need thousands of inter-related gene trees, tracing a whole collection of regions across the genome.

The regions of my genome that I share in common with my sister represent regions with a coalescence event in our parent’s generation. Regions that we don’t share directly in common from the last generation have a much deeper most recent common ancestor; when looking just at these latter regions, my sister is most likely not my closest genetic relative.

In the previous simulations, offspring inherit an identical version of the sequence of their parent's DNA during reproduction. There’s no genetic variation in these simulations, and so if we were to sequence the individuals in the current day, there’d be no way of reconstructing the gene tree. However, we know that there can be errors in DNA replication (mutations) leading to the development of new alleles, or variants, at a given locus. We can include these mutations in our simulations with two assumptions:

* Mutations are neutral, meaning that they do not affect the relative reproductive success of individuals in the population
* There’s no repeat or back mutation, i.e. mutations that affect precisely the same base pair (this is technically called the infinite sites model)

<br>

```{r fig.align="center"}
set.seed(3)

pop_size <- 3
ploidy <- 2
generations <- 5
tracked <- c(sample(((pop_size*ploidy*generations-pop_size*ploidy)+1):(pop_size*ploidy*generations), pop_size*ploidy))

parent_matrix <- generate_parent_matrix(pop_size=pop_size, ploidy=ploidy, generations=generations)
full_lineages <- lapply(1:(pop_size*ploidy*generations), find_lineage, parent_matrix=parent_matrix)
flat_lineages <- do.call(rbind, lapply(1:length(full_lineages), flatten_lineage_list, lineage_list=full_lineages))
flat_lineages$pedigree <- (flat_lineages$ancestor+(ploidy-1)) %/% ploidy
coalescent_events <- locate_coalescent_events(lineages=flat_lineages, tracked=tracked, pop_size=pop_size, ploidy=ploidy, generations=generations)
coalescent_tree <- generate_tree(coalescent_events=coalescent_events, tracked=tracked, generations=generations)
labels <- coalescent_tree[which(coalescent_tree$generation==generations),]

mutation_rate <- 0.10
polymorphisms <- scatter_mutations(parent_matrix=parent_matrix, mutation_rate=mutation_rate)

ggplot() +
  theme_void() +
  theme(legend.position="none") +
  coord_fixed() +
  geom_circle(
    data=data.frame("point"=1:(pop_size*generations)),
    aes(
      x0=(point-1) %/% (pop_size),
      y0=(pop_size) - ((point-1) %% (pop_size)),
      r=0.25
    ),
    size=1,
    color="black",
    fill="white"
  ) +
  geom_point(
    data=data.frame("point"=1:(pop_size*ploidy*generations)),
    aes(
      (point-1) %/% (pop_size*ploidy),
      (pop_size) - (((point-1) %/% ploidy) %% (pop_size)) - ((point-1) %% ploidy - (ploidy-1)/2)/5
    ),
    size=2
  ) +
  geom_text(
    data=data.frame("point"=((pop_size*ploidy*generations-pop_size*ploidy)+1):(pop_size*ploidy*generations)),
    aes(
      (point-1) %/% (pop_size*ploidy) + 0.5,
      (pop_size) - (((point-1) %/% ploidy) %% (pop_size)) - ((point-1) %% ploidy - (ploidy-1)/2)/5,
      label=point-(pop_size*ploidy*generations-pop_size*ploidy)
    )
  ) +
  geom_line(
    data=flat_lineages,
    aes(
      (ancestor-1) %/% (pop_size*ploidy),
      (pop_size) - (((ancestor-1) %/% ploidy) %% (pop_size)) - ((ancestor-1) %% ploidy - (ploidy-1)/2)/5,
      group=starting_indiv
    ),
    size=1,
    color="grey"
  ) +
  geom_point(
    data=polymorphisms,
    aes(
      (indiv-1) %/% (pop_size*ploidy),
      (pop_size) - (((indiv-1) %/% ploidy) %% (pop_size)) - ((indiv-1) %% ploidy - (ploidy-1)/2)/5,
      size=2,
      color=allele
    )
  )
```

```{r fig.align="center", fig.cap="<br><br>*A different simulation of a diploid Wright-Fisher population of three individuals with each variant given a unique color. All alleles in the starting population are the same. Over time, mutations generate new alleles. The mutation rate is set to 0.1 mutations/transmission. The mutation events are added along the gene tree to generate a pseudo-sequence for the gene; only sites that are polymorphic in the present day population are shown.*"}
associated_tree <- data.frame("lineage"=tracked, "tree"=sapply(tracked, FUN=identify_deepest_ancestor, flat_lineages=flat_lineages))
if (max(as.numeric(polymorphisms$allele))-1 > 0) {
  mutation_events <- data.frame("mutation"=1:(max(as.numeric(polymorphisms$allele))-1), "indiv"=sapply(1:(max(as.numeric(polymorphisms$allele))-1), FUN=locate_mutation_events, polymorphisms=polymorphisms))
  tracked_mutations <- merge(mutation_events, flat_lineages[which(flat_lineages$starting_indiv %in% tracked),], by.x="indiv", by.y="ancestor")
  if (nrow(tracked_mutations) > 0) {
    tracked_mutations$generation <- ((tracked_mutations$indiv-1) %/% (pop_size*ploidy)) + 1
    tracked_mutations <- merge(coalescent_tree, tracked_mutations, by.x=c("lineage", "generation"), by.y=c("starting_indiv", "generation"))
    tracked_mutations <- merge(tracked_mutations, associated_tree, by="lineage")
    tracked_mutation_events <- tracked_mutations %>%
      group_by(mutation, generation, value, tree) %>%
      summarize(
        lineages = paste(lineage, collapse=" "),
        alleles = paste(sample(c("A", "T", "C", "G"), 2), collapse="/"),
        .groups = 'drop'
      ) %>%
      separate(alleles, c("ancestral_allele", "derived_allele")) %>%
      arrange(value)
    tracked_mutation_events$mutation_order <- 1:nrow(tracked_mutation_events)
  } else {
    tracked_mutation_events <- data.frame("mutation"=factor(), "generation"=numeric(), "value"=numeric(), "tree"=factor(), "lineages"=character(), "ancestral_allele"=character(), "derived_allele"=character(), "mutation_order"=numeric())
  }
} else {
  tracked_mutation_events <- data.frame("mutation"=factor(), "generation"=numeric(), "value"=numeric(), "tree"=factor(), "lineages"=character(), "ancestral_allele"=character(), "derived_allele"=character(), "mutation_order"=numeric())
}
if (nrow(tracked_mutation_events) > 0) {
  alleles <- as.data.frame(apply(tracked_mutation_events, FUN=identify_allele_seq, MARGIN=1, tracked=tracked, tree=associated_tree))
  allele_seq <- data.frame("lineage"=rownames(alleles), unite(alleles, "allele_seq", sep=""))
  allele_seq <- merge(allele_seq, coalescent_tree[which(coalescent_tree$generation == generations),], by="lineage")
  ancestral_seq <- paste(tracked_mutation_events$ancestral_allele, collapse="")
  allele_seq <- rbind(data.frame("lineage"="anc","allele_seq"=ancestral_seq,"generation"=0,"value"=nrow(allele_seq)/2+0.5), allele_seq)
} else {
  allele_seq <- data.frame("lineage"=factor(), "allele_seq"=character(), "generation"=numeric(), "value"=numeric())
}



palette <- hue_pal()(4)[3:4]

tree <- ggplot() +
  theme_void() +
  theme(legend.position="none") +
  coord_fixed() +
  geom_line(data=coalescent_tree, aes(generation, value, group=lineage), size=1)

if (max(as.numeric(polymorphisms$allele))-1 > 0 & nrow(tracked_mutation_events > 0)) {
  tree +
    geom_label(data=tracked_mutation_events, aes(generation, value, label=paste(paste(mutation_order, ".", sep=""), paste(ancestral_allele, ">", derived_allele, sep=""))), color=palette) +
    geom_text(data=allele_seq[1,], aes(generation+0.5, value, label=allele_seq)) +
    geom_text(data=allele_seq[2:nrow(allele_seq),], aes(generation+1, value, label=paste0(as.numeric(lineage)-(pop_size*ploidy*generations-pop_size*ploidy), " : ", allele_seq)))
} else {
  tree
}
```

<br>

Each allele is given a unique color; there are four gene variants represented in our simulation. Red is the ancestral allele, whereas green, teal, and purple are derived through mutation events in the second, fourth, and fifth generations, respectively. The green allele is lost from the population within two generations whereas the teal and purple alleles persist to the present. These mutation events can be laid on top of the genealogical tree allowing us to more clearly see which individuals are affected by each mutation. Because the blue allele is not represented in the present day population, it is not included on the tree. To generate a pseudo-sequence of this simulated gene, we assign the ancestral and derived variants for each mutation event to nucleotide bases (A, T, G, or C). Every individual then has a sequence of bases corresponding to the mutation events along its ancestral branches. We expect that lineages that are genealogically closely related will also have more similar sequences.

Given the same pedigree, different loci can have very different patterns of transmission, and so very different coalescent trees and very different mutational patterns in their sequences. This is one of the central insights of population genetics. There’s no one tree; a genome is a whole forest.

To wrap up this post, I've created an interactive tool that gives you control over many parameters which affect the simulation. Test it out for yourself! All of the code for this post can be found altogether as an R Markdown file [here](https://www.github.com/kitchensjn).

<br>

```{r}
run_coalescent_simulation <- function() {
  parent_matrix <- reactive({
    generate_parent_matrix(pop_size=input$pop_size, ploidy=input$ploidy, generations=input$generations, force_coal=input$force_coal)
  })
  
  tracked <- reactive({
    c(sample(((input$pop_size*input$ploidy*input$generations-input$pop_size*input$ploidy)+1):(input$pop_size*input$ploidy*input$generations), input$tracked_lineages)) # random lineages in modern generation to be tracked
  })
  
  polymorphisms <- reactive({
    scatter_mutations(parent_matrix=parent_matrix(), mutation_rate=input$mutation_rate)
  })
  
  flat_lineages <- reactive({
    full_lineages <- lapply(1:(input$pop_size*input$ploidy*input$generations), find_lineage, parent_matrix=parent_matrix())
    do.call(rbind, lapply(1:length(full_lineages), flatten_lineage_list, lineage_list=full_lineages))
  })
  
  
  output$coalescent_app2 <- renderPlot({
    
    if ((input$generations < 20) && (input$pop_size < 7)) {
      coalescent <- ggplot() +
        coord_fixed() +
        geom_circle(
          data=data.frame("point"=1:(input$pop_size*input$generations)),
          aes(
            x0=(point-1) %/% (input$pop_size),
            y0=(input$pop_size) - ((point-1) %% (input$pop_size)),
            r=0.25
          ),
          size=1,
          color="black",
          fill="white"
        )
    } else {
      coalescent <- ggplot()
    }
    
    coalescent <- coalescent +
      theme_void() +
      theme(legend.position="none") +
      geom_point(
        data=data.frame("point"=1:(input$pop_size*input$ploidy*input$generations)),
        aes(
          (point-1) %/% (input$pop_size*input$ploidy),
          (input$pop_size) - (((point-1) %/% input$ploidy) %% (input$pop_size)) - ((point-1) %% input$ploidy - (input$ploidy-1)/2)/5
        ),
        size=1,
        color="grey"
      ) +
      geom_line(
        data=flat_lineages(),
        aes(
          (ancestor-1) %/% (input$pop_size*input$ploidy),
          (input$pop_size) - (((ancestor-1) %/% input$ploidy) %% (input$pop_size)) - ((ancestor-1) %% input$ploidy - (input$ploidy-1)/2)/5,
          group=starting_indiv
        ),
        size=1,
        color="grey"
      )
    
    if (input$lineage_color) {
      coalescent  <- coalescent +
        geom_line(
          data=flat_lineages()[which(flat_lineages()$starting_indiv %in% tracked()),],
          aes(
            (ancestor-1) %/% (input$pop_size*input$ploidy),
            (input$pop_size) - (((ancestor-1) %/% input$ploidy) %% (input$pop_size)) - ((ancestor-1) %% input$ploidy - (input$ploidy-1)/2)/5,
            group=starting_indiv,
            color=starting_indiv
          ),
          size=1
        )
    } else {
      coalescent  <- coalescent +
        geom_line(
          data=flat_lineages()[which(flat_lineages()$starting_indiv %in% tracked()),],
          aes(
            (ancestor-1) %/% (input$pop_size*input$ploidy),
            (input$pop_size) - (((ancestor-1) %/% input$ploidy) %% (input$pop_size)) - ((ancestor-1) %% input$ploidy - (input$ploidy-1)/2)/5,
            group=starting_indiv
          ),
          size=1,
          color="#4285f4ff"
        )
    }
    
    if (input$polymorphism_color) {
      if ((input$generations < 20) && (input$pop_size < 7)) {
        coalescent +
          geom_point(
            data=polymorphisms(),
            aes(
              (indiv-1) %/% (input$pop_size*input$ploidy),
              (input$pop_size) - (((indiv-1) %/% input$ploidy) %% (input$pop_size)) - ((indiv-1) %% input$ploidy - (input$ploidy-1)/2)/5,
              size=2,
              color=allele
            )
          )
      } else {
        coalescent +
          geom_point(
            data=polymorphisms(),
            aes(
              (indiv-1) %/% (input$pop_size*input$ploidy),
              (input$pop_size) - (((indiv-1) %/% input$ploidy) %% (input$pop_size)) - ((indiv-1) %% input$ploidy - (input$ploidy-1)/2)/5,
              color=allele
            )
          )
      }
    } else {
      coalescent
    }
  })
    
  coalescent_tree <- reactive({
    coalescent_events <- locate_coalescent_events(lineages=flat_lineages(), tracked=tracked(), pop_size=input$pop_size, ploidy=input$ploidy, generations=input$generations)
    generate_tree(coalescent_events=coalescent_events, tracked=tracked(), generations=input$generations)
  })
  
  associated_tree <- reactive({
    data.frame("lineage"=tracked(), "tree"=sapply(tracked(), FUN=identify_deepest_ancestor, flat_lineages=flat_lineages()))
  })
  
  tracked_mutation_events <- reactive({
    if (max(as.numeric(polymorphisms()$allele))-1 > 0) {
      mutation_events <- data.frame("mutation"=1:(max(as.numeric(polymorphisms()$allele))-1), "indiv"=sapply(1:(max(as.numeric(polymorphisms()$allele))-1), FUN=locate_mutation_events, polymorphisms=polymorphisms()))
      tracked_mutations <- merge(mutation_events, flat_lineages()[which(flat_lineages()$starting_indiv %in% tracked()),], by.x="indiv", by.y="ancestor")
      if (nrow(tracked_mutations) > 0) {
        tracked_mutations$generation <- ((tracked_mutations$indiv-1) %/% (input$pop_size*input$ploidy)) + 1
        tracked_mutations <- merge(coalescent_tree(), tracked_mutations, by.x=c("lineage", "generation"), by.y=c("starting_indiv", "generation"))
        tracked_mutations <- merge(tracked_mutations, associated_tree(), by="lineage")
        tracked_mutation_events <- tracked_mutations %>%
          group_by(mutation, generation, value, tree) %>%
          summarize(
            lineages = paste(lineage, collapse=" "),
            alleles = paste(sample(c("A", "T", "C", "G"), 2), collapse="/"),
            .groups = 'drop'
          ) %>%
          separate(alleles, c("ancestral_allele", "derived_allele")) %>%
          arrange(value)
        tracked_mutation_events$mutation_order <- 1:nrow(tracked_mutation_events)
        tracked_mutation_events
      } else {
        data.frame("mutation"=factor(), "generation"=numeric(), "value"=numeric(), "tree"=factor(), "lineages"=character(), "ancestral_allele"=character(), "derived_allele"=character(), "mutation_order"=numeric())
      }
    } else {
      data.frame("mutation"=factor(), "generation"=numeric(), "value"=numeric(), "tree"=factor(), "lineages"=character(), "ancestral_allele"=character(), "derived_allele"=character(), "mutation_order"=numeric())
    }
  })
  
  allele_seq <- reactive({
    if (nrow(tracked_mutation_events()) > 0) {
      alleles <- as.data.frame(apply(tracked_mutation_events(), FUN=identify_allele_seq, MARGIN=1, tracked=tracked(), tree=associated_tree()))
      allele_seq <- data.frame("lineage"=rownames(alleles), unite(alleles, "allele_seq", sep=""))
      merge(allele_seq, coalescent_tree()[which(coalescent_tree()$generation == input$generations),], by="lineage")
    } else {
      data.frame("lineage"=factor(), "allele_seq"=character(), "generation"=numeric(), "value"=numeric())
    }
  })
  
  ancestral_seq <- reactive({
    if (nrow(allele_seq()) > 0) {
      data.frame("lineage"="anc", "allele_seq"=paste(tracked_mutation_events()$ancestral_allele, collapse=""), "generation"=0, "value"=nrow(allele_seq())/2+0.5)
    } else {
      data.frame("lineage"=factor(), "allele_seq"=character(), "generation"=numeric(), "value"=numeric())
    }
  })
  
  output$tree_app2 <- renderPlot({
    if (length(tracked())>1) {
      tree <- ggplot() +
      theme_void() +
      theme(legend.position="none") +
      geom_line(data=coalescent_tree(), aes(generation, value, group=lineage), size=1)
    
      if (max(as.numeric(polymorphisms()$allele))-1 > 0 & nrow(tracked_mutation_events()) > 0) {
        tree +
          geom_label(data=tracked_mutation_events(), aes(generation, value, label=paste(paste(mutation_order, ".", sep=""), paste(ancestral_allele, ">", derived_allele, sep="")))) +
          geom_text(data=ancestral_seq(), aes(generation-(length(allele_seq[1])*input$generations/20), value, label=allele_seq)) +
          geom_text(data=allele_seq(), aes(generation+(length(allele_seq[1])*input$generations/20), value, label=paste0(as.numeric(lineage)-(input$pop_size*input$ploidy*input$generations-input$pop_size*input$ploidy), " : ", allele_seq)))
      } else {
        tree
      } 
    } else {
      ggplot()
    }
  }) 
}

run_coalescent_simulation()

observeEvent(input$rerun_coalescent_simulation, {
  run_coalescent_simulation()
})

wellPanel(style="background-color: white; border: 3px solid #000000; height: auto;", align='center',
  fluidRow(
    column(6, sliderInput(inputId="generations", label="Generations", min=10, max=50, value=10)),
    column(6, sliderInput(inputId="pop_size", label="Population Size", min=5, max=10, value=5))
  ),
  fluidRow(
    column(2, sliderInput(inputId="ploidy", label="Ploidy", min=1, max=2, value=1, step=1)),
    column(5, sliderInput(inputId="mutation_rate", label="Mutation Rate", min=0, max=0.25, value=0.05, step=0.05)),
    column(5, sliderInput(inputId="tracked_lineages", label="Tracked Lineages", min=1, max=5, value=5)),
  ),
  fluidRow(
    column(6, h3("Colors"), fluidRow(
      column(6, checkboxInput("lineage_color", "Lineages", value = FALSE)),
      column(6, checkboxInput("polymorphism_color", "Mutations", value = FALSE)),
    )),
    column(6, checkboxInput("force_coal", "Force Coalescence In Ancestral Generation", value = FALSE), actionButton("rerun_coalescent_simulation", "Rerun With Same Parameters"))
  ),
  hr(style = "border-top: 1px solid #000000;"),
  plotOutput("coalescent_app2", height="50vh"),
  fluidRow(style="padding-top: 20px;", plotOutput("tree_app2", height="50vh"))
)
```

<br><br><br><br>