---
title: "Linking SVs to annotated genes"
output: html_document
---
  
```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE, message = FALSE, warning = FALSE,
  cache = TRUE, cache.path = '.cache_knitr/link_sv_to_genes/',
  fig.width = 7, fig.height = 4, dpi = 350,
  fig.path='figure/2b-link_sv_to_genes/',
  dev = c('png', 'pdf')
)
```

In this document, we associate SVs to nearby genes. I have some reservations about whether I should shrink the SVs to only their midpoints for this comparison. I feel like for most analyses, it is okay to use them as they are, and I tried running both versions to see if the results are any different. Most are the same, except for:

* Calculating SV density in upstream and downstream regions, due to their short length.
* Gene Ontology for flanking regions gives different results for **flanking regions**. Midpoint gives nothing basically.
* I notice that the exon intersections are fewer. 1573 vs 3754 (density: 0.038 vs 0.090).
* Distance from relative start and end (not doing the midpoint will cause a skew to upstream).

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

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

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

updown_distance <- 2000  # distance considered upstream / downstream of a gene

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

## Import pangenome bubble coordinates
We import the data frame for SVs that has been preprocessed in `1-preprocess.Rmd`. We then use some of the relevant columns to create a GRanges object.
```{r}
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", ".")
)

chrom_sizes <- read_delim(
  glue(Sys.getenv("HOME"), "/igv/genomes/{bb_name}/chrom.sizes"),
  col_names = c('chr', 'len')
)
chrom_sizes <- Seqinfo(chrom_sizes$chr, chrom_sizes$len)

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

# resize the coordinates so it is only the midpoint of bubbles
#gr_sv <- resize(gr_sv_raw, width=1, fix="center")
gr_sv <- gr_sv_raw
```

## Obtain gene annotations
Load annotation for the `Astatotilapia calliptera` genome.
```{r}
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)
```

We also get a list of genes that might be fake putative TEs (see `4-pav_plots.Rmd` for more info). 
```{r}
df_gene_trueness <- readr::read_csv(
  glue(Sys.getenv("HOME"), "/code/malawi_transposon/cloud/genome/{bb_name}/gene_trueness.csv")
)

glue("{length(x)} genes overlap significantly with TE",
  x = df_gene_trueness %>%
    filter((maxOverlapSingleTE > threshGeneBodyOverlapSingleTE | totalOverlapMultipleTE > threshGeneBodyOverlapMultipleTE) ) %>% .$gene_id
)

glue("{length(x)} genes are associated with the DNA transposition term",
  x = df_gene_trueness %>%
    filter(hasTranspositionTerm == TRUE) %>% .$gene_id
)

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)
```


## Create GRanges for gene features
From this, we obtain GRanges for various type of genomic features, filtering out the false genes.
```{r}
gr_gene <- genes(txdb, columns = c('gene_id', 'tx_name'))
gr_gene <- gr_gene[!gr_gene$gene_id %in% list_of_dodgy_genes]

# exons in all coding and non-coding genes
gr_exon <- exons(txdb, columns = c('exon_name', "GENEID"))
gr_exon$GENEID <- sapply(gr_exon$GENEID, unlist)
gr_exon <- gr_exon[!gr_exon$GENEID %in% list_of_dodgy_genes]

# CDS: translated regions of the exons in protein-coding genes
gr_cds <- cds(txdb, columns = c('cds_name', "GENEID"))
gr_cds$GENEID <- sapply(gr_cds$GENEID, unlist)
gr_cds <- gr_cds[!gr_cds$GENEID %in% list_of_dodgy_genes]

glue("number of genes: {length(gr_gene)}")
glue("number of PROTEIN-CODING genes: {unlist(gr_cds$GENEID) %>% unique() %>% length()}")

# introns
gr_intron <- GenomicRanges::setdiff(gr_gene, gr_exon)

# 5-UTRs
gr_5utr <- unlist(fiveUTRsByTranscript(txdb, use.names=TRUE))
gr_5utr <- gr_5utr[gr_5utr$exon_name %in% gr_exon$exon_name]

# 3-UTRs
gr_3utr <- unlist(threeUTRsByTranscript(txdb, use.names=TRUE))
gr_3utr <- gr_3utr[gr_3utr$exon_name %in% gr_exon$exon_name]

# modified version of exon GRanges
# contains the CDS of coding genes and retains the exons of non coding genes
gr_5and3_utr <- GenomicRanges::union(gr_5utr, gr_3utr)
gr_exon_2 <- GenomicRanges::setdiff(gr_exon, gr_5and3_utr, ignore.strand=TRUE)

# upstream and downsream of genes
gr_upstream <- GenomicRanges::trim(GenomicRanges::flank(gr_gene, start=T, width=updown_distance))
gr_downstream <- GenomicRanges::trim(GenomicRanges::flank(gr_gene, start=F, width=updown_distance))

# intergenic regions (where strand information is dropped)
gr_gene_unstranded <- gr_gene
strand(gr_gene_unstranded) <- '*'
gr_intergenic <- gaps(gr_gene_unstranded)[strand(gaps(gr_gene_unstranded))=='*',]
```

Note: an alternative method to retrieve intronic regions is the `intronsByTranscripts` function. I chose not to use it, because they are parts of gene which are (for some reason not annotated as an exon, intron or anything: in particular `1:13917706-13919380, 1:670912-690561, 1:710725-713473`). I would want those to be classified as introns, thus I used my method above.

## Find overlap of SVs with genic features
To resolve overlapping features for a particular region, we use a hierachical logic where hits with certain features are reported with higher priority. This means that an intersection of an SV with a particular genomic region will only be reported once, even if that region has multiple associated gene annotations. Currently, an intersection is counted even if it's just 1bp overlap. This might need to be changed.

This `for` loop also calculates the total length of each feature.
```{r}
gr_sv$feature <- ""
gr_explored <- GRanges(seqinfo = chrom_sizes)

gr_list_ordered <- list(gr_exon_2, gr_5utr, gr_3utr, gr_intron, gr_upstream, gr_downstream, gr_intergenic)
feature_names_ordered <- c("exon", "5UTR", "3UTR", "intron", "upstream", "downstream", "intergenic")

feature_lengths_ordered <- rep(0, length(feature_names_ordered))
names(feature_lengths_ordered) <- feature_names_ordered

for (i in seq(1, length(gr_list_ordered)) ) {
  gr_to_add <- gr_list_ordered[[i]]
  feature_name <- feature_names_ordered[i]
    
  gr_subject <- GenomicRanges::setdiff(gr_to_add, gr_explored, ignore.strand=TRUE)
  feature_lengths_ordered[i] <- sum(lengths(gr_subject))
  
  idx <- findOverlaps(gr_sv, gr_subject)
  gr_sv$feature[idx@from][(gr_sv$feature=="")[idx@from]] <- feature_name
  gr_explored <- GenomicRanges::union(gr_explored, gr_subject)
}

head(as.data.frame(gr_sv))

feature_lengths_ordered
```

## Obtaining nearest feature
Distance to nearest gene for each SV.
```{r}
idx_gene <- distanceToNearest(gr_sv, gr_gene, ignore.strand=TRUE)

# get gene id
gr_sv$gene <- NA
gr_sv$gene[idx_gene@from] <- gr_gene$gene_id[idx_gene@to]

# get distance
gr_sv$geneDistance <- NA
gr_sv$geneDistance[idx_gene@from] <- mcols(idx_gene)$distance
```

Retrieve the distance to the nearest exon (we don't care as much about which actual exon it is).
```{r}
idx_exon <- distanceToNearest(gr_sv, gr_exon, ignore.strand=TRUE)
gr_sv$exonDistance <- NA
gr_sv$exonDistance[idx_exon@from] <- mcols(idx_exon)$distance
```

### To nearest gene start and gene end
First, we create GRanges to represent the start and end of genes. The width of the objects in the GRanges are set to 0.
```{r}
df_gene <- as.data.frame(gr_gene) %>% 
  as_tibble() %>% 
  dplyr::select(gene_id, seqnames, start, end, strand) %>% 
  mutate(
    true_start = if_else(strand == '+', start, end),
    true_end = if_else(strand == '+', end, start)
  )

gr_start <- GRanges(
  seqnames = df_gene$seqnames, ranges = IRanges(start = df_gene$true_start, width=0),
  strand = df_gene$strand, gene_id = df_gene$gene_id
)

gr_end <- GRanges(
  seqnames = df_gene$seqnames, ranges = IRanges(start = df_gene$true_end, width=0),
  strand = df_gene$strand, gene_id = df_gene$gene_id
)
```

For this, I would like to compute directional distance, where negative indicates that the SV is upstream of the gene start, while positive means that it is downstream. This is done by means of a custom made function called `compute_shortest_directional_distance`, because the `distancetoNearest` function from GenomicFeatures does not achieve this.

What this function does is to check whether the SV intersects with the gene start/end, in which case the distance is set to 0. Otherwise, it compares the upstream and downstream distances, and choose the one with the smallest absolute value. 
```{r}
source("modules/compute_shortest_directional_distance.R")
compute_shortest_directional_distance(gr_sv[1:10], gr_start)  # example
```

This is computed for both `gr_start` and `gr_end`, and merged with the main SV dataframe.
```{r}
tmp_df <- as_tibble(mcols(gr_sv))

tmp_df <- left_join(tmp_df, 
  compute_shortest_directional_distance(gr_sv, gr_start) %>% 
    dplyr::rename(startDistance = distance, startGene = gene), by='id')

tmp_df <- left_join(tmp_df,
  compute_shortest_directional_distance(gr_sv, gr_end) %>% 
    dplyr::rename(endDistance = distance, endGene = gene), by='id')

df_sv2 <- left_join(df_sv, tmp_df, by='id') %>% 
  relocate(feature:endGene, .after = lenbio_minmaxratio_categorised) %>% 
  mutate(feature = forcats::fct_relevel(feature, feature_names_ordered)
)
```

## Plots
Genomic feature where SVs are overlapping with.
```{r feature_overlap, fig.width=3.5, fig.height=3}
df_sv2 %>% 
  ggplot(aes(feature, fill=simple)) + 
  geom_bar(stat="count") +
  xlab(NULL) +
  ylab("Number of bubbles") +
  labs(fill="simple bubble") +
  theme(legend.position = 'top')

df_sv2 %>% 
  ggplot(aes(feature, fill = feature)) + 
  geom_bar(stat="count", color='black') +
  scale_y_continuous(expand = expansion(mult = c(.05, .1))) +
  geom_text(aes(label = ..count..), vjust = -0.5, stat = "count") +
  scale_fill_brewer('Genomic feature', palette = 'Pastel1') +
  xlab(NULL) +
  ylab("Number of bubbles") +
  theme(legend.position = 'none', axis.text.x = element_text(angle = 20, hjust = 1) )
```

Version where we normalise by feature length to calculate SV density. **WORTH COLOR AND GROUPING BY SIZE OF BUBBLE**.
```{r feature_overlap_density, fig.width=3.5, fig.height=3}
df_sv2 %>%
  group_by(feature) %>% 
  summarise(count = n()) %>% 
  left_join(tibble(feature=names(feature_lengths_ordered), feature_length=feature_lengths_ordered), by='feature') %>% 
  mutate(feature = forcats::fct_relevel(feature, feature_names_ordered)) %>% 
  ggplot(aes(x = feature, y = count / feature_length * 1000, fill = feature)) + 
  geom_bar(stat="identity", color='black') +
  scale_y_continuous(expand = expansion(mult = c(.05, .1))) +
  geom_text(aes(label = format(round(count / feature_length * 1000, 4), n_small=4) ), size=3.3, vjust = -0.5) +
  scale_fill_brewer('Genomic feature', palette = 'Pastel2') +
  xlab(NULL) +
  ylab("Bubble density per kbp of feature") +
  theme(legend.position = 'none', axis.text.x = element_text(angle = 20, hjust = 1) )
```

```{r feature_overlap_density-poster, fig.width=4.2, fig.height=2.4}
df_sv2 %>%
  group_by(feature) %>% 
  summarise(count = n()) %>% 
  left_join(tibble(feature=names(feature_lengths_ordered), feature_length=feature_lengths_ordered), by='feature') %>% 
  mutate(feature = forcats::fct_relevel(feature, feature_names_ordered)) %>% 
  ggplot(aes(x = feature, y = count / feature_length * 1000, fill = feature)) + 
  geom_bar(stat="identity", color='black') +
  scale_y_continuous(expand = expansion(mult = c(.05, .15))) +
  geom_text(aes(label = format(round(count / feature_length * 1000, 3), n_small=3) ), vjust = -0.5) +
  scale_fill_brewer('Genomic feature', palette = 'Pastel2') +
  xlab(NULL) +
  ylab("Bubble density per kbp") +
  theme_linedraw() +
  theme(
    legend.position = 'none', 
    axis.text.x = element_text(angle = 35, hjust = 1)
  )
```

Distribution of simple bubbles by feature.
```{r feature_overlap_facet}
df_sv2 %>%
  ggplot(aes(simple, fill=simple)) + 
  geom_bar(stat="count") +
  facet_wrap(.~feature, scales = 'free_y', nrow=2) +
  scale_x_discrete("Simple bubble", limits=rev) +
  ylab("Number of bubbles") +
  theme(legend.position = 'none')
```

Distance to nearest gene.
```{r distance_nearest_gene, fig.width=6, fig.height=2.8}
df_sv2 %>%
  ggplot(aes(geneDistance + 1, fill=feature)) + 
  geom_histogram(bins=50, color='black') +
  scale_fill_brewer('Genomic feature', palette = 'Pastel2') +
  scale_x_log10("Distance to nearest gene (bp)",
                breaks = trans_breaks("log10", function(x) 10^x), 
                labels = trans_format("log10", math_format(10^.x)) ) +
  ylab("Number of bubbles") +
  theme(
    legend.title = element_text(size = rel(0.8)),
    legend.text = element_text(size = rel(0.8)),
    legend.key.size=unit(0.025, "npc")
  )
```

Print some statistics
```{r}
glue(
  "SVs located in genes: {sum(df_sv2$geneDistance == 0, na.rm=TRUE)}\n",
  "SVs located in genes and within 2000bp of start/end: {sum(df_sv2$geneDistance < 2000, na.rm=TRUE)}\n",
  "Number of SVs in total: {nrow(df_sv)}"
)
```

Distance to nearest exon.
```{r distance_nearest_exon, fig.width=6, fig.height=2.8}
df_sv2 %>%
  ggplot(aes(exonDistance + 1, fill=feature)) + 
  stat_bin(breaks = sapply(seq(0,6,0.1), function(x) 10**x), color='black') +
  scale_fill_brewer('Genomic feature', palette = 'Pastel2') +
  scale_x_log10("Distance to nearest exon (bp)",
    breaks = trans_breaks("log10", function(x) 10^x), 
    labels = trans_format("log10", math_format(10^.x)) ) +
  ylab("Number of bubbles") +
  theme(
    legend.title = element_text(size = rel(0.8)),
    legend.text = element_text(size = rel(0.8)),
    legend.key.size=unit(0.025, "npc")
  )
```

Distance to nearest gene start. Note that we need to filter for the cases where the SV is in the gene and the start of the gene also happens to be the closest to the SV. Otherwise, the results would be misleading. 
```{r distance_nearest_geneStart, fig.width=4, fig.height=2.6}
df_sv2 %>% 
  filter(gene == startGene) %>%
  drop_na(startDistance, startGene) %>% 
  ggplot(aes(startDistance/1000)) +
  geom_histogram(breaks=seq(-500, 500, 25), color='black', fill='gray80') +
  geom_vline(xintercept = 0, color='black', linetype='dashed') +
  scale_x_continuous('Distance to nearest gene start (kbp)', labels = comma) +
  ylab("Number of bubbles") +
  theme(aspect.ratio=1)

# same as above, but colored based on what feature the SV is in
df_sv2 %>% 
  filter(gene == startGene) %>%
  drop_na(startDistance, startGene) %>% 
  ggplot(aes(startDistance/1000, fill=feature)) +
  geom_histogram(breaks=seq(-500, 500, 25), color='black') +
  geom_vline(xintercept = 0, color='black', linetype='dashed') +
  scale_fill_brewer('Genomic feature', palette = 'Pastel2') +
  scale_x_continuous('Distance to nearest gene start (kbp)', labels = comma) +
  ylab("Number of bubbles") +
  theme(aspect.ratio=1)
```

Distance to nearest gene end. As above, we need to filter for the cases where the SV is in the gene and the end of the gene also happens to be the closest to the SV. 
```{r distance_nearest_geneEnd, fig.width=4, fig.height=2.6}
df_sv2 %>% 
  filter(gene == endGene) %>% 
  drop_na(endDistance, endGene) %>% 
  ggplot(aes(endDistance/1000)) +
  geom_histogram(breaks=seq(-500, 500, 25), color='black', fill='gray80') +
  geom_vline(xintercept = 0, color='black', linetype='dashed') +
  scale_x_continuous('Distance to nearest gene end (kbp)', labels = comma) +
  ylab("Number of bubbles") +
  theme(aspect.ratio=1)

df_sv2 %>% 
  filter(gene == endGene) %>% 
  drop_na(endDistance, endGene) %>% 
  ggplot(aes(endDistance/1000)) +
  geom_histogram(aes(fill=feature), breaks=seq(-500, 500, 25), color='black') +
  geom_vline(xintercept = 0, color='black', linetype='dashed') +
  scale_fill_brewer('Genomic feature', palette = 'Pastel2') +
  scale_x_continuous('Distance to nearest gene end (kbp)', labels = comma) +
  ylab("Number of bubbles") +
  theme(aspect.ratio=1)
```

Here are some useful tables with figures we quote in the thesis:
```{r}
# counts of bubbles and density by feature
df_sv2 %>%
  group_by(feature) %>% 
  summarise(count = n()) %>% 
  left_join(tibble(feature=names(feature_lengths_ordered), feature_length=feature_lengths_ordered), by='feature') %>% 
  mutate(density = round(count / feature_length * 1000, 4)) %>% 
  reactable()

# distance from nearest exon (all)
mean(df_sv2$exonDistance, na.rm=TRUE)

# distance from nearest exon (intergenic)
mean(df_sv2 %>% filter(!feature %in% c("exon", "5UTR", "3UTR", "intron")) %>%.$exonDistance, na.rm=TRUE)

# distance from nearest gene start (intergenic)
mean(df_sv2 %>% filter(!feature %in% c("exon", "5UTR", "3UTR", "intron")) %>% .$startDistance %>% abs(), na.rm=TRUE)

# distance from nearest gene end (intergenic)
mean(df_sv2 %>% filter(!feature %in% c("exon", "5UTR", "3UTR", "intron")) %>% .$endDistance %>% abs(), na.rm=TRUE)
```


## More gene-centric analysis
The analyses above are one-to-one: each SV is associated to its closest gene. Here we take an all-in-all approach where we consider every possible hit. This is useful for some hypotheses: which genes have multiple SVs in them, and vice versa.

Find overlaps of SV with upstream, downstream regions followed by gene body. For each hit, we also compute the distance of the SV from the start of that gene.
```{r}
# SV located before a gene start, up to -2000bp away
hits <- findOverlaps(gr_upstream, gr_sv)
tmp_df1 <- gr_sv[hits@to, "id"] %>% as_tibble() %>% 
  add_column(
    geneID = gr_start[hits@from]$gene_id,
    distanceToGeneStart = -distance(gr_start[hits@from], gr_sv[hits@to]),
    distanceToGeneEnd = -distance(gr_end[hits@from], gr_sv[hits@to])) %>% 
  left_join(
    tibble(geneID = gr_gene$gene_id, geneLength = lengths(gr_gene)), by = "geneID") %>% 
  filter(abs(distanceToGeneStart) <= updown_distance) %>% 
  mutate(
    relDistanceToGeneStart = distanceToGeneStart / geneLength,
    relDistanceToGeneEnd = distanceToGeneEnd / geneLength
  )

# SV located after a gene end, up to 2000bp away
hits <- findOverlaps(gr_downstream, gr_sv)
tmp_df2 <- gr_sv[hits@to, "id"] %>% as_tibble() %>% 
  add_column(
    geneID = gr_start[hits@from]$gene_id,
    distanceToGeneStart = distance(gr_start[hits@from], gr_sv[hits@to]),
    distanceToGeneEnd = distance(gr_end[hits@from], gr_sv[hits@to])) %>% 
  left_join(
    tibble(geneID = gr_gene$gene_id, geneLength = lengths(gr_gene)), by = "geneID") %>% 
  filter(distanceToGeneStart < geneLength + updown_distance) %>% 
  mutate(
    relDistanceToGeneStart = distanceToGeneStart / geneLength,
    relDistanceToGeneEnd = distanceToGeneEnd / geneLength
  )

# SV located in the gene body
hits <- findOverlaps(gr_gene, gr_sv)
tmp_df3 <- gr_sv[hits@to, "id"] %>% as_tibble() %>% 
  add_column(
    geneID = gr_start[hits@from]$gene_id,
    distanceToGeneStart = distance(gr_start[hits@from], gr_sv[hits@to]),
    distanceToGeneEnd = -distance(gr_end[hits@from], gr_sv[hits@to])) %>% 
  left_join(
    tibble(geneID = gr_gene$gene_id, geneLength = lengths(gr_gene)), by = "geneID") %>% 
  filter(distanceToGeneStart < geneLength + updown_distance) %>% 
  mutate(
    relDistanceToGeneStart = distanceToGeneStart / geneLength,
    relDistanceToGeneEnd = distanceToGeneEnd / geneLength
  )

# combine the dataframes
df_svHitsByGene <- bind_rows(tmp_df1, tmp_df2, tmp_df3) %>% 
  rename(chr=seqnames) %>% 
  arrange(chr, start, geneID) %>% 
  distinct()  # need this because of complete overlap
rm(tmp_df1, tmp_df2, tmp_df3, hits)

df_svHitsByGene
```

Note that in the code above when the temporary dataframes were combined, there is a `distinct()` command which is needed to remove multiple hits of a very long bubble that completely covers a gene. Such situations will give three hits to the upstream, downstream and gene body, but the distanceToGeneStart will all be 0. This would be misleading. An example of such an SV is `3:28234379-28319387`.

Below we plot the relative distance of SV hits to gene. SVs do not seem to be biased to any side of the gene.
```{r all2all_relDistanceToGeneStart, fig.width=3, fig.height=2.5}
# relative distance of SV compared to gene
# there is not more upstream, artifact of analysis since we prioritise start of gene
df_svHitsByGene %>% 
  ggplot(aes(relDistanceToGeneStart)) +
  geom_histogram(fill='gray80', color='black', binwidth=0.1) +
  coord_cartesian(xlim=c(-1, 2)) +
  xlab("Relative distance from gene start") +
  ylab("Number of bubbles")
```

Aggregate counts of SV hits by gene, and compute density of SVs per 1kbp. 
```{r}
# tally up all SVs by gene
# this is useful to find genes with SVs enriched in them
# the next step would be to filter for meaningful SVs
genes_w_noSV <- names(gr_gene)[!names(gr_gene) %in% unique(df_svHitsByGene$geneID)] # genes with no SV
df_svHitsByGene_agg <- df_svHitsByGene %>%
  group_by(geneID) %>% 
  summarise(
    n_sv_all = n(), 
    n_sv_up = sum(relDistanceToGeneStart<=0),
    n_sv_body = sum(relDistanceToGeneStart>0 & relDistanceToGeneStart<1),
    n_sv_down = sum(relDistanceToGeneStart>=1) ) %>% 
  bind_rows(tibble(geneID = genes_w_noSV, n_sv_all = 0, n_sv_up = 0, n_sv_body = 0, n_sv_down = 0) ) %>%  # genes with no SVs
  arrange(geneID) %>% 
  left_join(tibble(geneID = names(gr_gene), geneLength = width(gr_gene)), by = 'geneID') %>% 
  mutate(
    density_all = n_sv_all / (geneLength + updown_distance) * 1000,
    density_body = n_sv_body / geneLength * 1000,
    density_up = n_sv_up / updown_distance * 1000,
    density_down = n_sv_down / updown_distance * 1000,
    density_flank = (n_sv_up + n_sv_down) / (2*updown_distance) * 1000
  )

glue("average SV density in genic regions\n",
  "overall    : {mean(df_svHitsByGene_agg$density_all)} ± {sd(df_svHitsByGene_agg$density_all) / sqrt(nrow(df_svHitsByGene_agg))} \n",
  "gene body  : {mean(df_svHitsByGene_agg$density_body)} ± {sd(df_svHitsByGene_agg$density_body) / sqrt(nrow(df_svHitsByGene_agg))} \n",
  "upstream   : {mean(df_svHitsByGene_agg$density_up)} ± {sd(df_svHitsByGene_agg$density_up) / sqrt(nrow(df_svHitsByGene_agg))} \n",
  "downstream : {mean(df_svHitsByGene_agg$density_down)} ± {sd(df_svHitsByGene_agg$density_down) / sqrt(nrow(df_svHitsByGene_agg))} \n",
  "flank      : {mean(df_svHitsByGene_agg$density_flank)} ± {sd(df_svHitsByGene_agg$density_flank) / sqrt(nrow(df_svHitsByGene_agg))} \n\n"
)

glue("SV density in genic region (arithmetic sum)\n",
  "overall    : { sum(df_svHitsByGene_agg$n_sv_all) / sum(df_svHitsByGene_agg$geneLength) * 1000 }\n",
  "gene body  : { sum(df_svHitsByGene_agg$n_sv_body) / sum(df_svHitsByGene_agg$geneLength) * 1000 }\n",
)
```

SV density distribution across gene body.
```{r}
summary(df_svHitsByGene_agg$density_body)

# summary without those that are zero
df_svHitsByGene_agg %>% filter(density_body!=0) %>% .$density_body %>% summary()

# those that have very high SV bubble densities (most because their genes are short)
sum(df_svHitsByGene_agg$density_body > 1)
```

Whole distribution
```{r all2all_svdensity_entire_unzoom, fig.width=7, fig.height=2}
df_svHitsByGene_agg %>% 
  ggplot(aes(density_body)) +
  geom_histogram(fill='black', binwidth=0.1) + 
  xlab("SV density per 1kbp (gene body)") +
  ylab(NULL) 
```

Zoomed version.
```{r all2all_svdensity_entire, fig.width=3.2, fig.height=2.5}
median_nonzero_density <- df_svHitsByGene_agg %>% 
  filter(density_body!=0) %>% .$density_body %>% median() %>% 
  round(3)

df_svHitsByGene_agg %>% 
  ggplot(aes(density_body)) +
  geom_vline(xintercept = median_nonzero_density, linetype=2, color='darkred') +
  geom_histogram(fill='gray80', color='black', binwidth=0.02) + 
  annotate(geom='text', x=median_nonzero_density+0.03, y=Inf, 
           label=glue("non-zero median = {round(median_nonzero_density, 3)}"), 
           size = rel(2.6), hjust=1, angle=90) +
  coord_cartesian(xlim=c(0,1)) +
  xlab("SV density per 1kbp (gene body)") +
  ylab("Number of genes") 
```

Hit density upstream and downstream of genes.
```{r all2all_svdensity_updown, fig.width=3.2, fig.height=2.5}
# count SV upstream of genes
df_svHitsByGene_agg %>% 
  ggplot(aes(x = cut(n_sv_up, breaks=c(-1,0,1,2,100), include.lowest=FALSE) )) +
  geom_bar() +
  geom_text(aes(label = scales::percent((..count..)/sum(..count..)) ), vjust = -0.5, stat = "count") +
  scale_x_discrete(glue("number of SVs {updown_distance}bp upstream"), labels=c("0", "1", "2", "more than 2")) +
  scale_y_continuous("Number of genes", expand = expansion(mult = c(.05, .1)) ) 

# count SV downstream of genes
df_svHitsByGene_agg %>% 
  ggplot(aes(x = cut(n_sv_down, breaks=c(-1,0,1,2,100), include.lowest=FALSE) )) +
  geom_bar() +
  geom_text(aes(label = scales::percent((..count..)/sum(..count..)) ), vjust = -0.5, stat = "count") +
  scale_x_discrete(glue("number of SVs {updown_distance}bp downstream"), labels=c("0", "1", "2", "more than 2")) +
  scale_y_continuous("Number of genes", expand = expansion(mult = c(.05, .1)) )
```

Print some statistics:
```{r}
glue(
  "ALL BUBBLES\n",
  "Genes with SVs            : {all}\n",
  "Genes with SVs gene body  : {body}\n",
  "Genes with SVs upstream   : {up}\n",
  "Genes with SVs downstream : {down}\n",
  "Genes with SVs in flanks  : {flank}\n",
  tmp_df = df_svHitsByGene_agg,
  all = sum(tmp_df$n_sv_all > 0),
  body = sum(tmp_df$n_sv_body > 0),
  up = sum(tmp_df$n_sv_up > 0),
  down = sum(tmp_df$n_sv_down > 0),
  flank = sum((tmp_df$n_sv_up + tmp_df$n_sv_down) > 0)
)
```

## Making gene lists for Gene Ontology
Get a summary of the widths of genes.
```{r}
# gene width
summary(width(gr_gene)) 

# estimated number of bubbles by gene gene width * bubble density (approx, uncorrected)
summary(width(gr_gene) *  (nrow(df_sv) / sum(seqlengths(chrom_sizes)))  )
```

### Straightforward yes/no SVs
```{r}
tmp_vector <- df_svHitsByGene_agg %>% filter(density_all > 0) %>% .$geneID
tibble(gene=tmp_vector) %>% arrange(gene) %>% write_delim("./local/gene_ontology/sv_bodyAndFlank.txt", col_names=FALSE)
glue("Genes with SVs in body + 2000bp flanking regions: {length(tmp_vector)}")

tmp_vector <- df_svHitsByGene_agg %>% filter(density_all == 0) %>% .$geneID
tibble(gene=tmp_vector) %>% arrange(gene) %>% write_delim("./local/gene_ontology/sv_bodyAndFlank-None.txt", col_names=FALSE)
glue("Genes with no bubbles in body + 2000bp flanking regions: {length(tmp_vector)}")

tmp_vector <- df_svHitsByGene_agg %>% filter(density_body > 0) %>% .$geneID
tibble(gene=tmp_vector) %>% arrange(gene) %>% write_delim("./local/gene_ontology/sv_bodyOnly.txt", col_names=FALSE)
glue("Genes with SVs in body only: {length(tmp_vector)}")

tmp_vector <- df_svHitsByGene_agg %>% filter(density_body == 0) %>% .$geneID
tibble(gene=tmp_vector) %>% arrange(gene) %>% write_delim("./local/gene_ontology/sv_bodyOnly-None.txt", col_names=FALSE)
glue("Genes with no bubbles: {length(tmp_vector)}")

tmp_vector <- df_svHitsByGene_agg %>% filter(density_flank > 0) %>% .$geneID
tibble(gene=tmp_vector) %>% arrange(gene) %>% write_delim("./local/gene_ontology/sv_flankOnly.txt", col_names=FALSE)
glue("Genes with SVs in flanking regions: {length(tmp_vector)}")

tmp_vector <- df_svHitsByGene_agg %>% filter(density_up > 0) %>% .$geneID
tibble(gene=tmp_vector) %>% arrange(gene) %>% write_delim("./local/gene_ontology/sv_flankUpOnly.txt", col_names=FALSE)
glue("Genes with SVs in upstream regions: {length(tmp_vector)}")

tmp_vector <- df_svHitsByGene_agg %>% filter(density_down > 0) %>% .$geneID
tibble(gene=tmp_vector) %>% arrange(gene) %>% write_delim("./local/gene_ontology/sv_flankDownOnly.txt", col_names=FALSE)
glue("Genes with SVs in downstream regions: {length(tmp_vector)}")
```

### More nuanced (I'll stick with the ones above)
Based on the plot earlier, we know that the median SV density is about 0.215. We can use this to find genes with under- and over-represented bubbles.
```{r}
bd_thresholds <- df_svHitsByGene_agg %>% 
  filter(density_body!=0) %>% .$density_body %>% quantile(c(.3, .50, .7))
bd_thresholds

# genes with fewer SVs than expected
tmp_vector <- df_svHitsByGene_agg %>% filter(density_body < bd_thresholds[1]) %>% .$geneID
tibble(gene=tmp_vector) %>% arrange(gene) %>% write_delim("./local/gene_ontology/sv_bodyOnly-RepresentUnder.txt", col_names=FALSE)
glue("Genes with under-represented bubbles: {length(tmp_vector)}")

# genes with more SVs than expected
tmp_vector <- df_svHitsByGene_agg %>% filter(density_body > bd_thresholds[3]) %>% .$geneID
tibble(gene=tmp_vector) %>% arrange(gene) %>% write_delim("./local/gene_ontology/sv_bodyOnly-RepresentOver.txt", col_names=FALSE)
glue("Genes with over-represented bubbles: {length(tmp_vector)}")

# genes with as much as SVs than expected
tmp_vector <- df_svHitsByGene_agg %>% filter(density_body >= bd_thresholds[1], density_body <= bd_thresholds[3]) %>% .$geneID
tibble(gene=tmp_vector) %>% arrange(gene) %>% write_delim("./local/gene_ontology/sv_bodyOnly-RepresentExpect.txt", col_names=FALSE)
glue("Genes with expected bubbles: {length(tmp_vector)}")
```

## Focusing on only simple bubbles (CAUTION)
Arbitrary criteria about whether the bubble is simple.
```{r}
df_svHitsByGene2 <- df_svHitsByGene %>% 
  left_join(df_sv2 %>% 
    mutate(simple_bubble = n_paths_bio <= 3 & lenbio_max > 100 & lenbio_minmaxratio < 0.05) %>% 
    select(id, simple_bubble)
  )
summary(df_svHitsByGene2$simple_bubble)
```

Here is a version where we compute the bubble density for only simple bubbles. However, I think this is **likely to be inaccurate, because it does not attempt to factor in the complex regions for normalisation.**
```{r}
# version where we only focus on simple bubbles
genes_w_noSV_simple <- names(gr_gene)[!names(gr_gene) %in% unique(filter(df_svHitsByGene2, simple_bubble==TRUE)$geneID)]
df_svHitsByGene_agg_simple <- df_svHitsByGene2 %>%
  filter(simple_bubble == TRUE) %>% 
  group_by(geneID) %>% 
  summarise(
    n_sv_all = n(), 
    n_sv_up = sum(relDistanceToGeneStart<=0),
    n_sv_body = sum(relDistanceToGeneStart>0 & relDistanceToGeneStart<1),
    n_sv_down = sum(relDistanceToGeneStart>=1) ) %>% 
  bind_rows(tibble(geneID = genes_w_noSV_simple, n_sv_all = 0, n_sv_up = 0, n_sv_body = 0, n_sv_down = 0) ) %>%  #add genes with no SVs whatsoever
  arrange(geneID) %>% 
  left_join(tibble(geneID = names(gr_gene), geneLength = width(gr_gene)), by = 'geneID') %>% 
  mutate(
    density_all = n_sv_all / (geneLength + updown_distance) * 1000,
    density_body = n_sv_body / geneLength * 1000,
    density_up = n_sv_up / updown_distance * 1000,
    density_down = n_sv_down / updown_distance * 1000,
    density_flank = (n_sv_up + n_sv_down) / (2*updown_distance) * 1000
  )

glue("mean simple SV density in genic regions\n",
  "overall    : {mean(df_svHitsByGene_agg_simple$density_all)}\n",
  "gene body  : {mean(df_svHitsByGene_agg_simple$density_body)}\n",
  "upstream   : {mean(df_svHitsByGene_agg_simple$density_up)}\n",
  "downstream : {mean(df_svHitsByGene_agg_simple$density_down)}\n",
  "flank      : {mean(df_svHitsByGene_agg$density_flank)}\n"
)
```
Gene counts.
```{r}
glue(
  "SIMPLE\n",
  "Genes with SVs            : {all}\n",
  "Genes with SVs gene body  : {body}\n",
  "Genes with SVs upstream   : {up}\n",
  "Genes with SVs downstream : {down}\n",
  "Genes with SVs in flanks  : {flank}\n",
  tmp_df = df_svHitsByGene_agg_simple,
  all = sum(tmp_df$n_sv_all > 0),
  body = sum(tmp_df$n_sv_body > 0),
  up = sum(tmp_df$n_sv_up > 0),
  down = sum(tmp_df$n_sv_down > 0),
  flank = sum((tmp_df$n_sv_up + tmp_df$n_sv_down) > 0)
)
```

I leave this here in case I wish to revisit this, but the above analysis could be repeated by just focusing on the simple bubbles, because including the more complex ones might not be entirely accurate. Writing the code for this is, however, much more involved especially with regards to differentiating no SV vs no complex SV. I decide that it is currently not worth the time to do, aside from the simple case of identifying genes have simple SVs near/in them.
```{r}
tmp_vector <- df_svHitsByGene_agg_simple %>% filter(density_all > 0) %>% .$geneID
tibble(gene=tmp_vector) %>% arrange(gene) %>% write_delim("./local/gene_ontology/svSimple_bodyAndFlank.txt", col_names=FALSE)
glue("Genes with SVs in body + 2000bp flanking regions: {length(tmp_vector)}")

tmp_vector <- df_svHitsByGene_agg_simple %>% filter(density_body > 0) %>% .$geneID
tibble(gene=tmp_vector) %>% arrange(gene) %>% write_delim("./local/gene_ontology/svSimple_bodyOnly.txt", col_names=FALSE)
glue("Genes with SVs in body only: {length(tmp_vector)}")

tmp_vector <- df_svHitsByGene_agg_simple %>% filter(density_flank > 0) %>% .$geneID
tibble(gene=tmp_vector) %>% arrange(gene) %>% write_delim("./local/gene_ontology/svSimple_flankOnly.txt", col_names=FALSE)
glue("Genes with SVs in flanking regions: {length(tmp_vector)}")
```


# astCal only stuff
## Choosing SVs for PCR Validation 
For PCR validation, we want to focus on SVs which are simple and ideally not too large (under 1500bp was a threshold set by my fellow collaborator in the lab). In addition to that, we want to choose bubbles where there are allelic difference between species we have samples for. Ideally the length difference is between the alleles are as large as possible, making it as close to a complete indel. Another criteria is that the SV be in close proximity to genes and coding exons, because these are easier to design PCR primers for.

In addition to the criteria above, there are two sets of SVs of interest to us.

1. Different between `troMau`, `otoArg`, `copChr`, `rhaChi` and `rhaChi2`. These are the species for which we have the original frozen samples from which the DNA reads for their genome assemblies were created. Running PCR validation would be a direct confirmation that our method works to detect SVs that are present in those individuals.

2. Different between `troMau` and `astCal`. While we no longer have the original samples for the `astCal` genome, both these species are currently grown in our lab, and we plan to detect the presence of the SV in using a sizable number of individuals. Unlike experiment 1, the philosophy is not to validate the presence of the SV, but to see how prevalent a particular SV in these populations. This is because most polymorphisms are not fixed in any species in Lake Malawi, so the presence of a SV in a species in the graph, does not mean that it is necessary present in all individuals of that species.

In the case of the second experiment, we also have RNAseq data from the gonads of male `AstCal` and `troMau` individuals. SVs that are near genes which are differentially expressed between them are also interesting for us, so we have that as an additional criteria.

```{r}
# import differential expression data
df_fe_raw <- read_delim(
  paste0(Sys.getenv("HOME"), "/code/malawi_transposon/cloud/other/rnaseq/gonad/AcalvsTma_DESeq2_wannotation_forFu.txt"),
  col_types = 'cddddddcccccc', na = c("", "NA", "")
)

# some genes have multiple entrezID, we need to deduplicate
df_fe <- df_fe_raw %>% 
  select(gene, baseMean, log2FoldChange, pvalue, padj, description, gene_biotype, external_gene_name) %>%  
  distinct() %>% 
  mutate(diffexpBool = abs(log2FoldChange) > 1.25 & -log10(padj) > 5.5) %>% # stricter than plot
  filter(!gene %in% list_of_dodgy_genes)

list_of_de_genes <- df_fe %>% filter(diffexpBool == TRUE & !gene %in% list_of_dodgy_genes) %>% .$gene

# how many genes are differentially expressed
table(df_fe$diffexpBool)
```


Print differentially expressed gene list, removing the dodgy genes too.
```{r}
df_fe %>% filter(diffexpBool == TRUE) %>% select(gene) %>% write_delim("./local/gene_ontology/diffExp-testes.txt", col_names=FALSE)

list_of_de_genes <- df_fe %>% filter(diffexpBool == TRUE & !gene %in% list_of_dodgy_genes) %>% .$gene
length(list_of_de_genes)
```

Check if our PCR-validated genes are in the gene list
```{r}
# sv pcr
genes_to_validate1 <- "ENSACLG00000020186,ENSACLG00000006581,ENSACLG00000024576,ENSACLG00000017196,ENSACLG00000025203,ENSACLG00000027349,ENSACLG00000022427,ENSACLG00000015018,ENSACLG00000017346,ENSACLG00000002767,ENSACLG00000014403,ENSACLG00000026652,ENSACLG00000027791,ENSACLG00000021983,ENSACLG00000015768,ENSACLG00000013807"
for (gene_item in unlist(strsplit(genes_to_validate1, ',')) ) {
  print(glue("{gene_item}: {gene_item %in% list_of_de_genes} & {gene_item %in% list_of_dodgy_genes}") )
}
print("")

# nice genes related to testis?
genes_to_validate2 <- "ENSACLG00000020324,ENSACLG00000027791,ENSACLG00000000568,ENSACLG00000012444"
for (gene_item in unlist(strsplit(genes_to_validate2, ',')) ) {
  print(glue("{gene_item}: {gene_item %in% list_of_de_genes} & {gene_item %in% list_of_dodgy_genes}") )
}
```

Volcano plot of the differential expression.
```{r volcano_diffExp, fig.width=6.6, fig.height=6}
library(EnhancedVolcano)

tmp_df <- df_fe %>% 
  select(gene, external_gene_name, log2FoldChange, padj) %>% 
  mutate(gene_unify = if_else(is.na(external_gene_name), gene, external_gene_name))

EnhancedVolcano(tmp_df, 
  lab = tmp_df$gene_unify, x = 'log2FoldChange', y = 'padj',
  title = NULL, subtitle = NULL, caption = NULL,
  FCcutoff = 1, pCutoff = 1e-5)

EnhancedVolcano(tmp_df, 
  lab = tmp_df$gene_unify, x = 'log2FoldChange', y = 'padj',
  title = NULL, subtitle = NULL, caption = NULL,
  FCcutoff = 1, pCutoff = 1e-5) +
  xlim(-14, 14) +
  ylim(0, 70)
```

```{r}
tmp_df %>% arrange(desc(abs(log2FoldChange))) %>% head()

tmp_df %>% arrange(padj) %>% head()
```

Apply the criteria described above to score SVs.
```{r}
# add DE information to sv dataframe
df_sv2_pcr <- df_sv2 %>%  
  left_join(select(df_fe, gene, description, log2FoldChange, pvalue), by='gene')

# add columns about whether the concerned species have different alleles.
# na.rm = FALSE asserts that all compared species must have alleles
species_set1 <- paste0(c('troMau', 'otoArg', 'copChr', 'rhaChi', 'rhaChi2'), '_len')
species_set2 <- paste0(c('astCal', 'troMau'), '_len')
df_sv2_pcr$score_diffPath1  <- rowSums(select(df_sv2_pcr, all_of(species_set1)) == df_sv2_pcr$troMau_len, na.rm = FALSE )
df_sv2_pcr$score_diffPath2  <- rowSums(select(df_sv2_pcr, all_of(species_set2)) == df_sv2_pcr$troMau_len, na.rm = FALSE )

# compute score for suitability for PCR
tmp_df_scorePCR <- df_sv2_pcr %>% 
  transmute(
    id,
    n_paths_bio <= 3,
    lenbio_minmaxratio < 0.05,
    lenbio_max < 1500,
    geneDistance < 5000,
    exonDistance < 5000,
    startDistance < 10000,
    endDistance < 10000
)
tmp_df_scorePCR$scorePCR <- select(tmp_df_scorePCR, -id) %>% rowSums(na.rm = TRUE) 
df_sv2_pcr <- left_join(df_sv2_pcr, select(tmp_df_scorePCR, id, scorePCR), by='id')

# compute score for first and second PCR experiments
df_sv2_pcr <- df_sv2_pcr %>% 
  mutate(
    scoreTotal1 = select(., scorePCR, score_diffPath1) %>% 
      transmute(scorePCR, score_diffPath1 = if_else(score_diffPath1!=5, 1, 0)) %>% 
      rowSums(na.rm = TRUE),
    scoreTotal2 = select(., scorePCR, score_diffPath2, log2FoldChange, pvalue) %>% 
      transmute(scorePCR, 
                score_diffPath2 = if_else(score_diffPath2!=2, 1, 0),
                score_diffExp = if_else(abs(log2FoldChange)>1 & -log10(pvalue)>10, 1, 0) ) %>% 
      rowSums(na.rm = TRUE)
  )
```

Plot distribution of total scores.
```{r total_score_distribution, fig.width=4, fig.height=2.8}
ggplot(df_sv2_pcr, aes(factor(scoreTotal1))) + 
  geom_bar() +
  geom_text(aes(label = ..count..), vjust = -0.5, stat = "count", size=rel(2.8)) +
  xlab("Total score (set 1)") +
  scale_y_continuous("Number of bubbles", expand = expansion(mult = c(.05, .1)) )

ggplot(df_sv2_pcr, aes(factor(scoreTotal2))) + 
  geom_bar() +
  geom_text(aes(label = ..count..), vjust = -0.5, stat = "count", size=rel(2.8)) +
  xlab("Total score (set 2)") +
  scale_y_continuous("Number of bubbles", expand = expansion(mult = c(.05, .1)) ) 
```

Write highest scoring entries into text files.
```{r}
# set 1
df_sv2_pcr %>%  filter(scoreTotal1 >= 8) %>%  
  select(id, ends_with('_len'), gene, log2FoldChange, description, ends_with('_path')) %>% 
  write_csv(file = "./local/sv_for_pcr_validation/sv_set1.csv")

# set 2
df_sv2_pcr %>%  filter(scoreTotal1 >= 9) %>% 
  select(id, ends_with('_len'), gene, log2FoldChange, description, ends_with('_path')) %>% 
  write_csv(file = "./local/sv_for_pcr_validation/sv_set2.csv")

# intersect between both, with additional criteria that it is singleton in set 1
df_sv2_pcr %>%
  filter(scoreTotal1 >= 8 & scoreTotal2 >= 9, score_diffPath1==1 | score_diffPath1==4) %>%  
  select(id, ends_with('_len'), gene, log2FoldChange, description, ends_with('_path')) %>% 
  write_csv(file = "./local/sv_for_pcr_validation/sv_set1&2intersect_singletons.csv")

# intersect between both, with additional criteria that it is not singleton in set 1
df_sv2_pcr %>%  
  filter(scoreTotal1 >= 8 & scoreTotal2 >= 9, score_diffPath1!=1, score_diffPath1!=4) %>%  
  select(id, gene, log2FoldChange, description, ends_with('_len'), ends_with('_path')) %>% 
  write_csv(file = "./local/sv_for_pcr_validation/sv_set1&2intersect_prevalent.csv")
```

Finally, I want to extract the unique paths in a format the ODGI can use to extract their FASTA sequence. For example, `<s411594>s543336<s411593` has to become `411594-,543336+,411593-`.
```{r}
extract_gfa_path <- function(x) {
  segment_list <- stringr::str_extract_all(x, "s[0-9]+") %>% 
    unlist() %>% gsub("s", "", .)
  direction_list <- stringr::str_extract_all(x, "[><]") %>% 
    unlist() %>% 
    gsub(">", "+", .) %>% 
    gsub("<", "-", .)
  gfa_path <- paste(paste0(segment_list, direction_list), collapse=',')
  return(gfa_path)
}

df_sv2_pcr %>%
  filter(scoreTotal1 >= 8 & scoreTotal2 >= 9) %>% 
  select(ends_with('_path')) %>% 
  pivot_longer(cols=everything(), values_to="path") %>% 
  select(path) %>% 
  distinct() %>% 
  filter(path != "*") %>% 
  arrange(path) %>%
  mutate(gfa_path = sapply(path, extract_gfa_path),
         P="P", strand='*') %>% 
  select(P, path, gfa_path, strand) %>% 
  write_delim(file = "./local/sv_for_pcr_validation/allele_paths_for_gfa.txt", delim='\t', col_names=FALSE)
```

## More specific Gene Ontology with differential expression
Get a set of genes with simple SV plus minus 2000bp from start (as above), and with additional criterion of differential expression. The code below is quite elaborate, but it is intended to extract genes with SVs in the flanking regions or within the gene body near the start and end.
```{r}
df_svHitsByGene_nearEdges <- df_svHitsByGene2 %>% 
  filter(simple_bubble == TRUE) %>% 
  filter(
    # upstream
    distanceToGeneStart <= 0 |
      
    # downstream  
    distanceToGeneEnd >= 0 |
      
    # in gene body, but near gene start
    (relDistanceToGeneStart>0 & relDistanceToGeneStart < 0.1 & distanceToGeneStart < 1000) | 
    
    # in gene body, but near gene end  
    (relDistanceToGeneEnd>-0.1 & relDistanceToGeneEnd < 0 & distanceToGeneEnd > -1000)
  )
```

Bubbles inside the gene might be misleading due to SVs affecting the RNAseq mapping.
```{r}
tmp_vector <- df_svHitsByGene_nearEdges %>% 
  filter(geneID %in% list_of_de_genes) %>% 
  select(geneID) %>% distinct() %>% .$geneID

tibble(gene=tmp_vector) %>% arrange(gene) %>% 
  write_delim("./local/gene_ontology/DiffExpTestis-svSimple_nearBoundaries.txt", col_names=FALSE)

glue("Genes with SVs in upstream, downstream and near gene boundaries, differentially expressed in testes: {length(tmp_vector)}")
```

Other versions, similar to what was done before.
```{r}
tmp_vector <- df_svHitsByGene_agg_simple %>% filter(density_all > 0, geneID %in% list_of_de_genes) %>% .$geneID
tibble(gene=tmp_vector) %>% arrange(gene) %>% write_delim("./local/gene_ontology/DiffExpTestis-svSimple_bodyAndFlank.txt", col_names=FALSE)
glue("Genes with SVs in body + 2000bp flanking regions: {length(tmp_vector)}")

tmp_vector <- df_svHitsByGene_agg_simple %>% filter(density_body > 0, geneID %in% list_of_de_genes) %>% .$geneID
tibble(gene=tmp_vector) %>% arrange(gene) %>% write_delim("./local/gene_ontology/DiffExpTestis-svSimple_bodyOnly.txt", col_names=FALSE)
glue("Genes with SVs in body only: {length(tmp_vector)}")

tmp_vector <- df_svHitsByGene_agg_simple %>% filter(density_flank > 0, geneID %in% list_of_de_genes) %>% .$geneID
tibble(gene=tmp_vector) %>% arrange(gene) %>% write_delim("./local/gene_ontology/DiffExpTestis-svSimple_flankOnly.txt", col_names=FALSE)
glue("Genes with SVs in flanking regions: {length(tmp_vector)}")

tmp_vector <- df_svHitsByGene_agg_simple %>% filter(density_up > 0, geneID %in% list_of_de_genes) %>% .$geneID
tibble(gene=tmp_vector) %>% arrange(gene) %>% write_delim("./local/gene_ontology/DiffExpTestis-svSimple_flankUp.txt", col_names=FALSE)
glue("Genes with SVs in upstream regions: {length(tmp_vector)}")

tmp_vector <- df_svHitsByGene_agg_simple %>% filter(density_down > 0, geneID %in% list_of_de_genes) %>% .$geneID
tibble(gene=tmp_vector) %>% arrange(gene) %>% write_delim("./local/gene_ontology/DiffExpTestis-svSimple_flankDown.txt", col_names=FALSE)
glue("Genes with SVs in downstream regions: {length(tmp_vector)}")
```

Create Venn diagram to intersect differentially expressed genes and those with SVs.
```{r venn_diffExp, fig.width=2.8, fig.height=2.6}
# genes with SVs in body and flanking regions
set1 <- list_of_de_genes
set2 <- df_svHitsByGene_agg_simple %>% filter(density_all > 0) %>% .$geneID
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("DE" , "with SV"),
  hyper.test=FALSE, total.population = total_pop_size,
  output=FALSE
)
grid.draw(venn_object)
```

## Hybrid data from Miguel
Here we have some differential expression data on hybrids of AstCal and TroMau in the germline.
```{r}
# filter out false genes from DE set
geneList_DEGhybrid_raw <- readLines(
  paste0(Sys.getenv("HOME"), "/code/malawi_transposon/cloud/other/rnaseq/miguel_hybrid/hybrids_DEG.txt")
)
geneList_DEGhybrid <- geneList_DEGhybrid_raw[!geneList_DEGhybrid_raw %in% list_of_dodgy_genes]
length(geneList_DEGhybrid)
```

Write filtered differentially expressed gene list.
```{r}
tibble(gene = geneList_DEGhybrid) %>% write_delim("./local/gene_ontology/diffExp-hybrid.txt", col_names=FALSE)
```

Gene lists.
```{r}
tmp_vector <- df_svHitsByGene_agg_simple %>% filter(geneID %in% geneList_DEGhybrid) %>% select(geneID) %>% distinct() %>% .$geneID
tibble(gene=tmp_vector) %>% arrange(gene) %>% write_delim("./local/gene_ontology/DiffExpHybrid-svSimple_bodyAndFlank.txt", col_names=FALSE)
glue("Genes with SVs in upstream, downstream and near gene boundaries, differentially expressed in hybrids: {length(tmp_vector)}")

tmp_vector <- df_svHitsByGene_agg_simple %>% filter(density_body > 0, geneID %in% geneList_DEGhybrid) %>% select(geneID) %>% distinct() %>% .$geneID
tibble(gene=tmp_vector) %>% arrange(gene) %>% write_delim("./local/gene_ontology/DiffExpHybrid-svSimple_bodyOnly.txt", col_names=FALSE)
glue("Genes with SVs in body only, differentially expressed in hybrids: {length(tmp_vector)}")

tmp_vector <- df_svHitsByGene_agg_simple %>% filter(density_flank > 0, geneID %in% geneList_DEGhybrid) %>% select(geneID) %>% distinct() %>% .$geneID
tibble(gene=tmp_vector) %>% arrange(gene) %>% write_delim("./local/gene_ontology/DiffExpHybrid-svSimple_flankOnly.txt", col_names=FALSE)
glue("Genes with SVs in flanking regions only, differentially expressed in hybrids: {length(tmp_vector)}")
```

## Conclusions
In general, I did not notice any significant enrichment of terms in any of the gene sets. The difficult lies in the fact that almost every gene has an SV associated with it. I don't expect that relating it to TEs will lead to any major improvement. Probably need more power.

