work_representative-non-coding-transcriptome_part-4.Rmd

---
title: "work_representative-non-coding-transcriptome_part-4.Rmd"
author: "KA"
email: kalavatt@fredhutch.org
output:
  html_document:
    toc: yes
    df_print: paged
  html_notebook:
    toc: yes
    toc_float: yes
---
<br />

## Create `data.frame` from various non-coding annotations
### Get situated
#### Code
<details>
<summary><i>Code: Get situated</i></summary>
```{r}
#!/usr/bin/env Rscript

library(ComplexUpset)
library(GenomicRanges)
library(ggplot2)
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()
```
</details>
<br />

### Read in and process annotation files
#### Code
<details>
<summary><i>Code: Read in and process annotation files</i></summary>
```{r}
#!/usr/bin/env Rscript

load_bed_CUTs_SUTs_XUTs <- function(x) {
    y <- readr::read_tsv(x, col_names = FALSE, show_col_types = FALSE)
    colnames(y) <- c("chr", "start", "end", "name_etc", "score", "strand")
    return(y)
}


load_xlsx_CUTs_4x <- function(x) {
    y <- readxl::read_xlsx(x)
    colnames(y) <- c("chr", "start", "end", "name_etc", "strand")
    y$chr <- gsub("chr", "", y$chr)
    y <- y %>% dplyr::arrange(chr, start)
    return(y)
}


load_txt_CUTs_HMM <- function(x) {
    y <- readr::read_tsv(x, show_col_types = FALSE)
    colnames(y) <- c(
        "chr", "start", "end", "name_etc", "strand",
        "avg_covg_WT", "avg_covg_rrp6", "avg_fold"
    )
    y$chr <- gsub("chr", "", y$chr)
    y <- y %>% dplyr::arrange(chr, start)
    return(y)
}


load_gtf_ncRNAs <- function(x) {
    y <- rtracklayer::import(x) %>%
        tibble::as_tibble() %>%
        dplyr::select(-c("score", "phase"))
    colnames(y) <- c(
        "chr", "start", "end", "width", "strand",
        "source", "type", "gene_id", "transcript_id",
        "feature", "orf_classification", "source_id"
    )
    y <- y %>%
        dplyr::arrange(chr, start)
    return(y)
}


load_gtf_NUTs <- function(x) {
    y <- rtracklayer::import(x) %>%
        tibble::as_tibble() %>%
        dplyr::select(-c("score", "phase"))
    colnames(y) <- c(
        "chr", "start", "end", "width", "strand",
        "source", "type", "gene_id", "transcript_id"
    )
    y$chr <- gsub("chr", "", y$chr)
    y <- y %>% dplyr::arrange(chr, start)
    return(y)
}


load_csv_SRATs <- function(x) {
    y <- readr::read_csv(x, skip = 1, show_col_types = FALSE)
    colnames(y) <- c("chr", "start", "end", "info_sense", "strand", "name_etc")
    y$chr <- gsub("chr", "", y$chr)
    y <- y %>% dplyr::arrange(chr, start)
    return(y)
}


p_base <- "./infiles_gtf-gff3/representation"

#  Check
# list.dirs(p_base)
# list.files(paste(p_base, "CUTs_SUTs", sep = "/"))
# list.files(paste(p_base, "CUTs-HMM_CUTs-4X", sep = "/"))
# list.files(paste(p_base, "ncRNAs", sep = "/"))
# list.files(paste(p_base, "NUTs", sep = "/"))
# list.files(paste(p_base, "SRATs", sep = "/"))
# list.files(paste(p_base, "XUTs", sep = "/"))

in_CUTs <- paste(p_base, "CUTs_SUTs/CUTs.coord-R64.bed", sep = "/")
in_CUTs_4x <- paste(p_base, "CUTs-HMM_CUTs-4X/CUTs-4x.xlsx", sep = "/")
in_CUTs_HMM <- paste(p_base, "CUTs-HMM_CUTs-4X/CUTs-HMM.txt.gz", sep = "/")
in_ncRNAs <- paste(p_base, "ncRNAs/ncRNAs.gtf", sep = "/")
in_NUTs <- paste(p_base, "NUTs/NUTs.gtf", sep = "/")
in_SRATs <- paste(p_base, "SRATs/SRATs.csv", sep = "/")
in_SUTs <- paste(p_base, "CUTs_SUTs/SUTs.coord-R64.bed", sep = "/")
in_XUTs <- paste(p_base, "XUTs/XUTs.coord-R64.bed", sep = "/")

#  Check
# file.exists(in_CUTs)
# file.exists(in_SUTs)
# file.exists(in_CUTs_4x)
# file.exists(in_CUTs_HMM)
# file.exists(in_ncRNAs)
# file.exists(in_NUTs)
# file.exists(in_SRATs)
# file.exists(in_XUTs)

CUTs <- load_bed_CUTs_SUTs_XUTs(in_CUTs)
CUTs_4x <- load_xlsx_CUTs_4x(in_CUTs_4x)
CUTs_HMM <- load_txt_CUTs_HMM(in_CUTs_HMM)
ncRNAs <- load_gtf_ncRNAs(in_ncRNAs)
NUTs <- load_gtf_NUTs(in_NUTs)
SRATs <- load_csv_SRATs(in_SRATs)
SUTs <- load_bed_CUTs_SUTs_XUTs(in_SUTs)
XUTs <- load_bed_CUTs_SUTs_XUTs(in_XUTs)
```
</details>
<br />

### Munge the `data.frame`s such that they have the same columns
#### Code
<details>
<summary><i>Code: Munge the `data.frame`s</i></summary>
```{r}
#!/usr/bin/env Rscript

CUTs <- CUTs %>%
    dplyr::mutate(
                 id = stringr::str_sub(CUTs$name_etc, 6, 11),
               type = stringr::str_sub(CUTs$name_etc, 1, 4),
             source = "Xu et al. (Steinmetz), Nature 2009",
        coordinates = "lifted from R56 to R64",
                etc = stringr::str_sub(CUTs$name_etc, 13, nchar(CUTs$name_etc))
    ) %>%
    dplyr::select(-c(name_etc, score)) %>%
    dplyr::arrange(chr, start, strand)

CUTs_4x <- CUTs_4x %>%
    dplyr::mutate(
                 id = name_etc,
               type = "CUTs_4x",
                etc = NA_character_,
             source = "Vera (Dowell), BMC Genomics 2016",
        coordinates = "R64 (no liftOver)"
    ) %>%
    dplyr::select(-name_etc)

CUTs_HMM <- CUTs_HMM %>%
    dplyr::mutate(
                 id = name_etc,
               type = "CUTs_HMM",
             source = "Vera (Dowell), BMC Genomics 2016",
        coordinates = "R64 (no liftOver)"
    ) %>%
    dplyr::select(-name_etc) %>%
    tidyr::unite(
        etc, c(avg_covg_WT, avg_covg_rrp6, avg_fold)
    ) %>%
    dplyr::arrange(chr, start, strand)

ncRNAs <- ncRNAs %>%
    dplyr::mutate(
                 id = gene_id,
               type = feature,
             source = "SGD R64-1-1",
        coordinates = "R64 (no liftOver)",
                etc = NA_character_
    ) %>%
    dplyr::select(-c(
        "width", "gene_id", "transcript_id", "feature", "orf_classification",
        "source_id"
    )) %>%
    dplyr::arrange(chr, start, strand)

NUTs <- NUTs %>%
    dplyr::select(-c("type", "source", "width")) %>%
    dplyr::mutate(
               type = "NUTs",
             source = "Schulz et al. (Cramer), Cell 2013",
        coordinates = "R64 (no liftOver)"
    ) %>%
    dplyr::rename(id = transcript_id, etc = gene_id)
        
SRATs <- SRATs %>%
    dplyr::mutate(
               type = "SRATs",
             source = "Venkatesh et al. (Workman), Nat Comm 2016",
        coordinates = "R64 (no liftOver)"
    ) %>%
    dplyr::rename(id = name_etc, etc = info_sense) %>%
    dplyr::arrange(chr, start, strand)

SUTs <- SUTs %>%
    dplyr::select(-"score") %>%
    dplyr::mutate(
                 id = stringr::str_sub(SUTs$name_etc, 6, 11),
               type = stringr::str_sub(SUTs$name_etc, 1, 4),
             source = "Xu et al. (Steinmetz), Nature 2009",
        coordinates = "lifted from R56 to R64",
                etc = stringr::str_sub(SUTs$name_etc, 13, nchar(SUTs$name_etc))
    ) %>%
    dplyr::select(-name_etc) %>%
    dplyr::arrange(chr, start, strand)

XUTs <- XUTs %>%
    dplyr::select(-"score") %>%
    dplyr::mutate(
               type = paste0(stringr::str_sub(XUTs$name_etc, 1, 3), "s"),
                 id = stringr::str_sub(XUTs$name_etc, 14, nchar(XUTs$name_etc)),
             source = "van Dijk et al. (Thermes, Morillon), Nature 2011",
        coordinates = "lifted from R63 to R64",
           name_etc = NA_character_
    ) %>%
    dplyr::rename(etc = name_etc) %>%
    dplyr::arrange(chr, start, strand)
```
</details>
<br />

### Bind `data.frame`s
#### Code
<details>
<summary><i>Code: Bind `data.frame`s</i></summary>
```{r}
#!/usr/bin/env Rscript

all_pa_ncRNAs <- dplyr::bind_rows(
    CUTs, CUTs_4x, CUTs_HMM, ncRNAs, NUTs, SRATs, SUTs, XUTs
) %>%
    dplyr::arrange(strand, chr, start)
```
</details>
<br />

### Determine, record, and format overlaps between previously annotated ncRNA (pa-ncRNA) features
#### Links records by overlaps
##### Code
<details>
<summary><i>Code: Links records by overlaps</i></summary>
```{r}
#!/usr/bin/env Rscript

# #  Do the linking stratified by 'chr'
# group_chr <- plyr::ddply(
#     ncRNAs,
#     "chr",
#     function(x) { 
#         #  Check if a record should be linked with the previous record
#         y <- c(NA, x$end[-nrow(x)])
#         z <- ifelse(is.na(y), 0, y)
#         z <- cummax(z)
#         z[is.na(y)] <- NA
#         x$previous_end <- z
#         
#         return(x)
#     }
# )
# 
# group_chr$new_group <- is.na(group_chr$previous_end) | 
#     (group_chr$start >= group_chr$previous_end)
# group_chr$group <- cumsum(group_chr$new_group)

#  Do the linking stratified by both 'chr' and 'strand'
group_chr_strand <- plyr::ddply(
    all_pa_ncRNAs,
    c("chr", "strand"),
    function(x) { 
        #  Check if a record should be linked with the previous record
        y <- c(NA, x$end[-nrow(x)])
        z <- ifelse(is.na(y), 0, y)
        z <- cummax(z)
        z[is.na(y)] <- NA
        x$previous_end <- z
        
        return(x)
    }
) #%>%
    #dplyr::arrange(chr, start, strand)
#NOTE This call to arrange was throwing things off

group_chr_strand$new_group <- is.na(group_chr_strand$previous_end) | 
    (group_chr_strand$start >= group_chr_strand$previous_end)
group_chr_strand$group <- cumsum(group_chr_strand$new_group)

# #  Check groups with and without grouping by strand
# group_chr_strand %>% dplyr::select(chr, start, end, strand, group)
# 
# #  Check numbers of groups with and without grouping by strand
# group_chr_strand$group %>% forcats::as_factor() %>% levels() %>% length()
```
</details>
<br />

#### For overlaps by '`chr`', identify new groups on '`previous_end`'
##### Code
<details>
<summary><i>Code: For overlaps by '`chr`', identify new groups on '`previous_end`'</i></summary>
```{r, eval=FALSE}
#!/usr/bin/env Rscript

# run <- FALSE
# if(base::isTRUE(run)) {
#     #  Aggregate the data
#     agg_chr <- plyr::ddply(
#         group_chr,
#         .(chr, group),
#         plyr::summarize,
#         start = min(start),
#         end = max(end),
#         details_type = paste0(type, collapse = "; "),
#         details_all = paste0(type, ": ", id, ", ", strand, collapse = "; ")
#     ) %>%
#         dplyr::arrange(chr, start) %>%
#         dplyr::mutate(strand = ".") %>%
#         dplyr::relocate(
#             c(chr, start, end, strand, details_type, details_all, group)
#         ) %>%
#         dplyr::rename(
#             details_type_init = details_type, details_all_init = details_all
#         )
#     
#     # #  Check
#     # agg_chr$details_type_init %>% head()
#     # agg_chr$details_all_init %>% head()
#     
#     #  stackoverflow.com/questions/49389792/how-to-collapse-all-character-vectors-in-a-list
#     details_type <- agg_chr$details_type %>%
#         stringr::str_remove_all("s") %>%
#         stringr::str_split("; ") %>%
#         purrr::map(sort) %>%
#         purrr::map(unique)
#     details_type <- sapply(
#         details_type[!sapply(details_type, purrr::is_empty)],
#         paste,
#         collapse = " "
#     ) %>%
#         tibble::as_tibble_col(column_name = "details_type")
#     
#     details_all <- agg_chr$details_all %>%
#         stringr::str_remove_all("s") %>%
#         stringr::str_split("; ") %>%
#         purrr::map(sort) %>%
#         purrr::map(unique)
#     details_all <- sapply(
#         details_all[!sapply(details_all, purrr::is_empty)],
#         paste,
#         collapse = "; "
#     ) %>%
#         tibble::as_tibble_col(column_name = "details_all")
#     
#     details_id <- agg_chr$details_all %>%
#         stringr::str_remove_all("s|,\\s\\+|,\\s\\-") %>%
#         stringr::str_replace_all(
#             "(CUT|CUT_4x|CUT_HMM|NUT|SRAT|SUT|XUT):\\s", ""
#         ) %>%
#         stringr::str_split("; ") %>%
#         purrr::map(sort) %>%
#         purrr::map(unique)
#     details_id <- sapply(
#         details_id[!sapply(details_id, purrr::is_empty)],
#         paste,
#         collapse = "; "
#     ) %>%
#         tibble::as_tibble_col(column_name = "details_id")
#     
#     details_strand <- agg_chr$details_all %>%
#         stringr::str_remove_all("s") %>%
#         stringr::str_remove_all(
#             "(CUT|CUT_4x|CUT_HMM|NUT|SRAT|SUT|XUT):\\s|\\;|\\,"
#         ) %>%
#         stringr::str_remove_all(
#             "(CUT|NUT|SRAT|SUT)\\d+\\s|[:digit:]+(F|R)\\-\\d+\\s"
#         ) %>%
#         tibble::as_tibble_col(column_name = "details_strand")
#     
#     n_features <- lengths(
#         gregexpr("(^[A-z]|[1-9])*;", details_id$details_id)
#     ) + 1L
#     
#     agg_chr <- dplyr::bind_cols(
#         agg_chr, details_type, details_id, details_strand, n_features,
#         details_all
#     ) %>%
#         dplyr::select(-group) %>%
#         dplyr::relocate(
#             c(details_type_init, details_all_init),
#             .after = details_all
#         ) %>%
#         dplyr::rename("n_features" = "...11")
#     
#     agg_chr$n_features <- ifelse(
#         (agg_chr$n_features == 2 & agg_chr$details_strand == "+") |
#         (agg_chr$n_features == 2 & agg_chr$details_strand == "-"),
#         1,
#         agg_chr$n_features
#     )
#     
#     rm(details_id, details_type, details_strand, n_features, details_all)
# }
```
</details>
<br />

#### For overlaps by '`chr`' and '`strand`, identify new groups on '`previous_end`'
##### Code
<details>
<summary><i>Code: For overlaps by '`chr`' and '`strand`, identify new groups on '`previous_end`'</i></summary>
```{r}
#!/usr/bin/env Rscript

run <- TRUE
if(base::isTRUE(run)) {
    #  Aggregate the data
    agg_chr_strand <- plyr::ddply(
        group_chr_strand,
        .(chr, strand, group),
        plyr::summarize,
        start = min(start),
        end = max(end),
        details_type = paste0(type, collapse = "; "),
        details_all = paste0(type, ": ", id, collapse = "; ")
    ) %>%
        dplyr::arrange(chr, start, strand) %>%
        dplyr::relocate(
            c(chr, start, end, strand, details_type, details_all, group)
        ) %>%
        dplyr::rename(
            details_type_init = details_type, details_all_init = details_all
        )
    
    # #  Check
    # agg_chr_strand$details_type_init %>% head()
    # agg_chr_strand$details_all_init %>% head()
    
    #  stackoverflow.com/questions/49389792/how-to-collapse-all-character-vectors-in-a-list
    details_type <- agg_chr_strand$details_type %>%
        stringr::str_remove_all("s") %>%
        stringr::str_split("; ") %>%
        purrr::map(sort) %>%
        purrr::map(unique)
    details_type <- sapply(
        details_type[!sapply(details_type, purrr::is_empty)],
        paste,
        collapse = " "
    ) %>%
        tibble::as_tibble_col(column_name = "details_type")
    
    details_all <- agg_chr_strand$details_all %>%
        stringr::str_remove_all("s") %>%
        stringr::str_split("; ") %>%
        purrr::map(sort) %>%
        purrr::map(unique)
    details_all <- sapply(
        details_all[!sapply(details_all, purrr::is_empty)],
        paste,
        collapse = "; "
    ) %>%
        tibble::as_tibble_col(column_name = "details_all")
    
    details_id <- agg_chr_strand$details_all %>%
        stringr::str_remove_all("s") %>%
        stringr::str_replace_all(
            "(CUT|CUT_4x|CUT_HMM|ncRNA|NUT|SRAT|SUT|XUT):\\s", ""
        ) %>%
        stringr::str_split("; ") %>%
        purrr::map(sort) %>%
        purrr::map(unique)
    details_id <- sapply(
        details_id[!sapply(details_id, purrr::is_empty)],
        paste,
        collapse = "; "
    ) %>%
        tibble::as_tibble_col(column_name = "details_id")
    
    n_types <- lengths(
        gregexpr("(^[A-z]|[1-9])* ", details_type$details_type)
    ) + 1L
    
    n_features <- lengths(
        gregexpr("(^[A-z]|[1-9])*;", details_id$details_id)
    ) + 1L
    
    agg_chr_strand <- dplyr::bind_cols(
        agg_chr_strand, details_type, details_id, n_types, n_features,
        details_all
    ) %>%
        dplyr::select(-group) %>%
        dplyr::relocate(
            c(details_type_init, details_all_init),
            .after = details_all
        ) %>%
        dplyr::rename(
            "n_types" = "...10",
            "n_features" = "...11"
        )
    
    #  With the above logic, values with only one feature or type are erroneously
    #+ counted as two, so fix that
    agg_chr_strand$n_types <- ifelse(
        grepl(" ", agg_chr_strand$details_type, fixed = TRUE),
        agg_chr_strand$n_types,
        1L
    )
    
    agg_chr_strand$n_features <- ifelse(
        grepl(";", agg_chr_strand$details_id, fixed = TRUE),
        agg_chr_strand$n_features,
        1L
    )
    
    rm(details_id, details_type, n_features, details_all)
}
```
</details>
<br />

### Examine overlaps between collapsed/merged pa-ncRNAs and SGD
#### Code
`#NOTE` No longer doing this here  
`#TODO` Link to where we are doing this
<details>
<summary><i>Code: Join `data.frame`s</i></summary>

```{r, eval=FALSE}
#!/usr/bin/env Rscript

# in_SGD <- paste("./infiles_gtf-gff3/already/Saccharomyces_cerevisiae.R64-1-1.108.gtf")
# # file.exists(in_SGD)
```
</details>
<br />
<br />

## Collect metrics for numbers of features/types
### Code
<details>
<summary><i>Code: Collect metrics for numbers of features/types</i></summary>
```{r}
#!/usr/bin/env Rscript

agg_chr_strand.bak <- agg_chr_strand
# agg_chr_strand <- agg_chr_strand.bak

agg_chr_strand <- agg_chr_strand %>%
    dplyr::mutate(
        CUT = stringr::str_count(details_type_init, "$CUTs|CUTs$|CUTs;|CUTs$"),
        CUT_4x = stringr::str_count(details_type_init, "CUTs_4x"),
        CUT_HMM = stringr::str_count(details_type_init, "CUTs_HMM"),
        ncRNA = stringr::str_count(details_type_init, "ncRNA"),
        NUT = stringr::str_count(details_type_init, "NUT"),
        SRAT = stringr::str_count(details_type_init, "SRAT"),
        SUT = stringr::str_count(details_type_init, "SUT"),
        XUT = stringr::str_count(details_type_init, "XUT"),
        
        #  Create columns of Boolean values
        b_CUT = as.logical(CUT),
        b_CUT_4x = as.logical(CUT_4x),
        b_CUT_HMM = as.logical(CUT_HMM),
        b_ncRNA = as.logical(ncRNA),
        b_NUT = as.logical(NUT),
        b_SRAT = as.logical(SRAT),
        b_SUT = as.logical(SUT),
        b_XUT = as.logical(XUT),
    )


evaluate_nos_types <- function(x) {
    y <- ifelse(
        rowSums(agg_chr_strand[, 20:27]) == x, TRUE, FALSE
    )
    return(y)
}


evaluate_nos_features <- function(x) {
    y <- ifelse(
        rowSums(agg_chr_strand[, 12:19]) == x, TRUE, FALSE
    )
    return(y)
}


#  Tally row-wise numbers of types
# rowSums(agg_chr_strand[, 20:27]) %>% min()  # [1] 1
# rowSums(agg_chr_strand[, 20:27]) %>% max()  # [1] 7
agg_chr_strand$types_1 <- evaluate_nos_types(1)
agg_chr_strand$types_2 <- evaluate_nos_types(2)
agg_chr_strand$types_3 <- evaluate_nos_types(3)
agg_chr_strand$types_4 <- evaluate_nos_types(4)
agg_chr_strand$types_5 <- evaluate_nos_types(5)
agg_chr_strand$types_6 <- evaluate_nos_types(6)
agg_chr_strand$types_7 <- evaluate_nos_types(7)

#  Tally row-wise numbers of features
# rowSums(agg_chr_strand[, 12:19]) %>% min()  # [1] 1
# rowSums(agg_chr_strand[, 12:19]) %>% max()  # [1] 12
agg_chr_strand$features_1 <- evaluate_nos_features(1)
agg_chr_strand$features_2 <- evaluate_nos_features(2)
agg_chr_strand$features_3 <- evaluate_nos_features(3)
agg_chr_strand$features_4 <- evaluate_nos_features(4)
agg_chr_strand$features_5 <- evaluate_nos_features(5)
agg_chr_strand$features_6 <- evaluate_nos_features(6)
agg_chr_strand$features_7 <- evaluate_nos_features(7)
agg_chr_strand$features_8 <- evaluate_nos_features(8)
agg_chr_strand$features_9 <- evaluate_nos_features(9)
agg_chr_strand$features_10 <- evaluate_nos_features(10)
agg_chr_strand$features_11 <- evaluate_nos_features(11)
agg_chr_strand$features_12 <- evaluate_nos_features(12)

# #  Check
# rowSums(agg_chr_strand[20:27]) %>% table()
# agg_chr_strand$n_types %>% table()
# (rowSums(agg_chr_strand[20:27]) == agg_chr_strand$n_types) %>% table()

all_pa_ncRNAs_cats <- colnames(agg_chr_strand)[20:27]

# #  Check
# agg_chr_strand[, all_pa_ncRNAs_cats] %>% head(15) %>% t()

agg_chr_strand$length <- agg_chr_strand$end - agg_chr_strand$start
agg_chr_strand <- agg_chr_strand %>%
    dplyr::relocate(length, .after = end)
```
</details>
<br />

## Run `ComplexUpset`
### Code
<details>
<summary><i>Code: Run `ComplexUpset`</i></summary>
```{r}
#!/usr/bin/env Rscript

#  Sort by intersection size
ComplexUpset::upset(
    agg_chr_strand,
    all_pa_ncRNAs_cats,
    annotations = list(
        "Length (bp)" = agg_chr_strand %>%
            ggplot2::ggplot(., aes(x = intersection, y = length)) +
            geom_violin(na.rm = TRUE),
        "No. features" = agg_chr_strand %>%
            ggplot2::ggplot(., aes(x = intersection, y = n_features)) +
            geom_violin(na.rm = TRUE),
        "Intersection ratio" = intersection_ratio(
            text_mapping = aes(label = !!upset_text_percentage())
        )
    ),
    # base_annotations =list(
    #     "Intersection size" = intersection_size(),
    #     "Intersection ratio" = intersection_ratio()
    # ),
    name = "Groups of pa-ncRNA categories",
    # sort_intersections_by = "ratio",
    n_intersection = 30,
    # n_intersection = 10,
    # min_size = 8,
    width_ratio = 0.1
)

#  Sort by which intersections are larger than expected by chance (assuming
#+ equal probabilities of belonging to multiple sets) via the intersection
#+ size/union ratio
ComplexUpset::upset(
    agg_chr_strand,
    all_pa_ncRNAs_cats,
    annotations = list(
        "Length (bp)" = agg_chr_strand %>%
            ggplot2::ggplot(., aes(x = intersection, y = length)) +
            geom_violin(na.rm = TRUE),
        "No. features" = agg_chr_strand %>%
            ggplot2::ggplot(., aes(x = intersection, y = n_features)) +
            geom_violin(na.rm = TRUE),
        "Intersection ratio" = intersection_ratio(
            text_mapping = aes(label = !!upset_text_percentage())
        )
    ),
    # base_annotations =list(
    #     "Intersection size" = intersection_size(),
    #     "Intersection ratio" = intersection_ratio()
    # ),
    name = "Groups of pa-ncRNA categories",
    sort_intersections_by = "ratio",
    n_intersection = 30,
    # n_intersection = 10,
    # min_size = 8,
    width_ratio = 0.1
)

#  Sort by set degree
ComplexUpset::upset(
    agg_chr_strand,
    all_pa_ncRNAs_cats,
    annotations = list(
        "Length (bp)" = agg_chr_strand %>%
            ggplot2::ggplot(., aes(x = intersection, y = length)) +
            geom_violin(na.rm = TRUE),
        "No. features" = agg_chr_strand %>%
            ggplot2::ggplot(., aes(x = intersection, y = n_features)) +
            geom_violin(na.rm = TRUE),
        "Intersection ratio" = intersection_ratio(
            text_mapping = aes(label = !!upset_text_percentage())
        )
    ),
    # base_annotations =list(
    #     "Intersection size" = intersection_size(),
    #     "Intersection ratio" = intersection_ratio()
    # ),
    name = "Groups of pa-ncRNA categories",
    sort_intersections_by = "degree",
    n_intersection = 35,
    # n_intersection = 10,
    # min_size = 8,
    width_ratio = 0.1
)
```
</details>
<br />
<br />

## Write out gtf files for individual and merged/collapsed ncRNA feature types
### Create outfile directory if it doesn't exist
#### Code
<details>
<summary><i>Create outfile directory if it doesn't exist</i></summary>
```{r}
#!/usr/bin/env Rscript

representation <- "outfiles_gtf-gff3/representation"
ifelse(
    dir.exists(representation),
    print("Directory exists"),
    dir.create(representation)
)
```
</details>
<br />

### Generate gtf for individual ncRNA feature types
#### Code
<details>
<summary><i>Generate gtf for individual ncRNA feature types</i></summary>
```{r}
#!/usr/bin/env Rscript

write_gtf <- function(x, y) {
    # :param x: tibble
    # :param y: outfile
    # :return: NA
    readr::write_tsv(
        x,
        y,
        col_names = FALSE,
        quote = "none",
        escape = "none"
    )
}


handle_Xu <- function(x) {
    y <- x %>%
        dplyr::rename(seqname = chr) %>%
        dplyr::mutate(feature = "feature") %>%
        dplyr::relocate(c(source, feature), .after = seqname) %>%
        dplyr::mutate(score = ".", frame = ".") %>%
        dplyr::relocate(score, .after = end) %>%
        dplyr::relocate(frame, .after = strand) %>%
        dplyr::mutate(
            attribute = paste(
                paste0("gene_id \"", id, "\""),
                paste0("transcript_id \"", id, "\""),
                paste0("type \"", stringr::str_remove(type, "s$"), "\""),
                paste0(
                    "liftOver \"",
                    stringr::str_replace_all(coordinates, " ", "-"),
                    "\""
                ),
                paste0("source_misc_1 \"", stringr::str_split_i(x$etc, "_", 1), "\""),
                paste0("source_misc_2 \"", stringr::str_split_i(x$etc, "_", 2), "\""),
                paste0("source_misc_3 \"", stringr::str_split_i(x$etc, "_", 3), "\""),
                sep = "; "
            )
        ) %>%
        dplyr::select(-c(id, type, etc, coordinates))
    
    return(y)
}


handle_Vera <- function(x) {
    y <- x %>%
        dplyr::rename(seqname = chr) %>%
        dplyr::mutate(feature = "feature") %>%
        dplyr::relocate(c(source, feature), .after = seqname) %>%
        dplyr::mutate(score = ".", frame = ".") %>%
        dplyr::relocate(score, .after = end) %>%
        dplyr::relocate(frame, .after = strand) %>%
        dplyr::mutate(
            attribute = paste(
                paste0("gene_id \"", id, "\""),
                paste0("transcript_id \"", id, "\""),
                paste0("type \"", stringr::str_remove(type, "s"), "\""),
                paste0(
                    "liftOver \"",
                    stringr::str_replace_all(coordinates, " ", "-"),
                    "\""
                ),
                sep = "; "
            )
        ) %>%
        dplyr::select(-c(id, type, etc, coordinates))
    
    return(y)
}


handle_Schulz <- function(x) {
    y <- x %>%
        dplyr::rename(seqname = chr) %>%
        dplyr::mutate(feature = "feature") %>%
        dplyr::relocate(c(source, feature), .after = seqname) %>%
        dplyr::mutate(score = ".", frame = ".") %>%
        dplyr::relocate(score, .after = end) %>%
        dplyr::relocate(frame, .after = strand) %>%
        dplyr::mutate(
            attribute = paste(
                paste0("gene_id \"", id, "\""),
                paste0("transcript_id \"", id, "\""),
                paste0("type \"", stringr::str_remove(type, "s"), "\""),
                paste0(
                    "liftOver \"",
                    stringr::str_replace_all(coordinates, " ", "-"),
                    "\""
                ),
                paste0("source_misc \"", etc, "\""),
                sep = "; "
            )
        ) %>%
        dplyr::select(-c(id, type, etc, coordinates))
    
    return(y)
}


`handle_SGD_R64-1-1` <- function(x) {
    y <- x %>%
        dplyr::rename(seqname = chr) %>%
        dplyr::mutate(feature = "feature") %>%
        dplyr::relocate(c(source, feature), .after = seqname) %>%
        dplyr::mutate(score = ".", frame = ".") %>%
        dplyr::relocate(score, .after = end) %>%
        dplyr::relocate(frame, .after = strand) %>%
        dplyr::mutate(
            attribute = paste(
                paste0("gene_id \"", id, "\""),
                paste0("transcript_id \"", id, "\""),
                paste0("type \"", type, "\""),
                paste0(
                    "liftOver \"",
                    stringr::str_replace_all(coordinates, " ", "-"),
                    "\""
                ),
                sep = "; "
            )
        ) %>%
        dplyr::select(-c(type, id, coordinates, etc))
    return(y)
}


handle_van_Dijk <- function(x) {
    y <- x %>%
        dplyr::rename(seqname = chr) %>%
        dplyr::mutate(feature = "feature") %>%
        dplyr::relocate(c(source, feature), .after = seqname) %>%
        dplyr::mutate(score = ".", frame = ".") %>%
        dplyr::relocate(score, .after = end) %>%
        dplyr::relocate(frame, .after = strand) %>%
        dplyr::mutate(
            attribute = paste(
                paste0("gene_id \"", id, "\""),
                paste0("transcript_id \"", id, "\""),
                paste0("type \"", stringr::str_remove(type, "s"), "\""),
                paste0(
                    "liftOver \"",
                    stringr::str_replace_all(coordinates, " ", "-"),
                    "\""
                ),
                sep = "; "
            )
        ) %>%
        dplyr::select(-c(id, type, etc, coordinates))
    
    return(y)
}


handle_Venkatesh <- function(x) {
    y <- x %>%
        dplyr::rename(seqname = chr) %>%
        dplyr::mutate(feature = "feature") %>%
        dplyr::relocate(c(source, feature), .after = seqname) %>%
        dplyr::mutate(score = ".", frame = ".") %>%
        dplyr::relocate(score, .after = end) %>%
        dplyr::relocate(frame, .after = strand) %>%
        dplyr::mutate(
            attribute = paste(
                paste0("gene_id \"", id, "\""),
                paste0("transcript_id \"", id, "\""),
                paste0("type \"", stringr::str_remove(type, "s$"), "\""),
                paste0(
                    "liftOver \"",
                    stringr::str_replace_all(coordinates, " ", "-"),
                    "\""
                ),
                paste0("source_misc_1 \"", stringr::str_split_i(x$etc, "_", 1), "\""),
                paste0("source_misc_2 \"", stringr::str_split_i(x$etc, "_", 2), "\""),
                paste0("source_misc_3 \"", stringr::str_split_i(x$etc, "_", 3), "\""),
                paste0("source_misc_4 \"", stringr::str_split_i(x$etc, "_", 4), "\""),
                sep = "; "
            )
        ) %>%
        dplyr::select(-c(id, type, etc, coordinates))
    
    return(y)
}


gtf_CUTs <- handle_Xu(CUTs)
gtf_CUTs_4x <- handle_Vera(CUTs_4x)
gtf_CUTs_HMM <- handle_Vera(CUTs_HMM)
gtf_ncRNAs <- `handle_SGD_R64-1-1`(ncRNAs)
gtf_NUTs <- handle_Schulz(NUTs)
gtf_SRATs <- handle_Venkatesh(SRATs)
gtf_SUTs <- handle_Xu(SUTs)
gtf_XUTs <- handle_van_Dijk(XUTs)

write_gtf(
    gtf_CUTs,
    paste(representation, "Greenlaw-et-al_CUTs.gtf", sep = "/")
)
write_gtf(
    gtf_CUTs_4x,
    paste(representation, "Greenlaw-et-al_CUTs-4x.gtf", sep = "/")
)
write_gtf(
    gtf_CUTs_HMM,
    paste(representation, "Greenlaw-et-al_CUTs-HMM.gtf", sep = "/")
)
write_gtf(
    gtf_ncRNAs,
    paste(representation, "Greenlaw-et-al_ncRNAs.gtf", sep = "/")
)
write_gtf(
    gtf_NUTs,
    paste(representation, "Greenlaw-et-al_NUTs.gtf", sep = "/")
)
write_gtf(
    gtf_SRATs,
    paste(representation, "Greenlaw-et-al_SRATs.gtf", sep = "/")
)
write_gtf(
    gtf_SUTs,
    paste(representation, "Greenlaw-et-al_SUTs.gtf", sep = "/")
)
write_gtf(
    gtf_XUTs,
    paste(representation, "Greenlaw-et-al_XUTs.gtf", sep = "/")
)
```
</details>
<br />

### Generate combined gtf for previously annotated ncRNA (pa-ncRNA) feature types (not merged/collapsed)
#### Code
<details>
<summary><i>Generate combined gtf for previously annotated ncRNA (pa-ncRNA) feature types (not merged/collapsed)</i></summary>
```{r}
#!/usr/bin/env Rscript

#  Combine tibbles of the various individual pa-ncRNA feature types -----------
gtf_combined <- dplyr::bind_rows(
    gtf_CUTs,
    gtf_CUTs_4x,
    gtf_CUTs_HMM,
    gtf_ncRNAs,
    gtf_NUTs,
    gtf_SRATs,
    gtf_SUTs,
    gtf_XUTs
)

#  Add an extra attribute that combines "type" with "gene_id" -----------------
#+ ...thereby avoiding duplicate feature names (e.g., CUT, CUT_HMM, CUT_4x have
#+ some identical feature names even though they are different features) in
#+ counts matrices, subsequent DE analyses, etc. Call the new attribute
#+ category "complete."

#  Parse "attribute" to make a column of feature-wise "type" strings
gtf_combined$type <- gtf_combined$attribute %>%
    stringr::str_split("; ") %>%
    unlist() %>%
    grep("type", ., value = TRUE) %>%
    stringr::str_remove_all("\"|type\ ")

#  Parse "attribute" to make a column of feature-wise "gene_id" strings
gtf_combined$gene_id <- gtf_combined$attribute %>%
    stringr::str_split("; ") %>%
    unlist() %>%
    grep("gene_id", ., value = TRUE) %>%
    stringr::str_remove_all("\"|gene_id\ ")

#  Concatenate the "type" and "gene_id" feature-wise strings, creating strings
#+ of "complete" names in a new "complete" column; then, remove the "type" and
#+ "gene_id" columns
gtf_combined <- gtf_combined %>%
    tidyr::unite(complete, type, gene_id, sep = "-")

#  Add the "complete" column info to "attribute", then remove the "complete"
#+ column
gtf_combined <- gtf_combined %>%
    dplyr::mutate(
        attribute = paste0(attribute, "; complete \"", complete, "\"")
    ) %>%
    dplyr::select(-"complete")

#  Check how the new feature-wise "attribute" entries look
# gtf_combined$attribute %>% head(100)


#  Finally, write out the gtf -------------------------------------------------
write_gtf(
    gtf_combined,
    paste(
        representation,
        "Greenlaw-et-al_non-collapsed-non-coding-transcriptome.gtf",
        sep = "/"
    )
)
```
</details>
<br />

### Generate gtf for merged/collapsed ncRNA feature types
#### Code
<details>
<summary><i>Generate gtf for merged/collapsed ncRNA feature types</i></summary>
```{r}
#!/usr/bin/env Rscript

agg_chr_strand[, 1:12] %>%
    dplyr::rename(seqname = chr) %>%
    dplyr::mutate(
        source = "Greenlaw et al., work_representative-non-coding-transcriptome_part-4.Rmd (chunks 1-7)",
        feature = "feature"
    ) %>%
    dplyr::relocate(c(source, feature), .after = seqname) %>%
    dplyr::mutate(score = ".", frame = ".") %>%
    dplyr::relocate(score, .after = end) %>%
    dplyr::relocate(frame, .after = strand) %>%
    dplyr::mutate(
        attribute = paste(
            paste0(
                "gene_id \"",
                details_all_init %>%
                    stringr::str_replace_all("\\_H", "H") %>%
                    stringr::str_replace_all("\\_4", "4") %>%
                    stringr::str_remove_all("s") %>%
                    stringr::str_replace_all(": ", "-") %>%
                    stringr::str_replace_all("; ", "_"),
                "\""
            ),
            paste0(
                "transcript_id \"",
                details_all_init %>%
                    stringr::str_replace_all("\\_H", "H") %>%
                    stringr::str_replace_all("\\_4", "4") %>%
                    stringr::str_remove_all("s") %>%
                    stringr::str_replace_all(": ", "-") %>%
                    stringr::str_replace_all("; ", "_"),
                "\""
            ),
            paste0(
                "details_type_alpha \"",
                details_type %>%
                    stringr::str_replace_all("\\_H", "H") %>%
                    stringr::str_replace_all("\\_4", "4") %>%
                    stringr::str_replace_all(" ", "_"),
                "\""
            ),
            paste0(
                "details_type \"",
                details_type_init %>%
                    stringr::str_replace_all("\\_H", "H") %>%
                    stringr::str_replace_all("\\_4", "4") %>%
                    stringr::str_remove_all("s") %>%
                    stringr::str_replace_all("; ", "_"),
                "\""
            ),
            paste0(
                "details_id \"",
                details_id %>%
                    stringr::str_replace_all("; ", "_"),
                "\""
            ),
            paste0(
                "details_all \"",
                details_all %>%
                    stringr::str_replace_all("\\_H", "H") %>%
                    stringr::str_replace_all("\\_4", "4") %>%
                    stringr::str_replace_all(": ", "-") %>%
                    stringr::str_replace_all("; ", "_"),
                "\""
            ),
            paste0("n_types \"", n_types, "\""),
            paste0("n_features \"", n_features, "\""),
            paste0("n_types_features \"", n_types, " ", n_features, "\""),
            paste0("length \"", length, "\""),
            sep = "; "
        )
    ) %>%
    dplyr::select(-c(
            length, details_type, details_id, n_types, n_features, details_all,
            details_type_init, details_all_init
    )) %>%
    write_gtf(
        .,
        paste(
            representation,
            "Greenlaw-et-al_representative-non-coding-transcriptome.gtf",
            sep = "/"
        )
    )
```
</details>
<br />

## Next step
Go to [`work_representative-non-coding-transcriptome_part-5.md`](./work_representative-non-coding-transcriptome_part-5.md)
<br />