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

In this document, we make some plots for the presence-absence variation analysis done with ODGI. We do this PAV for the entire gene body, exons only, upstream (-2000 to +500bp of gene start) and downstream (-500 to +2000bp of gene end).

PAV data has been precalculated with `odgi`, but the first part of this notebook focuses on identifying which genes may not be reliable for this kind of analysis.

We need to restrict the analysis to genes where the gene body plus 2000bp flanking regions are located in regions which have coverage, and are not located in complex structural variation. We also need to do a preliminary filtering step to drop genes which are misannotated as transposable elements.

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

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/get_chrom_sizes.R")
source("modules/import_transposons_into_granges.R")
source("modules/misannotation_thresholds.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))

##########
subdir <- "malawi_haplochromines"
bb_name <- "astCal1.2_ensembl"
bb_name_sqlite <- "astCal_v1.2"
bb_name_abrv <- "astCal"
chr1_id <- "1"
```

## Data import
Import coordinates for genes on the backbone.
```{r}
# get coordinates for genes
path_to_annotation_sqlite <- glue(
  Sys.getenv("HOME"), "/igv/genomes/{bb_name}/{bb_name_sqlite}.gff.sqlite"
)
txdb <- AnnotationDbi::loadDb(file = path_to_annotation_sqlite)
gr_gene <- genes(txdb, columns = c('gene_id', 'tx_name'))
```

Import coverage values for each graph segment. **Should consider plotting histogram for the core segments?**
```{r}
df_coverage <- read_csv(
  glue(Sys.getenv("HOME"), 
    "/code/malawi_transposon/cloud/graph_genome/{subdir}/pseudoreference/pav/graph_bbsegments_coverage.csv"),
  col_types = 'cciiliiiiiiiiii')
df_coverage
```

Import coordinates of bubbles. We do some preprocessing to remove some bubbles (note that the criteria is slightly different from the `preprocess.Rmd`) and add columns required for the filtering of complex bubbles.
```{r}
df_sv <- read_csv(
  glue(Sys.getenv("HOME"), "/code/malawi_transposon/cloud/graph_genome/{subdir}/{subdir}-summary.csv"),
  col_types = paste0('ciiiiiiic', paste0(rep('ciicic', num_of_genomes), collapse=''), collapse=''),
  na = c("", "NA", ".") ) %>% 
  mutate(id=paste0(chr, ':', start, '-', end)) %>% # add identifier column 
  filter(n_paths != 1) %>% # drop bubbles with only one paths 
  relocate(id)

# count number of uncalled alleles
df_sv['n_uncalled'] <- df_sv %>% select(ends_with('_path')) %>% is.na() %>% rowSums() 

# remove bubbles where everything had the same allele
df_sv$allSameAlleleBool <- (
  rowSums(select(df_sv, ends_with('_path')) == df_sv[[glue("{bb_name}_path")]], na.rm = TRUE) 
  == (num_of_genomes - df_sv$n_uncalled)
)
df_sv <- filter(df_sv, allSameAlleleBool == FALSE)

# number of biological paths, as well as the shortest and longest lengths
df_sv <- left_join(df_sv, 
  df_sv %>% 
    select(id, n_paths, ends_with(c('_path', '_len'))) %>% 
    pivot_longer(cols = ends_with(c('_path', '_len')),
                 names_to = c('species', '.value'),
                 names_pattern = "(.*)_(.*)") %>% 
    drop_na(path, len) %>%   # some species have no variants called
    group_by(id, path, len) %>% 
    count(path) %>% 
    group_by(id) %>% 
    summarise(
      n_paths_bio = n(), lenbio_min = min(len), lenbio_max = max(len)
    )
)

df_sv
```


## Preprocessing gene coordinates
### Transposable element and gene description
First we calculate how much the TE multi and single overlaps agree with each other.
```{r misannotation_histogram, fig.width = 3.2, fig.height = 2.8}
df_gene_trueness <- readr::read_csv(
  glue(Sys.getenv("HOME"), "/code/malawi_transposon/cloud/genome/{bb_name}/gene_trueness.csv")
)

ggplot(df_gene_trueness, aes(maxOverlapSingleTE*100, fill = totalOverlapMultipleTE>threshGeneBodyOverlapMultipleTE)) +
  geom_histogram(bins=50, color='black') +
  xlab("% overlap with single TE fragment") +
  ylab("Number of genes") +
  coord_cartesian(ylim=c(0, 8500)) +
  scale_fill_manual(glue("overlap w\nmultiple > {threshGeneBodyOverlapMultipleTE*100}%"), values=c('gray80', 'gray20')) +
  theme(legend.position = 'top', legend.title=element_text(size=rel(0.8)), legend.text=element_text(size=rel(0.7)) )

ggplot(df_gene_trueness, aes(totalOverlapMultipleTE*100, fill = maxOverlapSingleTE>threshGeneBodyOverlapSingleTE)) +
  geom_histogram(bins=50, color='black') +
  xlab("% overlap with multiple TE fragments") +
  ylab("Number of genes") +
  coord_cartesian(ylim=c(0, 8500)) +
  scale_fill_manual(glue("overlap w\nsingle > {threshGeneBodyOverlapSingleTE*100}%"), values=c('gray80', 'gray20')) +
  theme(legend.position = 'top', legend.title=element_text(size=rel(0.8)), legend.text=element_text(size=rel(0.7)) )
```

We need to do filtering. Genes that have an associated description are kept, and for those that done, we remove them if one of the following apply:

* Overlap more than **`r threshGeneBodyOverlapSingleTE*100`%** across its gene body with a single TE.
* Overlap more than **`r threshGeneBodyOverlapMultipleTE*100`%** in total with any TEs.
* Associated with the Gene Ontology term "DNA transposition"

```{r}
list_of_dodgy_genes <- df_gene_trueness %>% filter(
    (maxOverlapSingleTE > threshGeneBodyOverlapSingleTE | 
       totalOverlapMultipleTE > threshGeneBodyOverlapMultipleTE |
       hasTranspositionTerm == TRUE ) 
  ) %>% .$gene_id

glue("{length(x)} genes are potential TEs to be removed", x = list_of_dodgy_genes)
```

Venn diagram of the overlap between these.
```{r venn_filter, fig.width=4.2, fig.height=2.2}
set1 <- df_gene_trueness %>% filter(maxOverlapSingleTE > threshGeneBodyOverlapSingleTE) %>% .$gene_id
set2 <- df_gene_trueness %>% filter(totalOverlapMultipleTE > threshGeneBodyOverlapMultipleTE) %>% .$gene_id
total_pop_size <- length(unique(c(set1, set2)))

library(VennDiagram)
futile.logger::flog.threshold(futile.logger::ERROR, name = "VennDiagramLogger")
grid.newpage()
venn_object <- venn.diagram(
  x = list(set1, set2), filename=NULL, 
  category.names = c("single" , "multiple"),
  hyper.test=TRUE, total.population = total_pop_size,
  output=FALSE
)
grid.draw(venn_object)
```

### Coverage
Create a GRanges object to represent segments where **at least 6 out of 8** samples have a coverage above the predefined threshold of **`r threshMinCoverageAcrossSpecies*100`%**. This should be similar to be the **ubiquitous** coverage calculated in previous notebooks (maybe even higher because it does not weigh things by coverage).
```{r}
get_GRanges_for_covered_segments <- function(cov_threshold = 80, n_samples_min = 6) {
  chosen_cov_threshold <- paste0("n_cov", cov_threshold)
  tmp_df <- df_coverage %>% filter(!!as.name(chosen_cov_threshold) >= n_samples_min)
  gr_covered <- GenomicRanges::GRanges(
    seqnames = tmp_df$chr,
    ranges = IRanges(start = tmp_df$start, end = tmp_df$end),
    segment = tmp_df$segment,
    seqinfo = Seqinfo(chrom_sizes_list[[bb_name_abrv]][[1]], chrom_sizes_list[[bb_name_abrv]][[2]]) ) %>% 
    GenomicRanges::trim()
  return(gr_covered)
}

gr_covered <- get_GRanges_for_covered_segments(threshMinCoverageAcrossSpecies * 100)
glue(round(sum(width(gr_covered)) / sum(seqlengths(gr_covered)) * 100, 3), "% of the backbone is covered in 6 of the species.")
```

Here we wish to retrieve a list of genes which are almost entirely in backbone segments with high degree coverage.

* Extend the genes upstream and downstream by **50 bp**.
* Calculate the percentage sequence of each gene + flanking regions located inside these high coverage segments.
* Keep genes overlapping above **`r threshGeneBodyOverlapHighCovSegments*100`%**.

*Note: originally, I flanked both sides of the gene with 2000~bp for the PAV with upstream and downstream regions. However, I decided that doing a PAV for those regions can be quite dodgy, because they tend to overlap a lot into other genes. I still put a little small flanking region equal to half the chop size of ODGI, because I think that can help.*

```{r}
# add flanking regions
gr_gene_flanked <- resize(gr_gene, width = width(gr_gene) + 50, fix="start") %>% 
  resize(width = width(.) + 50, fix="end")

# overlap gene with covered segments
hits <- findOverlaps(gr_gene_flanked, gr_covered, ignore.strand=TRUE)
overlaps <- pintersect(gr_gene_flanked[queryHits(hits)], gr_covered[subjectHits(hits)], ignore.strand=TRUE)

# turn into dataframe to do a grouped by sum of coverage
df_percentGeneOverlap_coverage <- tibble(
  gene = gr_gene_flanked[queryHits(hits)]$gene_id,
  uncovered_range_idx = subjectHits(hits),
  percentGeneOverlap = width(overlaps) / width(gr_gene_flanked[queryHits(hits)]) ) %>% 
  group_by(gene) %>% 
  summarise(percentGeneOverlap = sum(percentGeneOverlap))

# filter through genes that overlap more than that threshold
geneList_goodCoverage <- df_percentGeneOverlap_coverage %>% 
  filter(percentGeneOverlap > threshGeneBodyOverlapHighCovSegments) %>% 
  .$gene %>% unique()

glue(length(geneList_goodCoverage), 
  " ({percent_covered}%) of genes overlap significantly with segments with good coverage",
   percent_covered = round(length(geneList_goodCoverage) / length(gr_gene) * 100, 2)
)
```

We also define another set of genes with a stricter coverage threshold. For certain statistics, we want to restrict this analysis to only the genes that actually have decent coverage across all the samples.
```{r}
gr_covered_Strict <- get_GRanges_for_covered_segments(90, n_samples_min = 8)
glue(round(sum(width(gr_covered_Strict)) / sum(seqlengths(gr_covered_Strict)) * 100, 3), "% of the backbone is covered in all the species.")

# overlap gene with covered segments
hits <- findOverlaps(gr_gene_flanked, gr_covered_Strict, ignore.strand=TRUE)
overlaps <- pintersect(gr_gene_flanked[queryHits(hits)], gr_covered_Strict[subjectHits(hits)], ignore.strand=TRUE)

# turn into dataframe to do a grouped by sum of coverage
df_percentGeneOverlap_coverage <- tibble(
  gene = gr_gene_flanked[queryHits(hits)]$gene_id,
  uncovered_range_idx = subjectHits(hits),
  percentGeneOverlap = width(overlaps) / width(gr_gene_flanked[queryHits(hits)]) ) %>% 
  group_by(gene) %>% 
  summarise(percentGeneOverlap = sum(percentGeneOverlap))

# filter through genes that overlap more than that threshold
geneList_goodCoverageStrict <- df_percentGeneOverlap_coverage %>% 
  filter(percentGeneOverlap > threshGeneBodyOverlapHighCovSegments) %>% 
  .$gene %>% unique()

glue(length(geneList_goodCoverageStrict), 
  " ({percent_covered}%) of genes overlap significantly with segments with good coverage",
   percent_covered = round(length(geneList_goodCoverageStrict) / length(gr_gene) * 100, 2)
)
```

Venn diagram to show that this stricter set of genes is in the superset of the looser one.
```{r vene_geneCoverage, fig.width=3.4, fig.height=3.4}
futile.logger::flog.threshold(futile.logger::ERROR, name = "VennDiagramLogger")
grid.newpage()
venn_object <- venn.diagram(
  x = list(geneList_goodCoverage, geneList_goodCoverageStrict),
  category.names = c("good coverage" , "strict coverage"),
  total.population = length(gr_gene),
  output=FALSE, filename=NULL
)
grid.draw(venn_object)
```

### Complex structural variation
In the flexible segments of complex bubbles, the PAV scores of such genes can be a bit unreliable, so we might need to exclude such genes at a later stage.
```{r}
# define new column for simple bubbles (stricter than original)
df_sv$simple2 <- (
  df_sv$n_paths_bio <= 5 &
  df_sv$lenbio_min < 1000 & 
  df_sv$lenbio_min / df_sv$lenbio_max < 0.2
)

# extract list of segments of complex bubbles (core + flexible)
list_of_complex_segments <- df_sv %>% 
  filter(simple2 == FALSE) %>% 
  separate_rows(segment_list) %>% .$segment_list

# extract flexible segments to create GRanges object
get_GRanges_for_complex_flex_segments <- function() {
  tmp_df <- df_coverage %>% 
    filter(segment %in% list_of_complex_segments, core_bool == FALSE) 
  gr_complex_flex <- GenomicRanges::GRanges(
    seqnames = tmp_df$chr,
    ranges = IRanges(start = tmp_df$start, end = tmp_df$end),
    segment = tmp_df$segment,
    seqinfo = Seqinfo(chrom_sizes_list[[bb_name_abrv]][[1]], chrom_sizes_list[[bb_name_abrv]][[2]])
  )
  return(gr_complex_flex)
}

gr_complex_flex <- get_GRanges_for_complex_flex_segments()
glue(round(sum(width(gr_complex_flex)) / sum(seqlengths(gr_complex_flex)) * 100, 3), "% of the backbone is in complex bubbles.")
```

Retrieve the list of genes of which **`r threshGeneBodyOverlapComplexRegion*100`%** of their length + flanking regions are in complex SVs.
```{r}
hits <- findOverlaps(gr_gene_flanked, gr_complex_flex, ignore.strand=TRUE)
overlaps <- pintersect(gr_gene_flanked[queryHits(hits)], gr_complex_flex[subjectHits(hits)], ignore.strand=TRUE)

df_percentGeneOverlap_complex <- tibble(
  gene = gr_gene_flanked[queryHits(hits)]$gene_id,
  uncovered_range_idx = subjectHits(hits),
  percentGeneOverlap = width(overlaps) / width(gr_gene_flanked[queryHits(hits)]) ) %>% 
  group_by(gene) %>% 
  summarise(percentGeneOverlap = sum(percentGeneOverlap))

# get list of genes that overlap more than that threshold with complex regions
geneList_complexSegments <- df_percentGeneOverlap_complex %>% 
  filter(percentGeneOverlap > threshGeneBodyOverlapComplexRegion) %>% 
  .$gene %>% unique()

glue(length(geneList_complexSegments), 
  " ({percent_covered}%) of genes overlap significantly with segments of complex SVs",
   percent_covered = round(length(geneList_complexSegments) / length(gr_gene) * 100, 2)
)
```

## Compiling TE misannotation, coverage and complex bubble information
Create a dataframe of all these information to be joined to PAV dataframes later.
```{r}
df_gene_quality <- tibble(gene = names(gr_gene)) %>% 
  mutate(
    notPutativeTE = !gene %in% list_of_dodgy_genes,
    goodCoverage = gene %in% geneList_goodCoverage,
    notInComplexRegions = !gene %in% geneList_complexSegments
  )

colSums(select(df_gene_quality, -gene))

# final set of reliable genes for us
glue(
  nrow(filter(df_gene_quality, notPutativeTE==T, goodCoverage==T, notInComplexRegions==T)),
  " genes are reliable for PAV analysis"
)
```

Venn Diagram comparing all three sets.
```{r venn_geneforPAV, fig.width=3.4, fig.height=3.4}
futile.logger::flog.threshold(futile.logger::ERROR, name = "VennDiagramLogger")
grid.newpage()
venn_object <- venn.diagram(
  x = list(
    filter(df_gene_quality, notPutativeTE == TRUE) %>% .$gene, 
    filter(df_gene_quality, goodCoverage == TRUE) %>% .$gene,
    filter(df_gene_quality, notInComplexRegions == TRUE) %>% .$gene
  ),
  category.names = c("not TE" , "good coverage", "not in complex segments"),
  total.population = length(gr_gene),
  output=FALSE, filename=NULL,
  fill = brewer.pal(3, "Pastel2")
)
grid.draw(venn_object)
```

## Coverage by species
Add information about covered genes by species into the gene quality data frame.
```{r}
df_covered_genes_by_species <- readr::read_csv(
    glue("./local/segment_coverage_by_species-{bb_name_abrv}.csv")
  ) %>% 
  mutate(covered = TRUE) %>% 
  pivot_wider(id_cols = "gene", names_from = sample, values_from = covered) %>% 
  replace(is.na(.), FALSE)

# quick count
colSums(df_covered_genes_by_species[,-1])
```

## Import PAV data
PAV data has been precalculated with `odgi`. Now we will import them.

PAV for whole genes.
```{r}
df_pav_gene <- read_delim(
  glue(Sys.getenv("HOME"), 
       "/code/malawi_transposon/cloud/graph_genome/{subdir}/pseudoreference/pav/",
       "txdb_gene_by_species.pav"), col_types = 'ciiccd') %>% 
  rename(gene = name, sample = group)
```

PAV for individual exons. This is then used calculate a weighted PAV.
```{r}
df_pav_exon <- read_delim(
  glue(Sys.getenv("HOME"), 
       "/code/malawi_transposon/cloud/graph_genome/{subdir}/pseudoreference/pav/",
       "txdb_exonparts_by_species.pav"), col_types = 'ciiccd') %>% 
  rename(gene = name, sample = group) %>% 
  separate_rows(gene, sep = ',')

# aggregate by species and gene
df_pav_exon_agg <- df_pav_exon %>% 
  mutate(width = end - start) %>% 
  group_by(sample, gene) %>% 
  summarise(covered_bases = sum(pav * width), total_bases = sum(width), n_exons = n()) %>% 
  mutate(pav = covered_bases / total_bases)
```

PAV for upstream and downstream gene regions.
```{r}
df_pav_upstream <- read_delim(
  glue(Sys.getenv("HOME"), 
       "/code/malawi_transposon/cloud/graph_genome/{subdir}/pseudoreference/pav/",
       "txdb_upstream_by_species.pav"), col_types = 'ciiccd') %>% 
  rename(gene = name, sample = group)

df_pav_downstream <- read_delim(
  glue(Sys.getenv("HOME"), 
       "/code/malawi_transposon/cloud/graph_genome/{subdir}/pseudoreference/pav/",
       "txdb_downstream_by_species.pav"), col_types = 'ciiccd') %>% 
  rename(gene = name, sample = group)
```

PAV for intronic regions.
```{r}
df_pav_intron <- read_delim(
  glue(Sys.getenv("HOME"), 
       "/code/malawi_transposon/cloud/graph_genome/{subdir}/pseudoreference/pav/",
       "txdb_intronparts_by_species.pav"), col_types = 'ciiccd') %>% 
  rename(gene = name, sample = group)

df_pav_intron_agg <- df_pav_intron %>% 
  mutate(width = end - start) %>% 
  group_by(sample, gene) %>% 
  summarise(covered_bases = sum(pav * width), total_bases = sum(width), n_exons = n()) %>% 
  mutate(pav = covered_bases / total_bases)
```

Combine all that information into a dataframe, and add information about the quality of the genes or PAV analysis.
```{r}
# combine all data so far
df_pav_master <- select(df_pav_gene, gene, sample, pav_gene = pav) %>% 
  full_join(select(df_pav_exon_agg, gene, sample, pav_exon = pav)) %>% 
  full_join(select(df_pav_upstream, gene, sample, pav_upstream = pav)) %>% 
  full_join(select(df_pav_downstream, gene, sample, pav_downstream = pav)) %>% 
  full_join(select(df_pav_intron_agg, gene, sample, pav_intron = pav)) %>% 
  full_join(df_gene_quality) %>% 
  mutate(geneQualScore = rowSums(select(., notPutativeTE, goodCoverage, notInComplexRegions))) %>%
  mutate(strictCoverage = gene %in% geneList_goodCoverageStrict) %>% 
  filter(sample != 'backbone') %>% 
  left_join(
    df_covered_genes_by_species %>% pivot_longer(cols = astCal:rhaChi2, names_to = "sample", values_to = "covered")
  ) %>% 
  replace_na(list(covered = FALSE))

head(df_pav_master, 10)
```

Pivot to longer format, where each gene + sample + feature combination is a row.
```{r}
df_pav_master_long <- df_pav_master %>% 
  filter(geneQualScore == max(geneQualScore)) %>% 
  select(gene, sample, pav_gene, pav_exon, pav_upstream, pav_downstream, pav_intron) %>% 
  pivot_longer(starts_with("pav_"), names_to='feature', names_pattern = "pav_(.+)", values_to='pav_score') %>% 
  mutate(feature = if_else(feature == 'gene', 'gene body', feature)) %>% 
  mutate(feature = forcats::fct_relevel(feature, 'gene body', 'exon', 'intron', 'upstream', 'downstream'))

head(df_pav_master_long, 10)
```


## PAV: Plots and stats
Here's the **total number of genes that have good quality metrics for PAV analysis**. This number is the same for all the species.
```{r}
highQualGeneCount <- df_pav_master %>% filter(geneQualScore == max(geneQualScore), sample == bb_name_abrv) %>% nrow()
highQualGeneCount

highQualGeneCount_Strict <- df_pav_master %>% filter(geneQualScore == max(geneQualScore), strictCoverage == TRUE, sample == bb_name_abrv) %>% nrow()
highQualGeneCount_Strict
```

This tabulates the number of genes which high PAV scores for each sample, considering all the genes that all possible to do PAV analysis for per species.
```{r}
# gene body
df_pav_master %>% 
  rename(pav_current = pav_gene) %>% 
  filter(notPutativeTE == TRUE, notInComplexRegions == TRUE, covered == TRUE, !is.na(pav_current)) %>% 
  group_by(sample) %>% 
  summarise(
    count = sum(pav_current > 0.9, na.rm = TRUE), 
    total = n(), 
    percent = sum(pav_current > 0.9, na.rm = TRUE) / n() * 100
  ) 

# exon
df_pav_master %>% 
  rename(pav_current = pav_exon) %>% 
  filter(notPutativeTE == TRUE, notInComplexRegions == TRUE, covered == TRUE, !is.na(pav_current)) %>% 
  group_by(sample) %>% 
  summarise(
    count = sum(pav_current > 0.9, na.rm = TRUE), 
    total = n(), 
    percent = sum(pav_current > 0.9, na.rm = TRUE) / n() * 100
  )
```

This tabulates the number of genes which high PAV scores for each sample, depending on which genic feature is considered. For this, we use the **strict coverage** gene list.
```{r}
# gene body
df_pav_master %>% 
  rename(pav_current = pav_gene) %>% 
  filter(geneQualScore == max(geneQualScore), strictCoverage == TRUE, !is.na(pav_current)) %>% 
  group_by(sample) %>% 
  summarise(
    count = sum(pav_current > 0.9, na.rm = TRUE), 
    total = n(), 
    percent = sum(pav_current > 0.9, na.rm = TRUE) / n() * 100
  )

# exons
df_pav_master %>% 
  rename(pav_current = pav_exon) %>% 
  filter(geneQualScore == max(geneQualScore), strictCoverage == TRUE, !is.na(pav_current)) %>% 
  group_by(sample) %>% 
  summarise(
    count = sum(pav_current > 0.9, na.rm = TRUE), 
    total = n(), 
    percent = sum(pav_current > 0.9, na.rm = TRUE) / n() * 100
  )
```

With a different threshold of 95%.
```{r}
threshold_cutoff <- 0.95

# gene body
df_pav_master %>% 
  rename(pav_current = pav_gene) %>% 
  filter(notPutativeTE == TRUE, notInComplexRegions == TRUE, covered == TRUE, !is.na(pav_current)) %>% 
  group_by(sample) %>% 
  summarise(
    count = sum(pav_current > threshold_cutoff, na.rm = TRUE), 
    total = n(), 
    percent = sum(pav_current > threshold_cutoff, na.rm = TRUE) / n() * 100
  ) 

# exon
df_pav_master %>% 
  rename(pav_current = pav_exon) %>% 
  filter(notPutativeTE == TRUE, notInComplexRegions == TRUE, covered == TRUE, !is.na(pav_current)) %>% 
  group_by(sample) %>% 
  summarise(
    count = sum(pav_current > threshold_cutoff, na.rm = TRUE), 
    total = n(), 
    percent = sum(pav_current > threshold_cutoff, na.rm = TRUE) / n() * 100
  )

# gene body (strict coverage)
df_pav_master %>% 
  rename(pav_current = pav_gene) %>% 
  filter(geneQualScore == max(geneQualScore), strictCoverage == TRUE, !is.na(pav_current)) %>% 
  group_by(sample) %>% 
  summarise(
    count = sum(pav_current > threshold_cutoff, na.rm = TRUE), 
    total = n(), 
    percent = sum(pav_current > threshold_cutoff, na.rm = TRUE) / n() * 100
  )

# exons (strict coverage)
df_pav_master %>% 
  rename(pav_current = pav_exon) %>% 
  filter(geneQualScore == max(geneQualScore), strictCoverage == TRUE, !is.na(pav_current)) %>% 
  group_by(sample) %>% 
  summarise(
    count = sum(pav_current > threshold_cutoff, na.rm = TRUE), 
    total = n(), 
    percent = sum(pav_current > threshold_cutoff, na.rm = TRUE) / n() * 100
  )
```

### PAV distribution: continuous
Shown below is the PAV score distribution for the entire gene body of selected genes across all species.
```{r pav_distribution_genebody_allspecies, fig.width=6.4, fig.height=2.2}
species_subset <- c("astCal", "mayZeb", "copChr", "rhaChi")

df_pav_master %>% 
  filter(geneQualScore == max(geneQualScore), strictCoverage == TRUE, sample %in% species_subset) %>% 
  ggplot(aes(pav_gene, fill = sample)) +
  geom_histogram(binwidth=0.05) +
  facet_wrap(.~factor(sample, levels=color_schema_species$species_abrv), nrow=1) +
  scale_fill_manual(values=color_schema_species$color, 
    breaks=color_schema_species$species_abrv, labels=color_schema_species$species_abrv) +
  xlab("Percentage conservation of gene body based on backbone annotations") +
  ylab("Number of genes") +
  theme(legend.position = 'none')
```

This suggests that overall, genes are relatively well conserved. Better conservation is observed especially for exonic regions.

### PAV distribution: discrete
The plots above have a lot of empty space in them. To make them more compact and easier to draw comparisons between, we discretise the PAV scores into uneven bins.
```{r pav_distribution_discrete, fig.width=7.7, fig.height=3}
df_pav_master_long %>% 
  filter(sample == 'aulStu', gene %in% geneList_goodCoverageStrict, !str_detect(feature, "stream")) %>% 
  ggplot(aes(cut(pav_score, breaks=c(-0.05,0.05,0.5,0.95,1.05), include.lowest = TRUE))) +
  geom_bar(stat="count", fill = '#a5c25d', color = 'black') +
  facet_wrap(.~feature, nrow = 1) +
  geom_text(aes(label = ..count..), vjust = -0.5, stat = "count") +
  scale_x_discrete(
    "Percentage conservation of feature by gene based on backbone annotations",
    breaks = c("[-0.05,0.05]", "(0.05,0.5]", "(0.5,0.95]","(0.95,1.05]", NA),
    labels = c("<5%", "5-50%", "50-95%", ">95%", "N/A")
  ) +
  scale_y_continuous("Number of genes", expand = expansion(mult = c(.05, .15)))
```

As we can see, some genes do not have introns. For interspecies comparisons, we drop those without introns.
```{r pav_distribution_discrete_allspecies, fig.width=7.2, fig.height=3.2}
genes_w_no_NA_pav <- unique(df_pav_master %>% drop_na() %>% .$gene)
sum(genes_w_no_NA_pav %in% geneList_goodCoverageStrict)

df_pav_master_long %>% 
  filter(sample != bb_name_abrv, gene %in% geneList_goodCoverageStrict, gene %in% genes_w_no_NA_pav, !str_detect(feature, "stream")) %>% 
  ggplot(aes(cut(pav_score, breaks=c(-0.05,0.05,0.5,0.95,1.05), include.lowest = TRUE),
             fill=factor(sample, levels=color_schema_species$species_abrv) )) +
  geom_bar(stat="count", position='dodge') +
  facet_wrap(.~feature, nrow = 1) +
  scale_x_discrete(
    "Percentage conservation of feature by gene based on backbone annotations",
    breaks = c("[-0.05,0.05]", "(0.05,0.5]", "(0.5,0.95]","(0.95,1.05]"),
    labels = c("<5%", "5-50%", "50-95%", ">95%")
  ) +
  ylab("Number of genes") +
  scale_fill_manual("Species", values=color_schema_species$color, 
    breaks=color_schema_species$species_abrv, labels=color_schema_species$species_abrv) +
  theme(legend.position = 'top', legend.key.size=unit(0.02, "npc"),
        legend.title = element_text(size=10), legend.text = element_text(size=9)) +
  guides(fill = guide_legend(nrow = 1))
```

```{r pav_distribution_discrete_allspecies-PAPER, fig.width=4, fig.height=3.2}
df_pav_master_long %>% 
  filter(sample != bb_name_abrv, gene %in% geneList_goodCoverageStrict, gene %in% genes_w_no_NA_pav, feature %in% c('exon','gene body')) %>% 
  ggplot(aes(cut(pav_score, breaks=c(-0.05,0.95,1.05), include.lowest = TRUE),
             fill=factor(sample, levels=rev(color_schema_species$species_abrv)) )) +
  geom_bar(stat="count", position='dodge', color='black') +
  facet_wrap(.~feature, ncol = 1) +
  scale_x_discrete(
    "A. calliptera gene conservation",
    breaks = c("[-0.05,0.95]","(0.95,1.05]"),
    labels = c("0-95%", ">95%")
  ) +
  ylab("Number of genes") +
  scale_fill_manual("Species", values=alpha(color_schema_species$color, 0.9), 
    breaks=color_schema_species$species_abrv, labels=color_schema_species$species_abrv) +
  theme(
    legend.position = 'top', 
    legend.key.size=unit(0.02, "npc"),
    legend.title = element_text(size=10), legend.text = element_text(size=9)) +
  guides(fill = guide_legend(nrow = 2)) +
  coord_flip()
```

### PAV matrix heatmap
I would like to build a heatmap-like matrix to show the PAV ratio for genes in the order they show up in the chromosome. For now, we focus on genes located on the first half of chromosome 1. The following code counts the number of genes present and retained after filtering out those in unreliable regions.
```{r}
genes_to_plot <- gr_gene[gr_gene$gene_id %in% genes_w_no_NA_pav & gr_gene$gene_id %in% geneList_goodCoverageStrict] %>% 
  subsetByOverlaps(ranges=GRanges(seqnames=chr1_id, IRanges(0, 25*10**6)))
length(genes_to_plot)

df_pav_master_long %>%  
  filter(gene %in% genes_to_plot$gene_id, sample == bb_name_abrv, feature == 'gene body') %>% 
  nrow()
```
  
PAV matrices for these genes across different genic features.
```{r pav_matrix, fig.width = 12, fig.height = 2.5}
df_pav_master_long %>%  
  filter(gene %in% genes_to_plot$gene_id, feature == 'gene body') %>% 
  ggplot(aes(x = gene, y = factor(sample, levels = color_schema_species$species_abrv), fill = pav_score*100)) +
  geom_tile(width=1, height=1, color = "black") +
  coord_fixed(ratio = 16/1) + #vertical / horizontal
  scale_fill_gradient("Proportion\nconserved", low = "white", high = "brown", limits=c(0,100)) +
  xlab("gene body") +
  scale_y_discrete("species", limits=rev) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.line.x = element_blank())

df_pav_master_long %>%  
  filter(gene %in% genes_to_plot$gene_id, feature == 'exon') %>% 
  ggplot(aes(x = gene, y = factor(sample, levels = color_schema_species$species_abrv), fill = pav_score*100)) +
  geom_tile(width=1, height=1, color = "black") +
  coord_fixed(ratio = 16/1) + #vertical / horizontal
  scale_fill_gradient("Proportion\nconserved", low = "white", high = "brown", limits=c(0,100)) +
  xlab("exon") +
  scale_y_discrete("species", limits=rev) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.line.x = element_blank())

df_pav_master_long %>%  
  filter(gene %in% genes_to_plot$gene_id, feature == 'intron') %>% 
  ggplot(aes(x = gene, y = factor(sample, levels = color_schema_species$species_abrv), fill = pav_score*100)) +
  geom_tile(width=1, height=1, color = "black") +
  coord_fixed(ratio = 16/1) + #vertical / horizontal
  scale_fill_gradient("Proportion\nconserved", low = "white", high = "brown", limits=c(0,100)) +
  xlab("intron") +
  scale_y_discrete("species", limits=rev) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.line.x = element_blank())

df_pav_master_long %>%  
  filter(gene %in% genes_to_plot$gene_id, feature == 'upstream') %>% 
  ggplot(aes(x = gene, y = factor(sample, levels = color_schema_species$species_abrv), fill = pav_score*100)) +
  geom_tile(width=1, height=1, color = "black") +
  coord_fixed(ratio = 16/1) + #vertical / horizontal
  scale_fill_gradient("Proportion\nconserved", low = "white", high = "brown", limits=c(0,100)) +
  xlab("upstream") +
  scale_y_discrete("species", limits=rev) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.line.x = element_blank())

df_pav_master_long %>%  
  filter(gene %in% genes_to_plot$gene_id, feature == 'downstream') %>% 
  ggplot(aes(x = gene, y = factor(sample, levels = color_schema_species$species_abrv), fill = pav_score*100)) +
  geom_tile(width=1, height=1, color = "black") +
  coord_fixed(ratio = 16/1) + #vertical / horizontal
  scale_fill_gradient("Proportion\nconserved", low = "white", high = "brown", limits=c(0,100)) +
  xlab("downstream") +
  scale_y_discrete("species", limits=rev) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.line.x = element_blank())
```

## Gene lists of interest
In this category, we do need to filter for strict, because the idea is to find private genes to the backbone. Considering we have a fairly generous criteria that at least 5 samples must have coverage, what we find as absent in the others is likely to be true.

### Dealing with markedly absent genes
Based on the gene body, there seem to be genes that are completely absent in many species. 
```{r prevalence_of_gene, fig.width=7.8, fig.height=3.2}
df_pav_master_long_summarised <- df_pav_master_long %>% 
  mutate(
    pav_score_is_very_high = pav_score > 0.9,
    pav_score_is_high = pav_score > 0.5 # less strict to detect truly missing genes
  ) %>% 
  filter(sample != bb_name_abrv) %>%
  group_by(gene, feature) %>% 
  summarise(
    high_pav_count = sum(pav_score_is_high) + 1,
    very_high_pav_count = sum(pav_score_is_very_high) + 1
  )

df_pav_master_long_summarised %>% 
  filter(!str_detect(feature, "stream")) %>% 
  ggplot(aes(very_high_pav_count)) +
  geom_bar(stat='count') +
  geom_text(aes(label = ..count..), vjust = -0.5, stat = "count", size=2.8) +
  scale_x_continuous("Number of samples with >90% gene sequence conserved", breaks = seq(1,8)) +
  scale_y_continuous("Number of genes", expand = expansion(mult = c(.05, .1)) ) +
  facet_wrap(.~feature, nrow = 1)
```

```{r prevalence_of_gene_bodyonly, fig.width=3, fig.height=2.7}
df_pav_master_long_summarised %>% 
  filter(feature == "gene body") %>% 
  ggplot(aes(high_pav_count)) +
  geom_bar(stat='count', fill = '#D65F5F', color = 'black') +
  geom_text(aes(label = ..count..), vjust = -0.5, stat = "count", size=2.8) +
  scale_x_continuous("Number of samples with\n>50% gene sequence conserved", breaks = seq(1,8)) +
  scale_y_continuous("Number of genes", expand = expansion(mult = c(.05, .1)) ) 

df_pav_master_long_summarised %>% 
  filter(feature == "gene body") %>% 
  ggplot(aes(very_high_pav_count)) +
  geom_bar(stat='count', fill = '#D65F5F', color = 'black') +
  geom_text(aes(label = ..count..), vjust = -0.5, stat = "count", size=2.8) +
  scale_x_continuous("Number of samples with\n>90% gene sequence conserved", breaks = seq(1,8)) +
  scale_y_continuous("Number of genes", expand = expansion(mult = c(.05, .1)) ) 
```

Print some stats. Count the number of genes present in all the species at different thresholds of sequence conservation.
```{r}
df_pav_master_long_summarised %>% 
  filter(feature %in% c("gene body", "exon")) %>% 
  group_by(feature, very_high_pav_count) %>% 
  filter(very_high_pav_count == 8) %>% 
  summarise(
    count = n(), 
    percent = n()/highQualGeneCount * 100
  ) %>% 
  as.data.frame()

df_pav_master_long_summarised %>% 
  filter(feature %in% c("gene body", "exon")) %>% 
  group_by(feature, high_pav_count) %>% 
  filter(high_pav_count == 8) %>% 
  summarise(
    count = n(), 
    percent = n()/highQualGeneCount * 100
  ) %>% 
  as.data.frame()
```

Write list of genes present in seven or eight samples.
```{r}
df_pav_master_long_summarised %>% 
  filter(feature == "gene body", high_pav_count >= 8) %>% 
  .$gene %>% writeLines(glue("./local/gene_ontology/geneList_Universal-{bb_name_abrv}.txt"))
```

Here we obtain the genes that are markedly absent from all the non-backbone samples. We also write to a CSV file their gene descriptions to better see what they might be.
```{r}
# load gene descriptions
gene_descriptions <- readr::read_csv(
  glue(Sys.getenv("HOME"), "/code/malawi_transposon/cloud/genome/{bb_name}/gene_description.csv"),
  skip = 1, col_names = c("gene", "gene_description", "gene_name") ) %>% 
  mutate(gene_description = gsub("\\[Source:.*?\\]", "", gene_description))

# get list of backbone genes completely absent from others
list_of_absent_genes <- df_pav_master_long_summarised %>% 
  filter(feature == "gene body", high_pav_count == 1) %>% 
  left_join(select(gene_descriptions, gene, gene_description))

# write to CSV file to be used for functional enrichment and word cloud
list_of_absent_genes %>% 
  select(gene, gene_description) %>% 
  write_csv(file = glue("./local/gene_ontology/geneList_markedlyAbsent-{bb_name_abrv}.csv") )

glue(
  "There are {nrow(list_of_absent_genes)} backbone genes not detected in other species. ",
  "{sum(is.na(list_of_absent_genes$gene_description))} do not have a description."
)
```

Stricter set where must have coverage.
```{r}
list_of_absent_genes_Strict <- list_of_absent_genes %>% filter(gene %in% geneList_goodCoverageStrict)

list_of_absent_genes_Strict %>% 
  filter(gene %in% geneList_goodCoverageStrict) %>% 
  select(gene, gene_description) %>% 
  write_csv(file = glue("./local/gene_ontology/geneList_markedlyAbsent-{bb_name_abrv}-strict.csv") )

glue(
  "There are {nrow(list_of_absent_genes_Strict)} backbone genes not detected in other species (strict). ",
  "{sum(is.na(list_of_absent_genes_Strict$gene_description))} do not have a description."
)
```


## Final notes
There are some further issues which are worth bearing in mind:

* I do not know how trustworthy the PAV conservation scores are for upstream and downstream regions, because they tend to overflow into another gene.
* It is worth doing a quick ortholog analysis for the genes that are markedly absent and have poor coverage between different backbone choices.

## Specific to calliptera
Here's a list of genes that are completely "*absent*" in the non-backbone species on the first half of chromosome 1.
```{r}
df_pav_master_long_summarised %>% 
  filter(feature == "gene body", high_pav_count == 1) %>% 
  filter(gene %in% genes_to_plot$gene_id)
```

From this, it is tempting to say that there are genes that are completely absent in certain species. Out of these three:

* ENSACLG00000007160 is a 65S ribosomal gene with a proper gene annotation. Probably correct!
* ENSACLG00000010070 actually shows still quite high PAV score of about 0.4 for the others. However, it has a large complex bubble within it. Might be worth adding another filter to say that the gene cannot overlap 50% of its gene body with a single complex bubble (just like for the TE).
* ENSACLG00000014925 is a novel gene. Could be a transposon? Should use BioMart annotations to filter out.

Simple boilerplate code to obtain gene coordinates for IGV inspection.
```{r}

gene_to_test <- "ENSACLG00000010070"
df_pav_master_long %>% filter(gene == gene_to_test, feature == "gene body")
gr_gene[gr_gene$gene_id == gene_to_test]

df_sv %>% filter(chr=='1', start> 24370000)
```

## Custom plots
```{r prevalence_of_gene_bodyonly-PAPER, fig.width=3.8, fig.height=2.8}
a1 <- c(378, 79, 92, 86, 77, 199, 2109, 21302)
a2 <- c(260, 78, 66, 63, 69, 242, 2396, 20964)

tibble(
  sample = c(rep('astCal', 8), rep('mayZeb', 8)), 
  gene_count = c(a1, a2),
  sample_freq = c(seq(1,8), seq(1,8)) ) %>% 
  ggplot(aes(x = sample_freq, y = gene_count, fill = sample)) +
  geom_bar(stat='identity', color = 'black', position = 'dodge') +
  geom_text(aes(label = gene_count, color = sample), 
            angle = 90, vjust = 0.5, hjust = -0.2, size = 2.8, fontface = 'bold',
            position = position_dodge(width = 1)
  ) +
  scale_x_continuous("Samples with over 50% gene conservation", breaks = seq(1,8)) +
  scale_y_continuous("Number of genes", expand = expansion(mult = c(.05, .25)) ) +
  scale_fill_manual("Backbone", values = c('#D65F5F', '#82C6E2')) +
  scale_color_manual("Backbone", values = c('#D65F5F', '#82C6E2')) +
  theme(
    legend.position = c(0.2, 0.9), legend.key.size=unit(0.05, "npc"),
    legend.title = element_text(size=10), legend.text = element_text(size=9)
  )
```