---
title: "Density plots across genome tracks"
output: html_document
---
  
```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE, message = FALSE, warning = FALSE,
  cache = TRUE, cache.path = '.cache_knitr/genome_tracks/',
  fig.width = 7, fig.height = 4, dpi = 350,
  fig.path='figure/5b-genome_tracks/',
  dev = c('png', 'pdf')
)
```

In this document, we make some plots to illustrate changes in graph size and pangenome for a wider context of species.

## Setup
```{r, cache=FALSE}
library(GenomicRanges)
library(scales)
library(dplyr)
library(tidyr)
library(ggplot2)
library(tibble)
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")
source("modules/get_chrom_sizes.R")
source("modules/import_transposons_into_granges.R")

# color schema for species
color_schema_species <- readr::read_csv(
  paste0(Sys.getenv("HOME"), 
         "/code/malawi_transposon/metadata/species_colorscheme.csv"), comment='#') %>% 
  mutate(color = paste0('#', color))

# chromosome sizes for calliptera
chrom_sizes <- read_delim(
  paste0(Sys.getenv("HOME"), "/igv/genomes/astCal1.2_ensembl/chrom.sizes"),
  col_names = c('chr', 'len')
)
chrom_sizes <- Seqinfo(chrom_sizes$chr, chrom_sizes$len)

# window sizes should be different based of features to aid resolution
window_size_cov <- 10**6 / 4
window_size_sv  <- 10**5 / 4
window_size_te  <- 10**5 / 4
```

## Data import
Import coverage data, as well as coordinates of SVs, TEs and genes.
```{r}
# import file with coordinates and coverage of backbone segments
file_dir <- paste0(Sys.getenv("HOME"), "/code/malawi_transposon/cloud/workflow/wider_context_comparisons/results/")
df_bb <- read_csv(paste0(file_dir, "malawi-segments.bb"))

# bubble coordinates
df_sv <- read_csv(
    "./local/malawi_haplochromines-summary-preprocessed.csv",
    col_types = paste0('cciiiiicciiicilliidc', 
                       paste0(rep('ciicic', num_of_genomes), collapse=''), collapse=''),
    na = c("", "NA", ".") ) %>% 
  select(id, chr, start, end)

gr_sv <- GRanges(
  seqnames = df_sv$chr, 
  ranges = IRanges(start = df_sv$start, end = df_sv$end), 
  id = df_sv$id, seqinfo = chrom_sizes
)

# TE coordinates
gr_te <- import_transposons_into_granges("astCal", annoType="calliptera")
```


## Coverage
Here on, we are plotting the coverage by genome position for each chromosome. This is done in a windowed approach, where we aggregate the number of bases without coverage. Flexible segments whose source and sink are covered above a threshold are given coverage of 1, this was done in some preprocessing.
```{r}
# calculate windows
bb_windows <- df_bb %>%
  select(segment, cum_start, cum_end) %>% 
  mutate(
    lower = floor(cum_start/window_size_cov),
    upper = ceiling(cum_end/window_size_cov)
  ) %>% 
  mutate(window = mapply(
      FUN = function(x,y) paste0(seq(x,y-1,1), collapse=","),
      x = lower, y = upper
    )
  ) %>% 
  separate_longer_delim(window, delim = ",") %>% 
  select(-lower, -upper)
  
# recalculate lower and upper boundaries based on defined windows
df_bb_windowed <- bb_windows %>% 
  mutate(
    start = as.integer(window) * window_size_cov + 1, 
    end = (as.integer(window)+1) * window_size_cov
  ) %>% 
  mutate(
    bp_in_window = if_else(cum_end > end, end, cum_end) - if_else(cum_start > start, cum_start, start)
  ) %>% 
  left_join(select(df_bb, segment, chr_name, starts_with("cov_")))

# aggregate over windows and species / type of coverage
df_bb_windowed_agg <- df_bb_windowed %>% 
  mutate(chr_name = gsub("astCal_", "", chr_name)) %>% 
  pivot_longer(starts_with("cov_"), names_to = "cov_type", names_pattern = "cov_(.*)") %>% 
  group_by(chr_name, start, cov_type) %>% 
  summarise(bases_covered = sum(value*bp_in_window), bases_total = sum(bp_in_window)) %>% 
  mutate(
    cov_perc = bases_covered / bases_total * 100,
    chr_name = factor(chr_name, levels = seqnames(chrom_sizes))
  )
```

Cumulative and ubiquitous coverage.
```{r tracks_coverage_cuml_ubiq, fig.width=7, fig.height=9}
df_bb_windowed_agg %>% 
  filter(cov_type %in% c("ubiq", "cuml")) %>% 
  filter(!stringr::str_detect(chr_name, "^L")) %>% 
  ggplot(aes(x=start, y=cov_perc, fill=cov_type)) +
  geom_area(position="identity") +
  geom_line(aes(x=start, y=0), color='black', linewidth=rel(0.6)) +
  facet_wrap(.~chr_name, ncol = 1, strip.position = "left") +
  scale_fill_manual(values=c("#92D050", "#018571")) +
  scale_x_continuous(expand = expansion(mult = c(.02, .05))) +
  ylim(-10, 105) +
  theme_void() +
  theme(legend.position = "none", strip.text.y = element_text(size = rel(1.6)))
```

Coverage of individual species.
```{r tracks_coverage_individual, fig.width=7, fig.height=9}
df_bb_windowed_agg %>% 
  filter(!cov_type %in% c("ubiq", "cuml")) %>% 
  filter(!stringr::str_detect(chr_name, "^L")) %>% 
  ggplot(aes(x=start, y=cov_perc, color=cov_type)) +
  geom_area(
    aes(y = max_cov_perc),
    data = . %>% group_by(chr_name, start) %>% summarise(max_cov_perc = max(cov_perc)),
    position="identity", color='white', fill='gray85'
  ) +
  geom_line(position="identity", linewidth=rel(0.4)) +
  geom_line(aes(x=start, y=0), color='black', linewidth=rel(0.6)) +
  facet_wrap(.~chr_name, ncol = 1, strip.position = "left") +
  scale_color_manual("sample", values=color_schema_species$color[2:8], labels=color_schema_species$species_abrv[2:8]) +
  scale_x_continuous(expand = expansion(mult = c(.02, .05))) +
  ylim(-5, 105) +
  theme_void() +
  guides(colour = guide_legend(nrow = 1)) +
  theme(legend.position = "top", 
        legend.title = element_text(size=rel(0.7)), legend.text=element_text(size=rel(0.7)),
        strip.text.y = element_text(size = rel(1.6)) )
```


## Bubble density
Taking the GRanges for the bubbles, and trim it down to the midpoint.
```{r}
gr_sv_resized <- resize(gr_sv, width=1, fix="center")
```

Create GRanges for the final windows, trimming as necessary.
```{r}
gr_windows <- GRanges(
  seqnames = seqnames(chrom_sizes),
  ranges = IRanges(start=rep(1, length(seqlengths(chrom_sizes))), end=seqlengths(chrom_sizes)),
  seqinfo = chrom_sizes) %>%
  tile(width = window_size_sv) %>% 
  unlist() 
```

Intersect the GRanges to obtain the number of bubbles per window
```{r}
idx <- findOverlaps(gr_windows, gr_sv_resized) %>% 
  as.data.frame() %>% 
  as_tibble() %>% 
  group_by(queryHits) %>% 
  summarise(n_bubbles = n())

# add this info into the GRanges
gr_windows$n_bubble <- 0
gr_windows$n_bubble[idx$queryHits] <- idx$n_bubbles

# convert to dataframe
df_sv_density <- as.data.frame(gr_windows) %>% 
  as_tibble() %>% 
  rename(chr_name = seqnames) %>% 
  select(chr_name, start, width, n_bubble)
```

This plot shows the bubble density calculated with a window size of `r window_size_sv` bp.
```{r tracks_bubble_density, fig.width=7, fig.height=9}
df_sv_density %>% 
  filter(!stringr::str_detect(chr_name, "^L")) %>% 
  ggplot(aes(x=start,y=n_bubble)) +
  geom_line(color='#b86500') +
  geom_line(aes(x=start, y=0), color='black', linewidth=rel(0.6)) +
  facet_wrap(.~chr_name, ncol = 1, strip.position = "left") +
  scale_x_continuous(expand = expansion(mult = c(.02, .05))) +
  theme_void() +
  theme(legend.position = "none", strip.text.y = element_text(size = rel(1.6)))
```


## TE density
Taking the GRanges for the transposons, trim it down to midpoint.
```{r}
gr_te_resized <- resize(gr_te, width=1, fix="center")
```

Create GRanges for the final windows, trimming as necessary
```{r}
gr_windows <- GRanges(
  seqnames = seqnames(chrom_sizes),
  ranges = IRanges(start=rep(1, length(seqlengths(chrom_sizes))), end=seqlengths(chrom_sizes)),
  seqinfo = chrom_sizes) %>%
  tile(width = window_size_te) %>% 
  unlist() 
```

Intersect the GRanges to obtain the number of TEs per window
```{r}
idx <- findOverlaps(gr_windows, gr_te_resized) %>% 
  as.data.frame() %>% 
  as_tibble() %>% 
  group_by(queryHits) %>% 
  summarise(n_te = n())

gr_windows$n_te <- 0
gr_windows$n_te[idx$queryHits] <- idx$n_te

df_te_density <- as.data.frame(gr_windows) %>% 
  as_tibble() %>% 
  rename(chr_name = seqnames) %>% 
  select(chr_name, start, width, n_te)
```

This plot shows the bubble density calculated with a window size of `r window_size_sv` bp.
```{r tracks_te_density, fig.width=7, fig.height=9}
df_te_density %>% 
  filter(!stringr::str_detect(chr_name, "^L")) %>% 
  ggplot(aes(x=start, y=n_te)) +
  geom_line(color='#5895db') +
  geom_line(aes(x=start, y=0), color='black', linewidth=rel(0.6)) +
  facet_wrap(.~chr_name, ncol = 1, strip.position = "left") +
  scale_x_continuous(expand = expansion(mult = c(.02, .05))) +
  theme_void() +
  theme(legend.position = "none", strip.text.y = element_text(size = rel(1.6)))
```


## Overlaying genomic tracks
```{r tracks_overlay_sv_coverage, fig.width=7, fig.height=9}
tmp_df1 <- df_bb_windowed_agg %>% filter(cov_type == "cuml", !stringr::str_detect(chr_name, "^L"))
tmp_df2 <- df_sv_density %>% 
  filter(!stringr::str_detect(chr_name, "^L")) %>% 
  mutate(n_bubble_scaled = n_bubble / max(n_bubble) * 100)

tmp_df1 %>% 
  ggplot(aes(x=start, y=cov_perc)) +
  geom_area(position="identity", fill="#92D050") +
  facet_wrap(.~chr_name, ncol = 1, strip.position = "left") +
  geom_line(aes(x=start,y=n_bubble_scaled), tmp_df2, color="#b86500") +
  geom_line(aes(x=start, y=0), color='black', linewidth=rel(0.6)) +
  scale_x_continuous(expand = expansion(mult = c(.02, .05))) +
  ylim(-10, 105) +
  theme_void() +
  theme(legend.position = "none", strip.text.y = element_text(size = rel(1.6)))

rm(tmp_df1, tmp_df2)
```

