---
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/intersect_properties/',
  fig.width = 6, fig.height = 3.5, dpi = 350,
  fig.path='figure/3-get_intersect_properties/'
)
```
 
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)
)


setwd(paste0(Sys.getenv("HOME"), "/code/malawi_transposon/script/R_SCRIPT/"))

source("modules/frequently_used_items.R")

##########
bb_name <- "astCal1.2_ensembl"
bb_name_abrv <- "astCal"

path_to_intersect <- glue(Sys.getenv("HOME"),
  "/code/malawi_transposon/cloud/workflow/te_permutation_shuffle/results/{bb_name_abrv}-intersect–te_as_whole.bed")
path_to_consensus_sizes <- glue(Sys.getenv("HOME"),
  "/code/malawi_transposon/cloud/repeatmodeler/calliptera_repbase_jun21/{bb_name}-families.fa.fai")
```

## Data import
```{r}
df_intersect <- readr::read_delim(
    path_to_intersect, 
    delim='\t', col_types='ciicd_ciici_i',
    col_names=unlist(strsplit("te_chr,te_start,te_end,te_id,te_age,sv_chr,sv_start,sv_end,sv_id,sv_len,overlap_len", ','))
  ) %>% 
  separate_wider_delim(te_id, delim=':', names=c('te_class','te_subclass','te_id','repmask_id') ) %>% 
  filter(!te_class %in% c('Simple_repeat', 'Low_complexity', 'tRNA', 'snRNA', 'rRNA', 'Satellite')) %>% 
  mutate(te_class = if_else(te_class == 'RC', 'Helitron', te_class)) %>% 
  mutate(te_class = if_else(te_class == 'SINE?', 'SINE', te_class))
  #mutate(te_class = factor(tmp_df$repeat_class_broad, levels = color_schema_te$te_family))

df_consensus_sizes <- readr::read_delim(
    path_to_consensus_sizes, delim='\t', col_type='ci___', col_names=c('repmask_id','consensus_length') 
  ) %>% 
  mutate(repmask_id = gsub("#.+", "", repmask_id) )
```

Add information about consensus family size and calculate statistics we need.
```{r}
df_intersect2 <- df_intersect %>% 
  left_join(df_consensus_sizes) %>% 
  mutate(te_len = te_end - te_start) %>% 
  mutate(
    overlap_perc_te = overlap_len / te_len,  # TE fragment inside the intersect
    overlap_perc_sv = overlap_len / sv_len,  # TE fragment inside the intersect
    intactness = te_len / consensus_length
  ) %>% 
  mutate(
    perfect_intersect = overlap_perc_te > 0.9 & overlap_perc_sv > 0.9,
    core_or_flex = if_else(overlap_perc_te > 0.9, "flexible", "core")
  )
```

## Plots
Count number of perfect intersections.
```{r}
df_perfect <- df_intersect2 %>% filter(perfect_intersect == TRUE)

nrow(df_perfect)

```


```{r}
glue(
  "Number of perfect intersections: {n} ({perc}%), measuring {n_bp} bp",
  perfect_intersections = df_intersect2 %>% filter(perfect_intersect == TRUE),
  n = nrow(perfect_intersections),
  n_bp = sum(perfect_intersections$te_end - perfect_intersections$te_start),
  perc = round(n / nrow(df_intersect2) * 100, 3)
)
```

Number of TE fragments in core and flexible regions.
```{r core_and_flex-n_fragments}
df_intersect2 %>% 
  filter(repmask_id != "rnd-1_family-37") %>% 
  ggplot(aes(te_age, fill = te_class)) +
  geom_histogram(breaks = seq(0, 40, 1) ) +
  facet_wrap(.~core_or_flex, scales = "free_y") +
  geom_vline(xintercept = 10, linetype = 'dotted') +
  xlab("Divergence from ancestral consensus (%)") +
  ylab("Number of fragments") +
  scale_fill_manual('Transposon Class', values = color_schema_te$color, breaks = color_schema_te$te_family) +
  theme(
    legend.position = 'top',
    legend.title = element_text(size = rel(0.8)),
    legend.text = element_text(size = rel(0.7)),
    legend.key.size=unit(0.03, "npc")
  )
```

Cumulative length of TE fragments in core and flexible regions.
```{r core_and_flex-n_bases}
df_intersect2 %>% 
  filter(repmask_id != "rnd-1_family-37") %>% 
  ggplot(aes(te_age, te_len, fill = te_class)) +
  stat_summary_bin(geom='bar', fun='sum', position='stack', breaks = seq(0, 40, 0.5)) +
  facet_wrap(.~core_or_flex, scales = "free_y") +
  geom_vline(xintercept = 10, linetype = 'dotted') +
  xlab("Divergence from ancestral consensus (%)") +
  ylab("Cumulative number of bases") +
  scale_fill_manual('Transposon Class', values = color_schema_te$color, breaks = color_schema_te$te_family) +
  theme(
    legend.position = 'top',
    legend.title = element_text(size = rel(0.8)),
    legend.text = element_text(size = rel(0.7)),
    legend.key.size=unit(0.03, "npc")
  )
```

Intactness measure (not reliable).
```{r core_and_flex-intactness}
df_intersect2 %>% 
  mutate(intactness = if_else(intactness > 1, 1, intactness)) %>% 
  ggplot(aes(x = intactness, fill = te_class)) +
  geom_histogram() +
  facet_wrap(.~core_or_flex, scales = "free_y") +
  xlab("Length of TE fragment / ancestral consensus") +
  ylab("Number of fragments") +
  scale_fill_manual('Transposon Class', values = color_schema_te$color, breaks = color_schema_te$te_family) +
  theme(
    legend.position = 'top',
    legend.title = element_text(size = rel(0.8)),
    legend.text = element_text(size = rel(0.7)),
    legend.key.size=unit(0.03, "npc")
  )
```