phase-subclonal-mutations.Rmd

---
title: "3. Phase Subclonal Mutations"
author:
    - "Will Hannon"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output: 
  html_document:
  highlight: tango
  number_sections: no
  theme: default
  toc: yes
  toc_depth: 3
  toc_float:
    collapsed: no
    smooth_scroll: yes
---

The goal of this notebook is to phase more clusters of variants using the same approach we used to define the 'genotypes' `Genome-1` and `Genome-2` in the previous notebook. **The output of this notebook is the list of variants annotated by whether they are part of any haplotypes/clusters**. 

```{r Setup, include=F}
require("knitr")
knitr::opts_chunk$set(echo = FALSE)
```

```{r Required Packages, message=F, warning=F, echo=F}

## ==== Install Required Packages ==== ##

## List of all packages needed
packages = c("tidyverse", "cluster", "ggrepel")

## Packages loading
invisible(lapply(c(packages), library, character.only = TRUE))

```

The inputs for this analysis are (1) the filtered and combined variants annotated by major genotypes from `02-determine-main-genotypes.Rmd` and (2) annotations of the positions of the main MeV open reading frames. 

```{r Inputs, echo=T}

## ==== File paths input data ==== ##

if (exists("snakemake")) {
  
  # Variant data annotated with major genotypes
  labeled.data = snakemake@params[['incsv']]
  
  # Annotations for the standard Measles ORFs  
  annotations.filepath = snakemake@params[['annotations']]
  
  # Path to the ancestral mutations
  ancestral.data = snakemake@params[['ancestral']]

} else {
  
  # Variant data annotated with major genotypes
  labeled.data = "../../results/variants/genotyped_variants.csv"
  
  # Annotations for the standard Measles ORFs  
  annotations.filepath = "../../config/annotations.csv"
  
  # Path to the ancestral mutations
  ancestral.data = "../../config/ref/annotated_SSPE_consensus_snps.csv"

}


```

```{r Outputs, echo=T}

## ==== File paths output data ==== ##

# Path to save the results
if (exists("snakemake")) {
  
  # Variants annotated by new haplotype clusters
  output.path = snakemake@params[['outcsv']]
  
  # Path to save figures
  figure.dir = paste0(snakemake@params[['figures']], "/")

} else {
  
  # Variants annotated by new haplotype clusters
  output.path = "../../results/variants/clustered_variants.csv"
  
  # Path to save figures
  figure.dir = "../../results/figures/"
  
}

# Make the figure directory if it doesn't exist
if (!file.exists(figure.dir)) {
  dir.create(figure.dir, recursive = TRUE)
  print("Directory created.")
} else {
  print("Directory already exists.")
}

```

## Identify Subclonal Haplotypes

Since the goal of this notebook is use the same technique that I used to identify `Genome-1`, `Genome-2`, and `Genome-1-1`, but for mutations that *aren't necessarily found in every single tissue*, I need to account for missing mutations. To do this, I'll spread the data frame of frequencies and assign an allele frequency of 0 if a variant is missing from a tissue. This will allow us to cluster all of the mutations, even if they weren't observed in a tissue sample.

```{r Process for Haplotyping, warning=F, message=F, echo=F}

# Read in the variants with G1, G2, and G1-1 labeled
labeled.df = read_csv(labeled.data, show_col_types = FALSE)

haplotypes.label = labeled.df %>% 
  select(SNP, Haplotype, Background) %>% 
  distinct()

# Expand the mutations to have a frequency for every tissue.
expanded.df = labeled.df %>% 
  select(SNP, Tissue, AF) %>% 
  pivot_wider(names_from = "Tissue", values_from = "AF", values_fill = 0) %>% 
  pivot_longer(cols = !SNP, names_to = "Tissue", values_to = "AF") %>% 
  left_join(., select(labeled.df, c("SNP", "Tissue", "DP")), by = c("SNP", "Tissue")) %>% 
  mutate(DP = if_else(is.na(DP), 0, DP)) %>% 
  left_join(., haplotypes.label, by = "SNP") 

# Get the mean frequency of the major genomes 
tissue.mean = labeled.df %>% 
  filter(Haplotype %in% c("genome-1", "genome-2")) %>% 
  group_by(Tissue, Haplotype) %>% 
  summarize(AF.mean = mean(AF, na.rm = TRUE), 
            SD = sd(AF, na.rm = TRUE),
            N = n()) %>% 
  mutate(SE = SD / sqrt(N),
         Lower.CI = qt(1 - (0.05 / 2), N - 1) * SE,
         Upper.CI = qt(1 - (0.05 / 2), N - 1) * SE) %>% 
  rename("AF" = AF.mean)

# Tissue order - roughly by the location in the brain 
tissue_order = c(
  "SSPE 1", 
  "SSPE 2",
  "Frontal Cortex 2", 
  "Frontal Cortex 1", 
  "Frontal Cortex 3", 
  "Parietal Lobe",
  "Temporal Lobe", 
  "Occipital Lobe",
  "Hippocampus",
  "Internal Capsule", 
  "Midbrain", 
  "UBS",
  "Brain Stem",
  "Cerebellum",
  "Cerebellum Nucleus"
)

```

I'll write the approach I used in `02-determine-main-genotypes.Rmd` into a reusable function. I can provide a list of SNPs, the full data frame, and the number of clusters I'm targeting.

```{r Function to Cluster SNPs by Frequency, echo=T}

cluster.snps <- function(list.of.snps, snp.df, n.clusters) {
  
  # SNP as column and Allele frequency as row
  frequency.by.snp = snp.df %>% 
    filter(SNP %in% list.of.snps) %>% 
    select(AF, SNP, Tissue) %>% 
    pivot_wider(names_from = SNP, values_from = AF, values_fill = 0) %>% 
    select(!Tissue) 
  
  # Calculate R between every pair of columns while handling NAs 
  snp.correlation = cor(frequency.by.snp, use="pairwise.complete.obs")
  
  # Convert to distance (positive corr is close to 0)
  snp.dist = as.dist(1 - snp.correlation)
  
  # Create k-medoids clustering with n clusters
  snp.kmedoids = pam(snp.dist, n.clusters) 
  kclusters = snp.kmedoids$cluster
  
  # Convert to a data frame
  kclusters.df = data.frame(SNP = names(kclusters), cluster = kclusters)

  # Assign the clusters to the original data frame
  kmedoids.SNPs = snp.df %>% 
    filter(SNP %in% list.of.snps) %>% 
    left_join(., kclusters.df, by = "SNP") %>% 
    mutate(cluster = if_else(is.na(cluster), "no cluster", as.character(cluster)))
  
  # Sort the cluster names by mean frequency 
  cluster.order = kmedoids.SNPs %>% 
    group_by(cluster) %>% 
    summarize(AF = mean(AF)) %>% 
    arrange(-AF) %>% 
    mutate(order = 1:n.clusters) %>% 
    select(cluster, order)
  kmedoids.SNPs = kmedoids.SNPs %>% 
    left_join(., cluster.order, by = "cluster") %>% 
    select(!cluster) %>% 
    rename(cluster = order) %>% 
    mutate(cluster = as.character(cluster))
  
  # Add the number of SNPs per cluster
  n.clusters.per.snp = kmedoids.SNPs %>% 
    select(SNP, cluster) %>%
    distinct() %>% 
    group_by(cluster) %>% 
    count() %>% 
    mutate(cluster_size = paste0(cluster, " (", n, " snps in cluster)"))
  
  return(left_join(kmedoids.SNPs, n.clusters.per.snp, by = "cluster"))
  
}

```

Although this approach works reasonably well, it's *very susceptible to noise*. This means that you can't just provide *all* the variants at once and expect the algorithm to partition them into correct clusters. Instead, the best approach I found was to break variants down into bins based on their frequency and cluster these separately. The higher average frequency, the less noise, and the better the algorithm is at partitioning the variants into meaningful clusters. 

### >= 25% Allele Frequency Mutations

First, I'll start with a 'high-frequency' bin. This includes all variants that reach greater than or equal to 25% in at least one tissue. 

```{r Twentyfive percent or more clusters, warning=F, message=F, fig.align='center', fig.width=19, fig.height=10, echo=F}

# SNPs that are above 25% at some point in some tissue
twentyfive.percent.or.more.snps = expanded.df %>% 
  filter(Haplotype == "subclonal") %>% # Not including the mutations we already haplotyped
  filter(AF >= 0.25) %>%
  pull(SNP) %>% 
  unique(.)

# Cluster the SNPs
twentyfive.percent.or.more.clusters.df = cluster.snps(list.of.snps = twentyfive.percent.or.more.snps, 
                                                      snp.df = expanded.df, 
                                                      n.clusters = 20)


# Mean of each tissue ~ for plotting
twentyfive.percent.or.more.clusters.mean = twentyfive.percent.or.more.clusters.df %>% 
  group_by(Tissue, cluster, cluster_size, n) %>% 
  summarize(AF = mean(AF)) 


# Plot the clusters
twentyfive.percent.or.more.clusters.df %>% 
    ggplot(aes(x = factor(Tissue, levels = tissue_order), y = AF)) +
      geom_line(aes(group = SNP)) +
      geom_line(data = twentyfive.percent.or.more.clusters.mean, aes(x = Tissue, y = AF, group = 1, col = (cluster)), size = 1) +
      geom_ribbon(data = tissue.mean, aes(x=Tissue, ymin=AF-SD, ymax=AF+SD, group = Haplotype, fill = Haplotype), alpha=0.2, colour = NA) +
      facet_wrap(~cluster_size) + 
      scale_fill_manual(values=c("#424ef5", "#cf1919")) + 
      xlab("Tissue") + 
      ylab("Allele Frequency") +
      labs(fill="Haplotype", col = "Cluster") +
      theme_bw(20) +
      theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
      theme(strip.text.x = element_text(size = 12))

```

Some clusters have only a single mutation in them. Clearly, some of these belong to `Genome-1`, `Genome-2`, and `Genome-1-1`, but are missing from a single tissue, perhaps due to low coverage in that particular tissue. 

```{r Resolving single clusters, warning=F, message=F, fig.align='center', fig.width=19, fig.height=10, echo=F}

high.frequency.not.clustered.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(n == 1) %>% # clusters that have only a single SNP in them
  pull(SNP)
  
expanded.df %>% 
  filter(SNP %in% high.frequency.not.clustered.snps) %>% 
  left_join(., distinct(select(labeled.df, SNP, Gene_Name, AA_Change))) %>% 
  mutate(mutation.type = paste(Gene_Name, SNP, sep = "-")) %>% 
    ggplot(aes(x = factor(Tissue, levels = tissue_order), y = AF)) +
      geom_line(aes(group = SNP)) +
      geom_ribbon(data = tissue.mean, aes(x=Tissue, ymin=AF-SD, ymax=AF+SD, group = Haplotype, fill = Haplotype), alpha=0.2, colour = NA) +
      facet_wrap(~mutation.type) + 
      scale_fill_manual(values=c("#424ef5", "#cf1919")) + 
      xlab("Tissue") + 
      ylab("Allele Frequency") +
      labs(fill="Haplotype", col = "Cluster") +
      theme_bw(20) +
      theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
      theme(strip.text.x = element_text(size = 12))

```

Cluster 5 (L-A15084G) is a mutation in `Genome-1-1` that was too low frequency to be called in the `Frontal Cortex 2` sample. Cluster 6 (M-C4502T) is mutation that's part of `Genome-2` but was too low coverage in the `Cerebellum`. Finally, Cluster 9 (M-C4573T) is a mutation in `Genome-1` that was also too low coverage in the `Cerebellum`. Generally, the `Cerebellum` has the lowest coverage of any tissue, so this isn't surprising. All-in-all, this adds three mutations to already existing haplotypes. I already addressed F-C7036T, N-G810A, and H-T7293C in the `identify-backgrounds.Rmd` notebook. Finally, M-C4532T might also belong as part of the reference?

I'll annotate these variants in the main data frame.

```{r Annotate single clusters, fig.align='center', fig.width=19, fig.height=10, echo=T}

# Clusters that are clearly genome-1 
genome.1.missing.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(4)) %>% 
  pull(SNP) %>% 
  unique(.)

# Clusters that are clearly genome-1-1
genome.1.1.missing.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(5)) %>% 
  pull(SNP) %>% 
  unique(.)

# Clusters that are clearly genome-2 
genome.2.missing.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(6)) %>% 
  pull(SNP) %>% 
  unique(.)

# Clusters that break the rules 
maybe.in.both.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(3, 7, 2, 1)) %>% 
  pull(SNP) %>% 
  unique(.)


# Annotate the labeled data frame with what's possible here
updated.haplotype.df = labeled.df %>%
  mutate(Background = case_when(SNP %in% genome.1.missing.snps ~ "genome-1",
                                SNP %in% genome.1.1.missing.snps ~ "genome-1",
                                SNP %in% genome.2.missing.snps ~ "genome-2",
                                SNP %in% maybe.in.both.snps ~ "both",
                                TRUE ~ Background)) %>% 
  mutate(Haplotype = case_when(SNP %in% genome.1.missing.snps ~ "genome-1",
                               SNP %in% genome.1.1.missing.snps ~ "genome-1-1",
                               SNP %in% genome.2.missing.snps ~ "genome-2",
                               SNP %in% maybe.in.both.snps ~ "both",
                               TRUE ~ Haplotype))


```

Now, lets look at the remainder of the clusters of SNPs. These are SNPs that actually clustered with other variants. 

```{r Multi SNP Clusters, fig.align='center', fig.width=19, fig.height=10, echo=F}

twentyfive.percent.or.more.clusters.df %>% 
  filter(n > 1) %>% 
    ggplot(aes(x = factor(Tissue, levels = tissue_order), y = AF)) +
      geom_line(data = filter(twentyfive.percent.or.more.clusters.df, n > 1),
                aes(x = factor(Tissue, levels = tissue_order), y = AF, group = SNP, col = (cluster)), size = 1) +
      geom_ribbon(data = tissue.mean, aes(x=factor(Tissue, levels = tissue_order), ymin=AF-SD, ymax=AF+SD, group = Haplotype, fill = Haplotype), alpha=0.2, colour = NA) +
      facet_wrap(~cluster_size) + 
      scale_fill_manual(values=c("#424ef5", "#cf1919")) + 
      xlab("Tissue") + 
      ylab("Allele Frequency") +
      labs(fill="Haplotype", col = "Cluster") +
      theme_bw(20) +
      theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
      theme(strip.text.x = element_text(size = 12))

```

All of these look like reasonable clusters, so I'll annotate the variant data frame with these new clusters (subclonal haplotypes). I'll also try to guess if these are on the background of `Genome-1` or `Genome-2` if it's possible.

```{r Annotate Multi SNP Clusters, fig.align='center', fig.width=19, fig.height=10, echo=T}

cluster.1.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(8)) %>% 
  pull(SNP) %>% 
  unique(.)

# The missing SNPs in Cerebellum are due to coverage
cluster.2.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(11, 14)) %>% 
  pull(SNP) %>% 
  unique(.)

cluster.3.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(10)) %>% 
  pull(SNP) %>% 
  unique(.)

cluster.4.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(13)) %>% 
  pull(SNP) %>% 
  unique(.)

cluster.5.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(15)) %>% 
  pull(SNP) %>% 
  unique(.)

cluster.6.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(16)) %>% 
  pull(SNP) %>% 
  unique(.)

cluster.7.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(17)) %>% 
  pull(SNP) %>% 
  unique(.)

cluster.8.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(18)) %>% 
  pull(SNP) %>% 
  unique(.)

cluster.9.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(19)) %>% 
  pull(SNP) %>% 
  unique(.)

cluster.10.snps = twentyfive.percent.or.more.clusters.df %>% 
  filter(cluster %in% c(20)) %>% 
  pull(SNP) %>% 
  unique(.)
  
# Annotate the labeled data frame with what's possible here
updated.haplotype.df = updated.haplotype.df %>%
  mutate(Haplotype = case_when(SNP %in% cluster.1.snps ~ "cluster 1",
                               SNP %in% cluster.2.snps ~"cluster 2",
                               SNP %in% cluster.3.snps ~"cluster 3",
                               SNP %in% cluster.4.snps ~ "cluster 4",
                               SNP %in% cluster.5.snps ~ "cluster 5",
                               SNP %in% cluster.6.snps ~ "cluster 6",
                               SNP %in% cluster.7.snps ~ "cluster 7",
                               SNP %in% cluster.8.snps ~ "cluster 8",
                               SNP %in% cluster.9.snps ~ "cluster 9",
                               SNP %in% cluster.10.snps ~ "cluster 10",
                               TRUE ~ Haplotype)) 

```

### Haplotypes >= 5% and < 25% 

Now, I'll move onto the next bin. These are mutations that reach at least 5% in one tissue but weren't part of the previous analysis. Most of these will be challenging to haplotype due to noise. I'll only try to specify the clusters that seem the most clear cut. 

```{r Twentyfive percent or fewer clusters, warning=F, message=F, fig.align='center', fig.width=19, fig.height=15, echo=F}

twentyfive.percent.or.fewer.snps = expanded.df %>% 
  filter(!SNP %in% twentyfive.percent.or.more.snps) %>% 
  filter(Haplotype == "subclonal") %>% 
  filter(AF < 0.25 & AF >= 0.05) %>%
  pull(SNP) %>% 
  unique(.)

# Cluster the SNPs
twentyfive.percent.or.fewer.clusters.df = cluster.snps(list.of.snps = twentyfive.percent.or.fewer.snps, 
                                      snp.df = expanded.df, 
                                      n.clusters = 40)

# Mean of each tissue
twentyfive.percent.or.fewer.clusters.mean = twentyfive.percent.or.fewer.clusters.df %>% 
  group_by(Tissue, cluster, cluster_size, n) %>% 
  summarize(AF = mean(AF))

# Plot the clusters
twentyfive.percent.or.fewer.clusters.df %>% 
  filter(n > 1) %>% 
    ggplot(aes(x = factor(Tissue, levels = tissue_order), y = AF)) +
      geom_line(aes(group = SNP)) +
      geom_line(data = filter(twentyfive.percent.or.fewer.clusters.mean, n > 1), aes(x = Tissue, y = AF, group = 1, col = (cluster)), size = 1) +
      geom_ribbon(data = tissue.mean, aes(x=Tissue, ymin=AF-SD, ymax=AF+SD, group = Haplotype, fill = Haplotype), alpha=0.2, colour = NA) +
      facet_wrap(~cluster_size) + 
      scale_fill_manual(values=c("#424ef5", "#cf1919")) + 
      xlab("Tissue") + 
      ylab("Allele Frequency") +
      labs(fill="Haplotype", col = "Cluster") +
      theme_bw(20) +
      theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
      theme(strip.text.x = element_text(size = 12))

```

Most of these aren't very obviously clusters that represent true haplotypes. They're more likely just an artifact of noise. I took the 8 that looked the most promising. It's probably possible to do this more systematically using the internal variation in a cluster, but the results will be the same and there aren't enough clusters to justify the need for this approach. 

```{r Promising low-frequency clusters, fig.align='center', fig.width=19, fig.height=10, echo=F}

promising.subclonal.clusters = c(15, 4, 14, 12)

twentyfive.percent.or.fewer.clusters.df %>% 
  filter(cluster %in% promising.subclonal.clusters) %>% 
    ggplot(aes(x = factor(Tissue, levels = tissue_order), y = AF)) +
      geom_line(aes(group = SNP, col = (cluster))) +
      geom_ribbon(data = tissue.mean, aes(x=factor(Tissue, levels = tissue_order), ymin=AF-SD, ymax=AF+SD, group = Haplotype, fill = Haplotype), alpha=0.2, colour = NA) +
      facet_wrap(~cluster_size) + 
      scale_fill_manual(values=c("#424ef5", "#cf1919")) + 
      xlab("Tissue") + 
      ylab("Allele Frequency") +
      labs(fill="Haplotype", col = "Cluster") +
      theme_bw(20) +
      theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
      theme(strip.text.x = element_text(size = 12))

```

```{r Annotate Remaining Clusters, echo = T}

cluster.11.snps = twentyfive.percent.or.fewer.clusters.df %>% 
  filter(cluster %in% c(15)) %>% 
  pull(SNP) %>% 
  unique(.)

cluster.12.snps = twentyfive.percent.or.fewer.clusters.df %>% 
  filter(cluster %in% c(4)) %>% 
  pull(SNP) %>% 
  unique(.)

cluster.13.snps = twentyfive.percent.or.fewer.clusters.df %>% 
  filter(cluster %in% c(14)) %>% 
  pull(SNP) %>% 
  unique(.)

cluster.14.snps = twentyfive.percent.or.fewer.clusters.df %>% 
  filter(cluster %in% c(12)) %>% 
  pull(SNP) %>% 
  unique(.)

# Annotate the labeled data frame with what's possible here
updated.haplotype.df = updated.haplotype.df %>%
  mutate(Haplotype = case_when(SNP %in% cluster.11.snps ~ "cluster 11",
                               SNP %in% cluster.12.snps ~ "cluster 12",
                               SNP %in% cluster.13.snps ~ "cluster 13",
                               SNP %in% cluster.14.snps ~ "cluster 14",
                               TRUE ~ Haplotype)) 

```

### What's left to cluster? 

We were able to cluster a reasonable number of the high-frequency mutations (>= 5% frequency). How many mutations haven't been haplotyped yet? 

```{r Whats remaining, message=F, warning=F, fig.align='center', fig.width= 5, fig.height=5, echo=F}

# Still labeled as 'subclonal' 
mutations.with.no.info = updated.haplotype.df %>% 
  filter(Background == 'subclonal') %>% 
  pull(SNP) %>% 
  unique()
print(paste("There are about", length(mutations.with.no.info), "that are still totally un-haplotyped."))


# Never above 5% - nearly impossible to haplotype 
mutations.never.above.5perc = updated.haplotype.df %>% 
  filter(Background == 'subclonal') %>% 
  group_by(SNP) %>% 
  summarise(Max = max(AF)) %>% 
  filter(Max < 0.05) %>% 
  pull(SNP) %>% 
  unique()

# Mutations that have some info
mutations.with.info = updated.haplotype.df %>% 
  filter(Background != 'subclonal') %>% 
  pull(SNP) %>% 
  unique()
print(paste("There are", length(mutations.with.info), "that have been clustered."))


# Mutations that are left, but are above 5% some of the time. 
mutations.left.to.haplotype = mutations.with.no.info[which(!mutations.with.no.info %in% mutations.never.above.5perc)]
print(paste("There are still about", length(mutations.left.to.haplotype), "that can reasonably be haplotyped (Above 5% in more at least one sample)."))

# What's the distribution of the frequency of these mutations?
updated.haplotype.df %>% 
  filter(SNP %in% mutations.left.to.haplotype) %>% 
  ggplot(aes(x = AF)) + 
    geom_histogram(bins = 30) + 
    theme_bw()


```

### All Current Haplotype Clusters 

Here are the frequency of all current haplotypes in the 15 tissue samples. 

```{r All Current Clusters, message=F, warning=F, fig.align='center', fig.width=19, fig.height=10, echo=F}

haplotypes.label = updated.haplotype.df %>% 
  select(SNP, Haplotype, Background) %>% 
  distinct()

haplotype.order = c("cluster 1",
                    "cluster 2",
                    "cluster 3",
                    "cluster 4",
                    "cluster 5",
                    "cluster 6",
                    "cluster 7",
                    "cluster 8",
                    "cluster 9",
                    "cluster 10",
                    "cluster 11",
                    "cluster 12",
                    "cluster 13",
                    "cluster 14", 
                    "genome-1-1")
  

# Expand the mutations to have a frequency for every tissue.
expanded.df = updated.haplotype.df %>% 
  select(SNP, Tissue, AF) %>% 
  pivot_wider(names_from = "Tissue", values_from = "AF", values_fill = 0) %>% 
  pivot_longer(cols = !SNP, names_to = "Tissue", values_to = "AF") %>% 
  left_join(., select(labeled.df, c("SNP", "Tissue", "DP")), by = c("SNP", "Tissue")) %>% 
  mutate(DP = if_else(is.na(DP), 0, DP)) %>% 
  left_join(., haplotypes.label, by = "SNP") 

# Get the mean frequency of the major genomes 
tissue.mean = updated.haplotype.df %>% 
  filter(Haplotype %in% c("genome-1", "genome-2")) %>% 
  group_by(Tissue, Haplotype) %>% 
  summarize(AF.mean = mean(AF, na.rm = TRUE), 
            SD = sd(AF, na.rm = TRUE),
            N = n()) %>% 
  mutate(SE = SD / sqrt(N),
         Lower.CI = qt(1 - (0.05 / 2), N - 1) * SE,
         Upper.CI = qt(1 - (0.05 / 2), N - 1) * SE) %>% 
  rename("AF" = AF.mean, "Genotype" = Haplotype)

haplotype.mean = expanded.df %>% 
  group_by(Tissue, Haplotype) %>% 
  summarize(AF = mean(AF))

expanded.df %>% 
    filter(!Haplotype %in% c("fixed", "both", "subclonal", "genome-1", "genome-2")) %>% 
    ggplot(aes(x = factor(Tissue, levels = tissue_order), y = AF)) +
      geom_line(aes(group = SNP)) +
      geom_line(data = filter(haplotype.mean, !Haplotype %in% c("fixed", "both", "subclonal", "genome-1", "genome-2")), 
                aes(x = factor(Tissue, levels = tissue_order), y = AF, group = 1, col = (Haplotype)), size = 1) +
      geom_ribbon(data = tissue.mean, aes(x=factor(Tissue, levels = tissue_order), ymin=AF-SD, ymax=AF+SD, group = Genotype, fill = Genotype), alpha=0.2, colour = NA) +
      facet_wrap(~factor(Haplotype, levels = haplotype.order)) + 
      scale_fill_manual(values=c("#424ef5", "#cf1919")) + 
      xlab("Tissue") + 
      ylab("Allele Frequency") +
      labs(fill="Genotype", col = "Cluster") +
      theme_bw(20) +
      theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
      theme(strip.text.x = element_text(size = 12))

```

Looking closely at the haplotypes, there is one amendment to be made. Cluster 10 and 14 are likely the same cluster. I think these weren't clustered together because the mean frequency of these mutations is split between the frequency bins that I used to partition mutations.

```{r, C14 C10, fig.align='center', fig.width=19, fig.height=10}

expanded.df %>% 
  filter(Haplotype %in% c("cluster 14", "cluster 10")) %>% 
    ggplot(aes(x = factor(Tissue, levels = tissue_order), y = AF)) +
      geom_line(aes(group = SNP)) +
      geom_ribbon(data = tissue.mean, aes(x=factor(Tissue, levels = tissue_order), ymin=AF-SD, ymax=AF+SD, group = Genotype, fill = Genotype), alpha=0.2, colour = NA) +
      facet_wrap(~factor(Haplotype, levels = haplotype.order)) + 
      scale_fill_manual(values=c("#424ef5", "#cf1919")) + 
      xlab("Tissue") + 
      ylab("Allele Frequency") +
      labs(fill="Genotype", col = "Cluster") +
      theme_bw(20) +
      theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
      theme(strip.text.x = element_text(size = 12))

```

For now, I'll combine these into a single cluster. I'll explore this further with bridging reads later in the analysis. 

```{r, Amend the Clusters, echo=T}

updated.haplotype.df = updated.haplotype.df %>%
  mutate(Haplotype = case_when(Haplotype %in% c("cluster 14") ~ "cluster 10",
                               TRUE ~ Haplotype)) 

```

I'll write out these newly phased clusters of mutations for further analysis. 

```{r, Write out, echo=F}

# Write these out to a file 
print(paste("Writing the results to", output.path))
write_csv(updated.haplotype.df, output.path)

```

## Visualize All Mutations

At this stage in the analysis, we've phased mutations into rough haplotypes. We've also assigned some mutations to the major 'genotypes' (`Genome 1` and `Genome 2`) that were missed in the previous notebook due to low coverage. 

Let's visualize all of these mutations in each tissue specimen. Mutations are annotated by whether they're `Genome 1`, `Genome 2`, or a 'Driver' mutations. These 'Drivers' are mutations with impacts similar to mutations that have been observed in previous cases of SSPE. 

```{r Figure 5, fig.align='center', fig.width=14, fig.height=7, echo=FALSE}

# For simplicity, rename the data frame
labeled.df = updated.haplotype.df

# Re-annotate with genome annotations
# Import genome annotations
annotations.df = read_csv(annotations.filepath, show_col_types = FALSE)

# Get the interval positions for these genes 
N = annotations.df %>% dplyr::filter(Locus == "N")
P.V.C = annotations.df %>% dplyr::filter(`Protein Name` == "phosphoprotein")
M = annotations.df %>% dplyr::filter(Locus == "M")
F. = annotations.df %>% dplyr::filter(Locus == "F")
H = annotations.df %>% dplyr::filter(Locus == "H")
L = annotations.df %>% dplyr::filter(Locus == "L")

# Add gene names to the variants
labeled.df = labeled.df %>% 
  mutate(Gene = case_when(POS < N$Start ~ "3'UTR",
                          POS > L$Stop ~ "5'UTR", 
                          TRUE ~ Gene_Name)) 

# Order of the tissues for the x-axis
expanded_tissue_order = c(
  "SSPE 1", 
  "SSPE 2",
  "Frontal Cortex 2", 
  "Frontal Cortex 1", 
  "Frontal Cortex 3", 
  "Parietal Lobe",
  "Temporal Lobe", 
  "Occipital Lobe",
  "Hippocampus",
  "Internal Capsule", 
  "Midbrain", 
  "Upper Brain Stem",
  "Brain Stem",
  "Cerebellum",
  "Cerebellum Nucleus"
)

# Format the data for plotting
figure_5.df = labeled.df %>% 
  mutate(Genotype = case_when(
    Haplotype == "genome-1" ~ "G-1", 
    Haplotype == "genome-1-1" ~ "G-1",
    Haplotype == "genome-2" ~ "G-2",
    TRUE ~ "subclonal"
  )) %>% 
  select(POS, REF, ALT, SNP, AF, Genotype, Tissue) %>% 
  mutate(Size = if_else(Genotype == "subclonal", "small", "normal")) %>% 
  mutate(Driver = case_when(
    SNP == "G3812A" ~ "M-W125*",
    SNP == "T6817A" ~ "F-L454M",
    SNP == "T3586C" ~ "M-F50S",
    # SNP == "T7293C" ~ "H-I8T",
    SNP == "C7036T" ~ "F-Q527*",
    TRUE ~ ""
  )) %>% 
  mutate(Tissue = case_when(
    Tissue == "UBS" ~ "Upper Brain Stem",
    TRUE ~ Tissue
  ))

core.snps = figure_5.df %>% 
  filter(Genotype != 'subclonal' & Driver == "")

subclonal.snps = figure_5.df %>% 
  filter(Genotype == 'subclonal' & Driver == "")

driver.snps = figure_5.df %>% 
  filter(Driver != "")

# Make figure 4b of the paper
ggplot(data = core.snps, aes(x = factor(Tissue, levels = (expanded_tissue_order)), y = AF, color = Genotype)) +
  geom_jitter(data = subclonal.snps, color = "#878484", size = 1.5, width = 0.25, alpha = 0.75) + 
  geom_jitter(size = 2, width = 0.2, alpha = 1) + 
  geom_jitter(data = driver.snps, aes(shape = Driver, fill=Genotype), color = "black", size = 3, alpha=1, width = 0, stroke = 2) + 
  xlab("Specimen") + 
  ylab("Allele Frequency") +
  scale_color_manual(values = c("#307EC2", "#B63F3E")) +
  scale_fill_manual(values = c("#307EC2", "#878484"), guide = "none") +
  scale_shape_manual(values=c(21,22,24))+
  theme_bw(20) +
  theme(axis.text.x = element_text(angle = 35, hjust=1)) 
  
ggsave(paste(figure.dir, "figure-5.tiff"), device = "tiff", width = 14, height = 7, dpi = 300)

```

What are the SNPs we're calling 'Genome 1 (G1)'. These will be different than the 'Candidate Genome 1' SNPs in Figure 2, because these were identified by their correlation in frequency between all tissue specimens. 

## What about Genome-1-1? 

In the previous notebook, we noticed a set of mutations that made up a haplotype that arises on the background of `Genome 1`, but it mostly missing from the Frontal Cortex 2 sample. We called this haplotype `Genome-1-1`, and it's very significant because it's indicative some early divergence in the frontal cortex. Let's take a look at this haplotype in every tissue.

```{r Figure 6 Without G-FC2, fig.align='center', fig.width=15, fig.height=7, echo=FALSE}

figure.6.df = labeled.df %>% 
  mutate(Genotype = case_when(
    Haplotype == "genome-1" ~ "G-01", 
    Haplotype == "genome-1-1" ~ "G-01 subcluster",
    Haplotype == "genome-2" ~ "G-2",
    TRUE ~ "subclonal"
  )) %>% 
  mutate(Tissue = case_when(
    Tissue == "UBS" ~ "Upper Brain Stem",
    TRUE ~ Tissue
  )) %>% 
  select(Tissue, AF, Genotype, POS) %>% 
  # There is one SNP missing only from the FC2 due to coverage, add this back in with AF as 0
  add_row(Tissue = "Frontal Cortex 2", AF=0, Genotype = "G-01 subcluster", POS = 15084)

subclonal.background = figure.6.df %>% 
  filter(Genotype == "subclonal")

main.genotypes = figure.6.df %>% 
  filter(Genotype != "subclonal")

main.genotypes %>% 
  ggplot(aes(x = POS, y = AF, col = Genotype)) +
    geom_point(data = subclonal.background, aes(x = POS, y = AF), col = "grey",  size = 1) + 
    geom_line() +
    geom_point(size = 2) + 
    xlab("Position") + 
    ylab("Allele Frequency") +
    facet_wrap(~factor(Tissue, levels = (expanded_tissue_order)), ncol=5)+ 
    scale_color_manual(values = c("#87CEEB", "#00008B", "#B63F3E")) +
    theme_bw(18) +
    theme(legend.position="bottom")
  

```