work_assess-process_R64-1-1-gff3_categorize-Trinity-transfrags_part-0.Rmd

---
title: "work_assess-process_R64-1-1-gff3_categorize-Trinity-transfrags_part-0.Rmd"
author: "KA"
email: "kalavatt@fredhutch.org"
output:
    html_notebook:
        toc: yes
        toc_float: true
---
<br />

Assess the differences between `saccharomyces_cerevisiae_R64-1-1_20110208.gff` [from SGD](http://sgd-archive.yeastgenome.org/sequence/S288C_reference/genome_releases/S288C_reference_genome_R64-1-1_20110203.tgz) and `Saccharomyces_cerevisiae.R64-1-1.108.gff3` [from Ensembl](https://ftp.ensembl.org/pub/release-108/gff3/saccharomyces_cerevisiae/), the file we have been using to generate counts matrices against mRNA features.

## Get situated
### Code
<details>
<summary><i>Code: Get situated</i></summary>

```{r Get situated, results='hide', message=FALSE, warning=FALSE}
#!/usr/bin/env Rscript

library(GenomicRanges)
library(IRanges)
library(plyr)
library(readxl)
library(rtracklayer)
library(tidyverse)

options(scipen = 999)
options(ggrepel.max.overlaps = Inf)

if(stringr::str_detect(getwd(), "kalavattam")) {
    p_local <- "/Users/kalavattam/Dropbox/FHCC"
} else {
    p_local <- "/Users/kalavatt/projects-etc"
}
p_wd <- "2022-2023_RRP6-NAB3/results/2023-0215"

setwd(paste(p_local, p_wd, sep = "/"))
getwd()

rm(p_local, p_wd)
```
</details>
<br />
<br />

## Load `gtf`, `gff3` files
### Code
<details>
<summary><i>Code: Load `gtf`, `gff3` files</i></summary>

```{r Load gtf gff3 files, results='hide', message=FALSE, warning=FALSE}
#!/usr/bin/env Rscript

p_SGD <- "./infiles_gtf-gff3/comprehensive/S288C_reference_genome_R64-1-1_20110203"
f_SGD <- "saccharomyces_cerevisiae_R64-1-1_20110208.gff"
gff3_SGD <- paste(p_SGD, f_SGD, sep = "/") %>%
    rtracklayer::import() %>%
    tibble::as_tibble() %>%
    dplyr::arrange(seqnames, start) %>%
    dplyr::mutate(seqnames = gsub("chr", "", seqnames))

rm(p_SGD, f_SGD)

p_Ens <- "./infiles_gtf-gff3/already"
f_Ens <- "Saccharomyces_cerevisiae.R64-1-1.108.gff3"
gff3_Ens <- paste(p_Ens, f_Ens, sep = "/") %>%
    rtracklayer::import() %>%
    tibble::as_tibble() %>%
    dplyr::arrange(seqnames, start)

rm(p_Ens, f_Ens)

p_Ens <- "./infiles_gtf-gff3/already"
f_Ens <- "Saccharomyces_cerevisiae.R64-1-1.108.gtf"
gtf_Ens <- paste(p_Ens, f_Ens, sep = "/") %>%
    rtracklayer::import() %>%
    tibble::as_tibble() %>%
    dplyr::arrange(seqnames, start)

rm(p_Ens, f_Ens)
```
</details>
<br />
<br />

## Perform comparisons
### Are the numbers of sn/snoRNA features per `gtf`/`gff3` the same?
#### Code
<details>
<summary><i>Code: Are the numbers of sn/snoRNA features per `gtf`/`gff3` the same?</i></summary>

```{r Are the numbers of sn/snoRNA features per gtf/gff3 the same?, results='hide', message=FALSE, warning=FALSE}
#!/usr/bin/env Rscript

cat("gff3_SGD colnames\n")
cat("-----------------\n")
gff3_SGD %>% colnames()
cat("\n")

cat("gff3_Ens colnames\n")
cat("-----------------\n")
gff3_Ens %>% colnames()
cat("\n")

cat("gtf_Ens colnames\n")
cat("----------------\n")
gtf_Ens %>% colnames()
cat("\n")

type_gff3_SGD <- gff3_SGD %>%
    dplyr::group_by(type) %>%
    dplyr::summarize(n = dplyr::n())

type_gff3_Ens <- gff3_Ens %>%
    dplyr::group_by(type) %>% 
    dplyr::summarize(n = dplyr::n())

type_gtf_Ens <- gtf_Ens %>%
    dplyr::group_by(type) %>% 
    dplyr::summarize(n = dplyr::n())

biotype_gff3_SGD <- gff3_SGD %>%
    dplyr::group_by(type) %>%
    dplyr::summarize(n = dplyr::n())

biotype_gff3_Ens <- gff3_Ens %>%
    dplyr::group_by(biotype) %>% 
    dplyr::summarize(n = dplyr::n())

biotype_gtf_Ens <- gtf_Ens %>%
    dplyr::group_by(gene_biotype) %>% 
    dplyr::summarize(n = dplyr::n())

type_gff3_SGD
type_gff3_Ens
type_gtf_Ens
biotype_gff3_SGD
biotype_gff3_Ens
biotype_gtf_Ens
```
</details>
<br />

### Find gene IDs in SGD not in Ensembl
#### Code
<details>
<summary><i>Code: Find gene IDs in SGD not in Ensembl</i></summary>

```{r Find gene IDs in SGD not in Ensembl, results='hide', message=FALSE, warning=FALSE}
#!/usr/bin/env Rscript

#  Initialize function --------------------------------------------------------
`%notin%` <- Negate(`%in%`)


#  For the SGD gff3, subset by column "type" entry "gene" ---------------------
z_gff3_SGD_genes <- gff3_SGD %>%
    dplyr::filter(type == "gene")

#  Tally number of genes per chromosome
# z_gff3_SGD_genes %>%
#     dplyr::group_by(seqnames) %>%
#     dplyr::summarize(n = dplyr::n())

z_gff3_Ens_genes_summary <- gff3_Ens %>%
    dplyr::filter(type == "gene") %>%
    dplyr::group_by(seqnames) %>%
    dplyr::summarize(n = dplyr::n())

z_gff3_Ens_genes_summary_sum <- sum(z_gff3_Ens_genes_summary$n)


#  For the Ensembl gff3, subset by column "type" entry "gene" -----------------
z_gff3_Ens_genes <- gff3_Ens %>%
    dplyr::filter(type == "gene") %>%
    dplyr::group_by(seqnames)

z_gff3_SGD_genes$ID %notin%
    stringr::str_remove(z_gff3_Ens_genes$ID, "gene:") %>%
table()


#  Find gene IDs in SGD not in Ensembl ----------------------------------------
z_gff3_SGD_genes[z_gff3_SGD_genes$ID %notin%
        stringr::str_remove(z_gff3_Ens_genes$ID, "gene:"), ]$ID
# [1] "YER109C" "YFL057C" "YFL056C" "YOR031W" "R0010W"  "R0020C"  "R0030W"  "R0040C"

#  (column ID) ------------------------
c("YER109C", "YFL057C", "YFL056C", "YOR031W") %in%
    stringr::str_remove(gff3_Ens$ID, "gene:|CDS:|transcript:")
# [1] FALSE FALSE  TRUE FALSE

gff3_Ens[stringr::str_remove(gff3_Ens$ID, "gene:|CDS:|transcript:") %in%
    c("YER109C", "YFL057C", "YFL056C", "YOR031W"), ]

#  (column protein_ID) ----------------
c("YER109C", "YFL057C", "YFL056C", "YOR031W") %in% gff3_Ens$protein_id
# [1] FALSE FALSE FALSE FALSE

# gff3_Ens[gff3_Ens$protein_id %in%
#     c("YER109C", "YFL057C", "YFL056C", "YOR031W"), ]

#  (column gene_id) -------------------
c("YER109C", "YFL057C", "YFL056C", "YOR031W") %in% gff3_Ens$gene_id
# [1] FALSE FALSE TRUE FALSE

gff3_Ens[gff3_Ens$gene_id %in%
    c("YER109C", "YFL057C", "YFL056C", "YOR031W"), ]

#  (column transcript_id) -------------
c("YER109C", "YFL057C", "YFL056C", "YOR031W") %in% gff3_Ens$transcript_id
# [1] FALSE FALSE TRUE FALSE

gff3_Ens[gff3_Ens$transcript_id %in%
    c("YER109C", "YFL057C", "YFL056C", "YOR031W"), ]


#  Find any rows associated with "YER109C", "YFL057C", "YFL056C", "YOR031W" ---
library(data.table)
gff3_Ens[apply(gff3_Ens, 1, function(x) any(x %like% "YER109C")), ]  # 0 rows
gff3_Ens[apply(gff3_Ens, 1, function(x) any(x %like% "YFL057C")), ]  # 1 row (in row for "YFL056C")
gff3_Ens[apply(gff3_Ens, 1, function(x) any(x %like% "YFL056C")), ]  # 3 rows
gff3_Ens[apply(gff3_Ens, 1, function(x) any(x %like% "YOR031W")), ]  # 0 rows


#  How many "gene" entries after filtering out "Mito" and "2-micron"? ---------
# z_gff3_SGD_genes %>%
#     dplyr::group_by(seqnames) %>%
#     dplyr::summarize(n = n())
# 
# z_gff3_SGD_genes %>%
#     nrow()  # [1] 6607
# 
# z_gff3_SGD_genes %>%
#     dplyr::filter(seqnames %notin% c("Mito", "2-micron")) %>%
#     nrow()  # [1] 6575
```
</details>
<br />

#### Notes
<details>
<summary><i>Notes: Find gene IDs in SGD not in Ensembl</i></summary>

- "YER109C", "YFL057C", "YFL056C", and "YOR031W" are on, respectively, chromosomes V, VI, VI, and XV; "R0010W", "R0020C", "R0030W", and "R0040C" are on "2-micron". Thus, "YER109C", "YFL057C", "YFL056C", and "YOR031W" are apparently valid genes that we *may* care about that are not&mdash;*for some reason*&mdash;associated with IDs in the Ensembl annotation.
- "YER109C", "YFL057C", "YFL056C", and "YOR031W" are not associated with `protein_ID` assignments in `gff3_Ens`; however, in `gff3_Ens`, "YFL056C" appears to be associated with `type` "pseudogene" and "pseudogene_transcript", `biotype` "pseudogene", and corresponding `gene_id` and `transcript_id` entries; the other genes&mdash;"YER109C", "YFL057C", and "YOR031W"&mdash;are not.
- "YFL056C" is associated with three rows in `gff3_Ens`: `type`s "pseudogene", "pseudogenic_transcript", and "exon".
- "YFL057C" is associated with the "YFL056C" `type` "pseudogene" row (it's in the `description` column)
</details>
<br />

### Examine repeat regions, telomeres, centromeres
#### Code
<details>
<summary><i>Code: Examine repeat regions, telomeres, centromeres</i></summary>

```{r Examine repeat regions telomeres centromeres, results='hide', message=FALSE, warning=FALSE}
#!/usr/bin/env Rscript

#  repeat region --------------------------------------------------------------
z_gff3_SGD_repeat_region <- gff3_SGD %>%
    dplyr::filter(type == "repeat_region")

z_rr_y_NA <- z_gff3_SGD_repeat_region[is.na(z_gff3_SGD_repeat_region$ID), ]
z_rr_n_NA <- z_gff3_SGD_repeat_region[!is.na(z_gff3_SGD_repeat_region$ID), ]

z_rr_y_NA$Parent %>% unlist()
#  [1] "TEL02L-YP"  "TEL04R-YP"  "TEL05L-YP"  "TEL05R-YP"  "TEL06L-YP"  "TEL07R-YP"  "TEL08L-YP"  "TEL08R-YP"  "TEL09L-YP"  "TEL10L-YP"  "TEL12L-YP1" "TEL12L-YP2" "TEL12R-YP2" "TEL12R-YP1"
# [15] "TEL13L-YP"  "TEL14L-YP"  "TEL15R-YP"  "TEL16L-YP"  "TEL16R-YP"

z_rr_n_NA$Parent %>% unlist()
# character(0)


#  telomere -------------------------------------------------------------------
z_gff3_SGD_telomere <- gff3_SGD %>%
    dplyr::filter(type == "telomere")

z_tel_y_NA <- z_gff3_SGD_telomere[is.na(z_gff3_SGD_telomere$ID), ]
z_tel_n_NA <- z_gff3_SGD_telomere[!is.na(z_gff3_SGD_telomere$ID), ]

table(z_tel_n_NA$end %in% z_rr_n_NA$end)
# TRUE
#   32

table(z_tel_n_NA$start %in% z_rr_n_NA$start)
# TRUE
#   32

gr_tel_n_NA <- GenomicRanges::GRanges(z_tel_n_NA)
gr_rr_n_NA <- GenomicRanges::GRanges(z_rr_n_NA)

IRanges::countOverlaps(gr_tel_n_NA, gr_rr_n_NA)
tmp <- IRanges::findOverlaps(gr_tel_n_NA, gr_rr_n_NA) %>% tibble::as_tibble()
rm(tmp)

IRanges::countOverlaps(gr_rr_n_NA, gr_tel_n_NA)
tmp <- IRanges::findOverlaps(gr_rr_n_NA, gr_tel_n_NA) %>% tibble::as_tibble()
rm(tmp)


#  centromere -----------------------------------------------------------------
z_gff3_SGD_centromere <- gff3_SGD %>%
    dplyr::filter(type == "centromere")

gr_cen <- GenomicRanges::GRanges(z_gff3_SGD_centromere)

IRanges::countOverlaps(gr_cen, gr_tel_n_NA)
IRanges::countOverlaps(gr_cen, gr_rr_n_NA)

# rm(z_rr_n_NA, z_rr_y_NA, z_tel_n_NA, z_tel_y_NA)
# rm(z_gff3_SGD_repeat_region)
# rm(z_gff3_SGD_telomere)
```
</details>
<br />

#### Notes
<details>
<summary><i>Notes: Examine repeat regions, telomeres, centromeres</i></summary>

- For `type` "repeat_region", go with `!is.na()`, which encompasses the ranges in `is.na()`.
- All of the ranges in `z_rr_n_NA` overlap the ranges in `z_tel_n_NA`. There are no overlaps between `z_gff3_SGD_centromere` and `z_rr_n_NA` (which might be expected), nor are there overlaps between `z_gff3_SGD_centromere` and `z_tel_n_NA` (which are not expected).
- `#CONCLUSION` Thus, I think I should exclude `type` "repeat_region" in favor of `type` "telomere". "Centromere" can be included without any problems with any features, I think.
</details>
<br />

### Examine ncRNA, noncoding_exon, tRNA, snoRNA, snRNA, rRNA
#### Code
<details>
<summary><i>Code: Examine ncRNA, noncoding_exon, tRNA, snoRNA, snRNA, rRNA</i></summary>

```{r Examine ncRNA noncoding_exon tRNA snoRNA snRNA rRNA, results='hide', message=FALSE, warning=FALSE}
#!/usr/bin/env Rscript

z_gff3_SGD_ncRNA <- gff3_SGD %>%
    dplyr::filter(type == "ncRNA")
z_gff3_SGD_noncoding_exon <- gff3_SGD %>%
    dplyr::filter(type == "noncoding_exon")
z_gff3_SGD_tRNA <- gff3_SGD %>%
    dplyr::filter(type == "tRNA")
z_gff3_SGD_snRNA_snoRNA <- gff3_SGD %>%
    dplyr::filter(type == "snRNA" | type == "snoRNA")
z_gff3_SGD_rRNA <- gff3_SGD %>%
    dplyr::filter(type == "rRNA")

z_gff3_noncoding_etc <- gff3_SGD %>%
    dplyr::filter(
        type == "ncRNA" |
        type == "noncoding_exon" |
        type == "tRNA" |
        type == "snoRNA" |
        type == "snRNA" |
        type == "rRNA"
    )

gr_noncoding_exon <- GenomicRanges::GRanges(z_gff3_SGD_noncoding_exon)
gr_noncoding_ncRNA <- GenomicRanges::GRanges(z_gff3_SGD_ncRNA)
gr_tRNA <- GenomicRanges::GRanges(z_gff3_SGD_tRNA)
gr_snRNA_snoRNA <- GenomicRanges::GRanges(z_gff3_SGD_snRNA_snoRNA)
gr_rRNA <- GenomicRanges::GRanges(z_gff3_SGD_rRNA)

sum_1 <- IRanges::countOverlaps(gr_noncoding_exon, gr_noncoding_ncRNA) %>% sum()
tmp_1 <- IRanges::findOverlaps(gr_noncoding_exon, gr_noncoding_ncRNA) %>%
    tibble::as_tibble()

sum_2 <- IRanges::countOverlaps(gr_noncoding_exon, gr_tRNA) %>% sum()
tmp_2 <- IRanges::findOverlaps(gr_noncoding_exon, gr_tRNA) %>%
    tibble::as_tibble()

sum_3 <- IRanges::countOverlaps(gr_noncoding_exon, gr_snRNA_snoRNA) %>% sum()
tmp_3 <- IRanges::findOverlaps(gr_noncoding_exon, gr_snRNA_snoRNA) %>%
    tibble::as_tibble()

sum_4 <- IRanges::countOverlaps(gr_noncoding_exon, gr_rRNA) %>% sum()
tmp_4 <- IRanges::findOverlaps(gr_noncoding_exon, gr_rRNA) %>%
    tibble::as_tibble()

z_gff3_SGD_noncoding_exon[
    z_gff3_SGD_noncoding_exon$Name %notin%
        z_gff3_SGD_noncoding_exon[c(tmp_1$queryHits, tmp_2$queryHits, tmp_3$queryHits), ]$Name, 
]

nrow(z_gff3_SGD_noncoding_exon) - (sum_1 + sum_2 + sum_3)

z_gff3_SGD_noncoding_exon[tmp_1$queryHits, ]
z_gff3_SGD_noncoding_exon[tmp_2$queryHits, ]
z_gff3_SGD_noncoding_exon[tmp_3$queryHits, ]
z_gff3_SGD_noncoding_exon[tmp_4$queryHits, ]

rm(
    z_gff3_SGD_ncRNA, z_gff3_SGD_noncoding_exon, z_gff3_SGD_tRNA,
    z_gff3_SGD_snRNA_snoRNA, z_gff3_SGD_rRNA
)
```
</details>
<br />

#### Notes
<details>
<summary><i>Notes: Examine ncRNA, noncoding_exon, tRNA, snoRNA, snRNA, rRNA</i></summary>

- `#CONCLUSION` We can exclude `type` "noncoding exon"&mdash;these will be included in "ncRNA", "tRNA", "snoRNA", "snRNA", "rRNA". I think they will not be included with anything else.
- `#CONCLUSION` So, here's what to include (so far) for `type`:
    - "centromere"
    - "genes" (`#MAYBE` stratified by `orf_classification`)
    - "ncRNA"
    - "rRNA"
    - "snRNA"
    - "snoRNA"
    - "telomere"
    - "tRNA"
</details>
<br />

### Examine ARS, LTR, LTR retrotransposon, TE gene
#### Code
<details>
<summary><i>Code: Examine ARS, LTR, LTR retrotransposon, TE gene</i></summary>

```{r Examine ARS LTR LTR-retrotransposon TE-gene, results='hide', message=FALSE, warning=FALSE}
#!/usr/bin/env Rscript

#  Initialize functions -------------------------------------------------------
convert_character_0_NA <- function(x) {
    z <- lapply(
        x, function(y) if(identical(y, character(0))) NA_character_ else y
    ) %>%
        unlist()
    
    return(z)
}


flatten_elements_to_one <- function(x) {
    # For character list elements with two or more subelements, collapse
    # ("flatten") the subelements into a single character element
    # 
    # :param x: <list>
    # :return: character vector of collapsed list elements (list e)
    
    l_collapsed <- x[lengths(x) >= 2] %>% length()
    collapsed <- vector(mode = "character", length = l_collapsed)
    for(i in 1:l_collapsed) {
        # i <- 1
        # cat(i, "\n")
        # cat(x[lengths(x) >= 2][[i]], "\n")
        collapsed[i] <- stringr::str_c(
            x[lengths(x) >= 2][[i]],
            collapse = ", "
        )
    }
    
    return(collapsed)
}


process_list_column <- function(x) {
    x[lengths(x) == 0] <- NA_character_
    if(length(x[lengths(x) >= 2]) != 0) {
        x[lengths(x) >= 2] <- x[lengths(x) >= 2] %>%
            flatten_elements_to_one()
    }
    y <- x %>% unlist()
    return(y)
}


#  Load and process tibbles for features --------------------------------------
z_gff3_SGD_ARS <- gff3_SGD %>%
    dplyr::filter(type == "ARS")
z_gff3_SGD_LTR <- gff3_SGD %>%
    dplyr::filter(type == "long_terminal_repeat")
z_gff3_SGD_LTR_retrotransposon <- gff3_SGD %>%
    dplyr::filter(type == "LTR_retrotransposon")
z_gff3_SGD_TE_gene <- gff3_SGD %>%
    dplyr::filter(type == "transposable_element_gene")


#  Unlist column Note -----------------
z_gff3_SGD_ARS$Note <- unlist(z_gff3_SGD_ARS$Note)
z_gff3_SGD_LTR$Note <- unlist(z_gff3_SGD_LTR$Note)
z_gff3_SGD_LTR_retrotransposon$Note <- unlist(z_gff3_SGD_LTR_retrotransposon$Note)
z_gff3_SGD_TE_gene$Note <- unlist(z_gff3_SGD_TE_gene$Note)


#  Unlist column Parent ---------------
z_gff3_SGD_ARS$Parent <- convert_character_0_NA(z_gff3_SGD_ARS$Parent)
z_gff3_SGD_LTR$Parent <- convert_character_0_NA(z_gff3_SGD_LTR$Parent)
z_gff3_SGD_LTR_retrotransposon$Parent <- convert_character_0_NA(z_gff3_SGD_LTR_retrotransposon$Parent)
z_gff3_SGD_TE_gene$Parent <- convert_character_0_NA(z_gff3_SGD_TE_gene$Parent)


#  Unlist column Alias ----------------
z_gff3_SGD_ARS$Alias <- process_list_column(z_gff3_SGD_ARS$Alias)
z_gff3_SGD_LTR$Alias <- process_list_column(z_gff3_SGD_LTR$Alias)
z_gff3_SGD_LTR_retrotransposon$Alias <- process_list_column(z_gff3_SGD_LTR_retrotransposon$Alias)
z_gff3_SGD_TE_gene$Alias <- process_list_column(z_gff3_SGD_TE_gene$Alias)


#  Unlist column Ontology_term --------
z_gff3_SGD_ARS$Ontology_term <- process_list_column(z_gff3_SGD_ARS$Ontology_term)
z_gff3_SGD_LTR$Ontology_term <- process_list_column(z_gff3_SGD_LTR$Ontology_term)
z_gff3_SGD_LTR_retrotransposon$Ontology_term <- process_list_column(z_gff3_SGD_LTR_retrotransposon$Ontology_term)
z_gff3_SGD_TE_gene$Ontology_term <- process_list_column(z_gff3_SGD_TE_gene$Ontology_term)


#  Check tibbles ----------------------
z_gff3_SGD_ARS
z_gff3_SGD_LTR
z_gff3_SGD_LTR_retrotransposon
z_gff3_SGD_TE_gene

z_gff3_SGD_ARS %>%
    dplyr::group_by(Note) %>%
    dplyr::summarize(n = dplyr::n())

z_gff3_SGD_LTR %>%
    dplyr::group_by(Note) %>%
    dplyr::summarize(n = dplyr::n())

z_gff3_SGD_LTR_retrotransposon %>%
    dplyr::group_by(Note) %>%
    dplyr::summarize(n = dplyr::n())

z_gff3_SGD_TE_gene %>%
    dplyr::group_by(Note) %>%
    dplyr::summarize(n = dplyr::n())


#  ------
gr_ARS <- GenomicRanges::GRanges(z_gff3_SGD_ARS)
gr_LTR <- GenomicRanges::GRanges(z_gff3_SGD_LTR)
gr_LTR_retrotransposon <- GenomicRanges::GRanges(z_gff3_SGD_LTR_retrotransposon)
gr_TE_gene <- GenomicRanges::GRanges(z_gff3_SGD_TE_gene)

#  1 ----
z_gff3_SGD_LTR %>% dplyr::arrange(seqnames, start)
z_gff3_SGD_LTR_retrotransposon %>% dplyr::arrange(seqnames, start)

sum_1 <- IRanges::countOverlaps(gr_LTR, gr_LTR_retrotransposon) %>% sum()
tmp_1 <- IRanges::findOverlaps(gr_LTR, gr_LTR_retrotransposon) %>%
    tibble::as_tibble()

z_gff3_SGD_LTR[tmp_1$queryHits, ]
z_gff3_SGD_LTR_retrotransposon[tmp_1$subjectHits, ]

z_gff3_SGD_LTR_retrotransposon[
    seq(1, nrow(z_gff3_SGD_LTR_retrotransposon)) %notin% unique(tmp_1$subjectHits), 
]

#  2 -----
z_gff3_SGD_LTR %>% dplyr::arrange(seqnames, start) # %>% nrow()
z_gff3_SGD_TE_gene %>% dplyr::arrange(seqnames, start)

sum_2 <- IRanges::countOverlaps(gr_LTR, gr_TE_gene) %>% sum()
tmp_2 <- IRanges::findOverlaps(gr_LTR, gr_TE_gene) %>%
    tibble::as_tibble()

z_gff3_SGD_LTR[tmp_2$queryHits, ]
z_gff3_SGD_TE_gene[tmp_2$subjectHits, ]

z_gff3_SGD_LTR[
    seq(1, nrow(z_gff3_SGD_LTR)) %notin% unique(tmp_2$subjectHits), 
] %>%
    nrow()

#  3 -----
z_gff3_SGD_LTR_retrotransposon %>% dplyr::arrange(seqnames, start) # %>% nrow()
z_gff3_SGD_TE_gene %>% dplyr::arrange(seqnames, start) # %>% nrow()

sum_3 <- IRanges::countOverlaps(gr_LTR_retrotransposon, gr_TE_gene) %>% sum()
tmp_3 <- IRanges::findOverlaps(gr_LTR_retrotransposon, gr_TE_gene) %>%
    tibble::as_tibble()

z_gff3_SGD_LTR_retrotransposon[tmp_3$queryHits, ]
z_gff3_SGD_TE_gene[tmp_3$subjectHits, ]

z_gff3_SGD_LTR_retrotransposon[
    seq(1, nrow(z_gff3_SGD_LTR_retrotransposon)) %notin% unique(tmp_3$subjectHits), 
]

#  4 -----
z_gff3_SGD_LTR
z_gff3_SGD_ARS

sum_4 <- IRanges::countOverlaps(gr_LTR, gr_ARS) %>% sum()
tmp_4 <- IRanges::findOverlaps(gr_LTR, gr_ARS) %>%
    tibble::as_tibble()

z_gff3_SGD_LTR[tmp_4$queryHits, ]
z_gff3_SGD_ARS[tmp_4$subjectHits, ]

z_gff3_SGD_LTR[
    seq(1, nrow(z_gff3_SGD_LTR)) %notin% unique(tmp_4$subjectHits), 
]

#  ------
z_gff3_SGD_LTR[tmp_4$queryHits, ]
z_gff3_SGD_ARS[tmp_4$subjectHits, ]
z_gff3_SGD_LTR[seq(1, nrow(z_gff3_SGD_LTR)) %notin% tmp_4$queryHits, ]
z_gff3_SGD_ARS[seq(1, nrow(z_gff3_SGD_ARS)) %notin% tmp_4$subjectHits, ]
```
</details>
<br />

#### Notes
<details>
<summary><i>Notes: Examine ARS, LTR, LTR retrotransposon, TE gene</i></summary>

- Almost all `z_gff3_SGD_LTR_retrotransposon` overlaps `z_gff3_SGD_LTR`&mdash;all except...
    - `IV 878303 884220 5918 - SGD LTR_retrotransposon ... YDRCTy1-2 SGD:S000006860 YDRCTy1-2 ...`
- 82 out of 89 TEs (`z_gff3_SGD_TE_gene`) at least partially overlap 301 out of 383 LTRs (`z_gff3_SGD_LTR`)
- All 89 TEs (`z_gff3_SGD_TE_gene`) at least partially overlap all but one row in `z_gff3_SGD_LTR_retrotransposon`
- 20 LTRs overlap ARSs
</details>
<br />

### Examine "mystery" categories in SGD `gff3`
#### Code
<details>
<summary><i>Code: Examine "mystery" categories in SGD `gff3`</i></summary>

```{r Examine "mystery" categories in SGD gff3, results='hide', message=FALSE, warning=FALSE}
#!/usr/bin/env Rscript

gff3_mystery <- gff3_SGD %>%
    dplyr::filter(
        type == "nucleotide_match" |
        type == "binding_site" |
        type == "region" |
        type == "five_prime_UTR_intron" |
        type == "gene_cassette" |
        type == "insertion" |
        type == "external_transcribed_spacer_region" |
        type == "internal_transcribed_spacer_region"
    )

#  Unlist column Note
gff3_mystery$Note <- process_list_column(gff3_mystery$Note)

#  Unlist column Parent
gff3_mystery$Parent <- process_list_column(gff3_mystery$Parent)

#  Unlist column Alias
gff3_mystery$Alias <- process_list_column(gff3_mystery$Alias)

#  Unlist column Ontology_term
gff3_mystery$Ontology_term <- process_list_column(gff3_mystery$Ontology_term)


gff3_SGD %>% dplyr::filter(dbxref == "SGD:S000029655")
gff3_SGD %>% dplyr::filter(ID == "HMR")
gff3_SGD %>% dplyr::filter(Name == "YJR027W")
gff3_SGD %>% dplyr::filter(Name == "YBL092W")

# gff3_SGD %>% dplyr::filter(type == "intron")
#  Yeah, let's exclude introns, which will be included along with gene or tRNA
#+ (or etc.) entries
```
</details>
<br />

#### Notes
<details>
<summary><i>Notes: Examine "mystery" categories in SGD `gff3`</i></summary>

- We can exclude `type` "nucleotide_match", which mostly seem to be parts of telomere and ARS features
- We can exclude `type` "binding_site", which are all within telomere features
- We can also exclude `type` "region", which seem to be encompassed in other features
- Retain `type` "five_prime_UTR_intron"
- ~~Retain `type` "gene_cassette"~~ *After discussion with Toshi, we realize that these features will be encompassed by other features (e.g., within `type` "gene"), so we can exclude `type` "gene_cassette"*
- ~~Retain `type` "external_transcribed_spacer_region" and "internal_transcribed_spacer_region"~~ *Actually, don't b/c these are already covered by `type` "rRNA"*
</details>
<br />
<br />

## Survey notes taken so far, then make decisions
### Notes
<details>
<summary><i>Notes: Survey notes taken so far, then make decisions</i></summary>

- Go with the SGD annotation for this work.
- We should not remove any entries from `type` "gene" despite there being an extra four chromosomal genes in the SGD annotation versus the Ensembl R64-1-1 annotation.
- We should exclude `type` "repeat_region" in favor of `type` "telomere". "Centromere" can be included without any problems with any features.
- We can exclude `type` "noncoding exon" since these will be included in variously "ncRNA", "tRNA", "snoRNA", "snRNA", "rRNA".
- `#DEKHO` `#IMPORTANT` So, here's what to include for `type`:
    - ~~"ARS" (both strands)~~
    - ~~"centromere" (both strands)~~
    - "gene" (antisense too)
    - "ncRNA" (antisense too)
    - "pseudogene" (those that are not collapsed/merged into "transposable_element" below; antisense too)
    - "rRNA" (antisense too)
    - "snRNA" (antisense too)
    - "snoRNA" (antisense too)
    - ~~"telomere" (both strands)~~
    - "transposable_element" (made from collapsing/merging "long_terminal_repeat", "LTR_retrotransposon", "transposable_element_gene" and "pseudogene"; antisense too)
    - "tRNA" (antisense too)
- And here's what not to include for `type`:
    - "binding_site"
    - "CDS"
    - "chromosome"
    - "external_transcribed_spacer_region" (already in `type` "rRNA")
    - "five_prime_UTR_intron" (already in `type` "gene" and `type` intergenic)
    - "gene_cassette" (already in `type` "gene" and `type` "intergenic")
    - "insertion"
    - "internal_transcribed_spacer_region" (already in `type` "rRNA")
    - "intron"
    - "noncoding_exon"
    - "nucleotide_match"
    - "region"
    - "repeat_region"
</details>
<br />
<br />

## Next step
Go to [`work_assess-process_R64-1-1_gff3_part-1.Rmd`](./work_assess-process_R64-1-1_gff3_part-1.Rmd)
<br />