---
title: "Getting intersect properties"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE, message = FALSE, warning = FALSE,
  cache = TRUE, cache.path = '.cache_knitr/permutation_plots/',
  fig.width = 6, fig.height = 3.5, dpi = 350,
  fig.path='figure/3b-permutation/',
  dev = c('png', 'pdf')
)
```
 
In this document, we calculate a few more statistics about the intersecting TEs to see what things are worth doing a permutation shuffle for.

I've also written this script to be flexible in case we need to rerun this with a different backbone.

## Setup
```{r, cache=FALSE}
library(scales)
library(factoextra)
library(ggplot2)
library(tibble)
library(dplyr)
library(tidyr)
library(readr)
library(glue)

theme_set(theme_classic())
theme_update(
  panel.background = element_rect(fill = "transparent", colour = NA),
  plot.background = element_rect(fill = "transparent", colour = NA),
  legend.background = element_rect(fill = "transparent", colour = NA)
)


source("modules/frequently_used_items.R")
```

## Data import
```{r}
filepath <- paste0(Sys.getenv("HOME"), "/code/malawi_transposon/script/te_permutation_shuffle/")

df <- read_csv(paste0(filepath, "results/astCal-stats_flex.csv"))

df_genome <- read_csv(paste0(filepath, "results/astCal-stats_genomewide.csv")) %>% 
  mutate(
    perm = "genomewide",
    te_class = factor(te_class, levels = c('All Transposons', 'All Transposons except Unknown', color_schema_te$te_family))
  )

df_flex <- filter(df, perm == 'default') %>% 
  mutate(
    perm = "flex", 
    te_class = factor(te_class, levels = c('All Transposons', 'All Transposons except Unknown', color_schema_te$te_family))
  )

df_perm <- filter(df, perm != 'default') %>% 
  mutate(te_class = factor(te_class, levels = c('All Transposons', 'All Transposons except Unknown', color_schema_te$te_family)) ) %>% 
  left_join(df_flex %>% select(te_class, age_10_flex = age_10, age_50_flex = age_50) )
```

## Summary statistics
```{r}
df_genome %>% mutate(across(where(is.double), ~ round(., 2)) ) %>% as.data.frame()

df_flex %>% mutate(across(where(is.double), ~ round(., 2)) ) %>% as.data.frame()

df_perm %>% 
  group_by(te_class) %>% 
  summarise(
    age_10 = round(mean(age_10), 2), 
    age_50 = round(mean(age_50), 2),
    pvalue10 = (sum(age_10 > age_10_flex) + 1) / n(),
    pvalue50 = (sum(age_10 > age_50_flex) + 1) / n(),
  ) %>% 
  as.data.frame()
```

## Permutation testing plots
```{r permutation_histograms, fig.height=4, fig.width=7}
df_perm %>% 
  ggplot(aes(age_10, fill = te_class)) +
  geom_histogram(bins=80) +
  facet_wrap(.~te_class, scales='free', ncol=3) +
  geom_point(aes(x=age_10, y=0), data=df_flex, shape=24) +
  geom_point(aes(x=age_10, y=0), data=df_genome, shape=21) +
  scale_fill_manual('Transposon Class', 
    values = c('gray70', 'gray50', color_schema_te$color_dark) ) +
  xlab("Percentage sequence of TE class") +
  ylab("Frequency") +
  theme(legend.position = 'none')

df_perm %>% 
  ggplot(aes(age_50, fill = te_class)) +
  geom_histogram(bins=80) +
  facet_wrap(.~te_class, scales='free', ncol=3) +
  geom_point(aes(x=age_50, y=0), data=df_flex, shape=24) +
  geom_point(aes(x=age_50, y=0), data=df_genome, shape=21) +
  scale_fill_manual('Transposon Class', 
    values = c('gray70', 'gray50', color_schema_te$color_dark) ) +
  xlab("Percentage sequence of TE class") +
  ylab("Frequency") +
  theme(legend.position = 'none')
```