work_assessment-processing_gtfs_part-1_Trinity.Rmd

---
title: "work_assessment-processing_gtfs_part-1_Trinity.Rmd"
author: "KA"
email: "kalavatt@fredhutch.org"
output:
    html_notebook:
        toc: yes
        toc_float: true
---
<br />
<br />

## Experiment gameplan
### Text
<details>
<summary><i>Text: Experiment gameplan</i></summary>

1. Use `GffRead` to filter `Trinity`-derived `GMAP`-aligned `gff3` files.
2. Use `htseq-count` to perform feature counting with respect to `type`s "locus", "mRNA", "exon", and "CDS" features from "base" calls to `GffRead`, and with respect to `type` "gene" from "intron-filtering-only" calls to `GffRead`.
3. Read in resulting counts matrices and `gtf` files using custom `R` code; using a mix of `bioconductor` and `tidyverse` programs, perform processing as necessary.
4. Regarding raw, non-normalized alignment-to-putative transcript/transcript fragment (transfrag), retain putative transcripts/transfrags with tallies greater than or equal to a given percentile; in this case, generate results for the 5th to 95th percentiles in iterations of five.
6. Use ~~`rtracklayer`~~ *custom code* to write the converted, collapsed, counts-filtered data objects as `gtf` files. *(Custom code is necessary as some files output by `rtracklayer` have formatting errors that break downstream analyses.)*
7. Send the processed `gtf` files to Alison for inspection, analyses.

For both Q and G1 cell states, Alison asked me to filter the assemblies generated with `Trinity` parameter `mkc` = 1, 4, 8. However, for the sake of completeness, and to understand clearly what "incorrectness" looks like, I will filter `mkc` = 2, 16, 32 as well. Remember that the `Trinity` parameters `mir`, `mg`, and `gf` were held at default levels for the specific assemblies to surveyed by Alison.
</details>
<br />
<br />

## Additional details
### Text, code
<details>
<summary><i>Text, code: Additional details</i></summary>

#### I. On calling `gffread`
```{bash, eval=FALSE}
#!/bin/bash

#  "Base" call to GffRead
gffread \
    -v \
    -g "${fasta_g}" \
    -i 1000 \
    -Z \
    -M -K -Q \
    -F -N -P \
    --force-exons --gene2exon \
    -o "${out}" \
    <(awk -F '\t' 'BEGIN {OFS = FS} { gsub(/chr/, "", $1); gsub(/M/, "Mito", $1); print }' "${in}") \
         > >(tee -a "${err_out%.}.stdout.txt") \
        2> >(tee -a "${err_out%.}.stderr.txt")

#  "Intron-filtering-only" call to GffRead
#+ 
#+ No collapsing, merging with these files; only filtering to exclude exonic
#+ features with introns greater than 1000 bp in length
gffread \
    -v -O \
    -i 1000 \
    -o "${out/.gff3/-intron-filtering-only.gff3}" \
    <(awk -F '\t' 'BEGIN {OFS = FS} { gsub(/chr/, "", $1); gsub(/M/, "Mito", $1); print }' "${in}") \
         > >(tee -a "${err_out%.}-intron-filtering-only.stdout.txt") \
        2> >(tee -a "${err_out%.}-intron-filtering-only.stderr.txt")
```

##### Meaning of parameters
###### "Base" call to `GffRead`
- `-v` expose (warn about) duplicate transcript IDs and other potential problems with the given GFF/GTF records
- `-g` full path to a fasta file with the genomic sequences for all input mappings (one per genomic sequence, with file names matching sequence names) (#NOTE i.e., the fasta - file output from running Trinity in genome-guided mode)
- `-i` discard transcripts having an intron larger than 1000 bp
- `-Z` merge very close exons into a single exon (when intron size<4)
- `-M` cluster the input transcripts into loci, discarding "redundant" transcripts (those with the same exact introns and fully contained or equal boundaries)
- `-K` for the -M option, also discard as redundant the shorter, fully contained transcripts (intron chains matching a part of the container)
- `-Q` for the -M option, no longer require boundary containment when assessing redundancy (can be combined with -K); only introns have to match for multi-exon transcripts, and >=80% overlap for single-exon transcripts
- `-F` keep all GFF attributes (for non-exon features)
- `-N` discard multi-exon mRNAs that have any intron with a non-canonical splice site consensus (i.e., not GT-AG, GC-AG or AT-AC)
- `-P` add transcript level GFF attributes about the coding status of each transcript, including partialness or in-frame stop codons (requires -g)
- `--force-exons` make sure that the lowest level GFF features are considered "exon" features (#NOTE This is standard in gff3 and gtf files)
- `-o` write the output records into <outfile> instead of stdout

###### "Intron-filtering-only" call to `GffRead`
- `-v` expose (warn about) duplicate transcript IDs and other potential problems with the given GFF/GTF records
- `-O` process other non-transcript GFF records (by default non-transcript records are ignored) (`#NOTE` This keeps 'gene' features (output by gmap) in the gff3/gtf, although "gene" is exactly the same as "mRNA" as output by gmap)
- `-i` discard transcripts having an intron larger than 1000 bp
- `-o` write the output records into <outfile> instead of stdout
<br />

#### II. On calling `htseq-count`
##### ...with respect to `type` "locus"
As a result of calling `gffread` with the `-Z` `-M` `-K` and `-Q` arguments, a new feature is added to the gff3/gtf files: "locus". "locus" is a putative parent feature made from collapsing and joining "children" features (hierarchically, "mRNA" and "exon" features) as described for `-Z` `-M` `-K` and `-Q`. This is an extension of the kind of collapsing performed by `gffcompare -C`. Thus, counts matrices were created with respect to these "locus" features.

```{bash, eval=FALSE}
#!/bin/bash

sbatch \
    --job-name="htseq-count-locus" \
    --nodes=1 \
    --cpus-per-task=8 \
    --error="${err_out}-locus.%A.stderr.txt" \
    --output="${err_out}-locus.%A.stdout.txt" \
    htseq-count \
        --order "pos" \
        --stranded "${hc_strd}" \
        --nonunique "all" \
        --type "locus" \
        --idattr "ID" \
        --nprocesses 8 \
        --counts_output "${out/.tsv/-locus.tsv}" \
        --with-header \
        ${bams[*]} \
        "${in}"
```

##### ...with respect to `type` "mRNA"
Create the counts matrix by looking at children "mRNA" features rather than parent "locus" features. This is the kind of collapsing performed by `gffcompare -C`.

```{bash, eval=FALSE}
#!/bin/bash

sbatch \
    --job-name="htseq-count-mRNA" \
    --nodes=1 \
    --cpus-per-task=8 \
    --error="${err_out}-mRNA.%A.stderr.txt" \
    --output="${err_out}-mRNA.%A.stdout.txt" \
    htseq-count \
        --order "pos" \
        --stranded "${hc_strd}" \
        --nonunique "all" \
        --type "mRNA" \
        --idattr "ID" \
        --nprocesses 8 \
        --counts_output "${out/.tsv/-mRNA.tsv}" \
        --with-header \
        ${bams[*]} \
        "${in}"
```

##### ...with respect to `type` "exon"
Create the counts matrix by looking at children "exon" features.

```{bash, eval=FALSE}
#!/bin/bash

sbatch \
    --job-name="htseq-count-exon" \
    --nodes=1 \
    --cpus-per-task=8 \
    --error="${err_out}-exon.%A.stderr.txt" \
    --output="${err_out}-exon.%A.stdout.txt" \
    htseq-count \
        --order "pos" \
        --stranded "${hc_strd}" \
        --nonunique "all" \
        --type "exon" \
        --idattr "Parent" \
        --nprocesses 8 \
        --counts_output "${out/.tsv/-exon.tsv}" \
        --with-header \
        ${bams[*]} \
        "${in}"
```

##### ...with respect to `type` "CDS"
`gffread` scans the `fasta`s supplied `-g` and roughly estimates start and stop codons from the nt sequences; it then adds these CDS estimates to the `gff3`/`gtf` files. Create the counts matrix by looking at these estimated "CDS" features.

```{bash, eval=FALSE}
#!/bin/bash

sbatch \
    --job-name="htseq-count-CDS" \
    --nodes=1 \
    --cpus-per-task=8 \
    --error="${err_out}-CDS.%A.stderr.txt" \
    --output="${err_out}-CDS.%A.stdout.txt" \
    htseq-count \
        --order "pos" \
        --stranded "${hc_strd}" \
        --nonunique "all" \
        --type "CDS" \
        --idattr "Parent" \
        --nprocesses 8 \
        --counts_output "${out/.tsv/-CDS.tsv}" \
        --with-header \
        ${bams[*]} \
        "${in}"
```

##### ...with respect to `type` "gene" ("intron-only-filtering" call to `GffRead`)
Create the counts matrix by looking at parent "gene" features from the files output in section <b>"On calling `gffread`"</b> above.

```{bash, eval=FALSE}
#!/bin/bash
sbatch \
    --job-name="htseq-count-gene" \
    --nodes=1 \
    --cpus-per-task=8 \
    --error="${err_out}-gene.%A.stderr.txt" \
    --output="${err_out}-gene.%A.stdout.txt" \
    htseq-count \
        --order "pos" \
        --stranded "${hc_strd}" \
        --nonunique "all" \
        --type "gene" \
        --idattr "ID" \
        --nprocesses 8 \
        --counts_output "${out/.tsv/-gene.tsv}" \
        --with-header \
        ${bams[*]} \
        "${in}"
```
<br />

#### III. The general workflow for making the filtered gtfs
1. Read in and process `gtf`/`gff3` files (as dataframes)
2. Read in and process counts matrices: "locus", "mRNA", "exon", "CDS", and "gene" (from "intron-only-filtering")
3. Categorize putative transcripts/transfrags by percentile: from percentile 5 to percentile 95 in steps of 5
4. Subset dataframes to retain percentile-filtered IDs: "locus", "mRNA", "exon", "CDS", and "gene" (from "intron-only-filtering")
5. Write out processed `gtf`s
<br />

#### IV. Copying over files to Alison
```{bash, eval=FALSE}
#!/bin/bash

mkdir -p /home/kalavatt/tsukiyamalab/alisong/KA.gtfs_quantile-filtered_Trinity-GG_G1-Q_2023-0410/{gtf-gff3,tsv}
# mkdir: created directory '/home/kalavatt/tsukiyamalab/alisong/KA.gtfs_quantile-filtered_Trinity-GG_G1-Q_2023-0410'
# mkdir: created directory '/home/kalavatt/tsukiyamalab/alisong/KA.gtfs_quantile-filtered_Trinity-GG_G1-Q_2023-0410/gtf-gff3'
# mkdir: created directory '/home/kalavatt/tsukiyamalab/alisong/KA.gtfs_quantile-filtered_Trinity-GG_G1-Q_2023-0410/tsv'

cd /home/kalavatt/tsukiyamalab/kalavatt/2022-2023_RRP6-NAB3/results/2023-0215

cp -r \
    outfiles_gtf-gff3/Trinity-GG \
    /home/kalavatt/tsukiyamalab/alisong/KA.gtfs_quantile-filtered_Trinity-GG_G1-Q_2023-0410/gtf-gff3

cp -r \
    outfiles_htseq-count/Trinity-GG \
    /home/kalavatt/tsukiyamalab/alisong/KA.gtfs_quantile-filtered_Trinity-GG_G1-Q_2023-0410/tsv

#  Manually copying in the 'filtered/' subdirectories from local to remote
```
</details>
<br />
<br />

## 0&ndash;2. Command-line work
See [`work_assessment-processing_gtfs_part-0.md`](./work_assessment-processing_gtfs_part-0.md), where *steps 0 through 2 (and an extra step "2.5")* are detailed. Below, we pick up with *step 3*.
<br />
<br />

## 3. Perform quantile filtering of putative transcripts/transfrags
### Get situated
#### Code
<details>
<summary><i>Code: Get situated</i></summary>

```{r}
#!/usr/bin/env Rscript

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

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

if(base::isTRUE(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 />

### Read in and process `gtf` files (as `data.frame`s)
#### Code
<details>
<summary><i>Code: Read in and process `gtf` files (as `data.frame`s)</i></summary>

```{r}
#!/usr/bin/env Rscript

#  Initialize functions ---------------
import_g <- function(file) rtracklayer::import(file)


read_in_gtfs <- function(vector_files) {
    list <- sapply(
        vector_files,
        import_g,
        simplify = FALSE,
        USE.NAMES = TRUE
    )
    names(list) <- paste0("k", c(1, 16, 2, 32, 4, 8))
    
    df <- sapply(
        list,
        as.data.frame,
        simplify = FALSE,
        USE.NAMES = TRUE
    )
    
    return(df)
}


#  G1 ---------------------------------
path_gtf_G <- "outfiles_gtf-gff3/Trinity-GG/G_N"
files_G <- list.files(
    path_gtf_G,
    pattern = "*.gffread.gff3",
    full.names = TRUE
)
l_G_gtf <- read_in_gtfs(files_G)

files_G_introns_filtered <- list.files(
    path_gtf_G,
    pattern = "*.gffread-intron-filtering-only.gff3",
    full.names = TRUE
)
l_G_gtf_introns_filtered <- read_in_gtfs(files_G_introns_filtered)


#  Q ----------------------------------
path_gtf_Q <- "outfiles_gtf-gff3/Trinity-GG/Q_N"
files_Q <- list.files(
    path_gtf_Q,
    pattern = "*.gffread.gff3",
    full.names = TRUE
)
l_Q_gtf <- read_in_gtfs(files_Q)

files_Q_introns_filtered <- list.files(
    path_gtf_Q,
    pattern = "*.gffread-intron-filtering-only.gff3",
    full.names = TRUE
)
l_Q_gtf_introns_filtered <- read_in_gtfs(files_Q_introns_filtered)


#  Comparisons ------------------------
run <- FALSE
if(base::isTRUE(run)) {
    l_G_gtf[["k1"]] %>% head()
    l_G_gtf[["k2"]] %>% head()

    l_G_gtf[["k1"]] %>% head()
    l_Q_gtf[["k1"]] %>% head()

    l_Q_gtf[["k1"]] %>% head()
    l_Q_gtf[["k2"]] %>% head()
}
```
</details>
<br />

### Read in and process counts matrices
#### Code
<details>
<summary><i>Code: Read in and process counts matrices</i></summary>

```{r}
#!/usr/bin/env Rscript

#  Initialize functions ---------------
read_in_mat <- function(vector_of_files) {
    out_list <- sapply(
        vector_of_files,
        readr::read_tsv,
        simplify = FALSE,
        USE.NAMES = TRUE
    )
    names(out_list) <- paste0("k", c(1, 16, 2, 32, 4, 8))
    
    return(out_list)
}


rename_columns <- function(list) {
    lapply(
        list,
        function(df) {
            names(df) <-
                gsub("bams_renamed\\/UT_prim_UMI\\/WT_", "", names(df)) %>%
                gsub("_day1_ovn_", "_", .) %>%
                gsub("_day7_ovn_", "_", .) %>% 
                gsub("_aux-F_tc-F_", "_", .) %>%
                gsub("_tech1\\.UT_prim_UMI\\.bam", "", .) %>%
                gsub("^\\.\\.\\.1", "id", .)
            
            return(df)
        }
    )
}


#  G1 ---------------------------------
path_mat_G <- "outfiles_htseq-count/Trinity-GG/G_N"

files_G_locus <- list.files(
    path_mat_G, pattern = "*-locus.tsv", full.names = TRUE
)
l_G_mat_locus <- read_in_mat(files_G_locus) %>%
    suppressMessages()

files_G_mRNA <- list.files(
    path_mat_G, pattern = "*-mRNA.tsv", full.names = TRUE
)
l_G_mat_mRNA <- read_in_mat(files_G_mRNA) %>%
    suppressMessages()

files_G_exon <- list.files(
    path_mat_G, pattern = "*-exon.tsv", full.names = TRUE
)
l_G_mat_exon <- read_in_mat(files_G_exon) %>%
    suppressMessages()

files_G_CDS <- list.files(
    path_mat_G, pattern = "*-CDS.tsv", full.names = TRUE
)
l_G_mat_CDS <- read_in_mat(files_G_CDS) %>%
    suppressMessages()

files_G_introns_filtered <- list.files(
    path_mat_G,
    pattern = "*.gffread-intron-filtering-only.hc-strd-eq-gene.tsv",
    full.names = TRUE
)
l_G_mat_introns_filtered <- read_in_mat(files_G_introns_filtered) %>%
    suppressMessages()


#  Q ----------------------------------
path_mat_Q <- "outfiles_htseq-count/Trinity-GG/Q_N"

files_Q_locus <- list.files(
    path_mat_Q, pattern = "*-locus.tsv", full.names = TRUE
)
l_Q_mat_locus <- read_in_mat(files_Q_locus) %>%
    suppressMessages()

files_Q_mRNA <- list.files(
    path_mat_Q, pattern = "*-mRNA.tsv", full.names = TRUE
)
l_Q_mat_mRNA <- read_in_mat(files_Q_mRNA) %>%
    suppressMessages()

files_Q_exon <- list.files(
    path_mat_Q, pattern = "*-exon.tsv", full.names = TRUE
)
l_Q_mat_exon <- read_in_mat(files_Q_exon) %>%
    suppressMessages()

files_Q_CDS <- list.files(
    path_mat_Q, pattern = "*-CDS.tsv", full.names = TRUE
)
l_Q_mat_CDS <- read_in_mat(files_Q_CDS) %>%
    suppressMessages()

files_Q_introns_filtered <- list.files(
    path_mat_Q,
    pattern = "*.gffread-intron-filtering-only.hc-strd-eq-gene.tsv",
    full.names = TRUE
)
l_Q_mat_introns_filtered <- read_in_mat(files_Q_introns_filtered) %>%
    suppressMessages()


#  Comparisons ------------------------
run <- FALSE
if(base::isTRUE(run)) {
    l_Q_mat_sam[["k1"]] %>% head()
    l_Q_mat_sam[["k2"]] %>% head()

    l_Q_mat_sam[["k1"]] %>% head()
    l_Q_mat_rev[["k1"]] %>% head()
}


#  Rename columns ---------------------
#+ ...in each dataframe composing each list
l_G_mat_locus <- rename_columns(l_G_mat_locus)
l_G_mat_mRNA <- rename_columns(l_G_mat_mRNA)
l_G_mat_exon <- rename_columns(l_G_mat_exon)
l_G_mat_CDS <- rename_columns(l_G_mat_CDS)
l_G_mat_introns_filtered <- rename_columns(l_G_mat_introns_filtered)

l_Q_mat_locus <- rename_columns(l_Q_mat_locus)
l_Q_mat_mRNA <- rename_columns(l_Q_mat_mRNA)
l_Q_mat_exon <- rename_columns(l_Q_mat_exon)
l_Q_mat_CDS <- rename_columns(l_Q_mat_CDS)
l_Q_mat_introns_filtered <- rename_columns(l_Q_mat_introns_filtered)
```
</details>
<br />

### Categorize putative transcripts/transcript fragments by percentile
#### Code
<details>
<summary><i>Code: Categorize putative transcripts/transcript fragments by percentile</i></summary>

```{r}
#!/usr/bin/env Rscript

#  Initialize function ----------------
filter_transfrags <- function(x, y) {
    df <- x
    
    if(y == "G1" | y == "G") {
        df <- x[, 1:3] %>% head(., -5)
    } else if(y == "Q") {
        df <- x[, c(1, 6, 7)] %>% head(., -5)
    } else if(y != "G1" | y != "G" | y != "Q") {
        stop("Parameter y must be either 'G1', 'G', or 'Q'")
    }
    
    df[, 4] <- df[, 2] + df[, 3]
    colnames(df)[4] <- "sum"

    ids <- list()

    ids$percentiles <- quantile(
        df$sum,
        probs = c(
            0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50,
            0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95
        )
    ) %>%
        round()
    
    df$below_05 <- ifelse(df$sum <= ids$percentiles[1], TRUE, FALSE)
    df$below_10 <- ifelse(df$sum <= ids$percentiles[2], TRUE, FALSE)
    df$below_15 <- ifelse(df$sum <= ids$percentiles[3], TRUE, FALSE)
    df$below_20 <- ifelse(df$sum <= ids$percentiles[4], TRUE, FALSE)
    df$below_25 <- ifelse(df$sum <= ids$percentiles[5], TRUE, FALSE)
    df$below_30 <- ifelse(df$sum <= ids$percentiles[6], TRUE, FALSE)
    df$below_35 <- ifelse(df$sum <= ids$percentiles[7], TRUE, FALSE)
    df$below_40 <- ifelse(df$sum <= ids$percentiles[8], TRUE, FALSE)
    df$below_45 <- ifelse(df$sum <= ids$percentiles[9], TRUE, FALSE)
    df$below_50 <- ifelse(df$sum <= ids$percentiles[10], TRUE, FALSE)
    df$below_55 <- ifelse(df$sum <= ids$percentiles[11], TRUE, FALSE)
    df$below_60 <- ifelse(df$sum <= ids$percentiles[12], TRUE, FALSE)
    df$below_65 <- ifelse(df$sum <= ids$percentiles[13], TRUE, FALSE)
    df$below_70 <- ifelse(df$sum <= ids$percentiles[14], TRUE, FALSE)
    df$below_75 <- ifelse(df$sum <= ids$percentiles[15], TRUE, FALSE)
    df$below_80 <- ifelse(df$sum <= ids$percentiles[16], TRUE, FALSE)
    df$below_85 <- ifelse(df$sum <= ids$percentiles[17], TRUE, FALSE)
    df$below_90 <- ifelse(df$sum <= ids$percentiles[18], TRUE, FALSE)
    df$below_95 <- ifelse(df$sum <= ids$percentiles[19], TRUE, FALSE)
    
    ids$below_05 <- df[df$below_05 == FALSE, ]$id
    ids$below_10 <- df[df$below_10 == FALSE, ]$id
    ids$below_15 <- df[df$below_15 == FALSE, ]$id
    ids$below_20 <- df[df$below_20 == FALSE, ]$id
    ids$below_25 <- df[df$below_25 == FALSE, ]$id
    ids$below_30 <- df[df$below_30 == FALSE, ]$id
    ids$below_35 <- df[df$below_35 == FALSE, ]$id
    ids$below_40 <- df[df$below_40 == FALSE, ]$id
    ids$below_45 <- df[df$below_45 == FALSE, ]$id
    ids$below_50 <- df[df$below_50 == FALSE, ]$id
    ids$below_55 <- df[df$below_55 == FALSE, ]$id
    ids$below_60 <- df[df$below_60 == FALSE, ]$id
    ids$below_65 <- df[df$below_65 == FALSE, ]$id
    ids$below_70 <- df[df$below_70 == FALSE, ]$id
    ids$below_75 <- df[df$below_75 == FALSE, ]$id
    ids$below_80 <- df[df$below_80 == FALSE, ]$id
    ids$below_85 <- df[df$below_85 == FALSE, ]$id
    ids$below_90 <- df[df$below_90 == FALSE, ]$id
    ids$below_95 <- df[df$below_95 == FALSE, ]$id
    
    return(ids)
}


#  G1 ---------------------------------
f_G_mat_locus <- sapply(
    l_G_mat_locus,
    filter_transfrags,
    simplify = FALSE,
    USE.NAMES = TRUE,
    y = "G1"
)
f_G_mat_mRNA <- sapply(
    l_G_mat_mRNA,
    filter_transfrags,
    simplify = FALSE,
    USE.NAMES = TRUE,
    y = "G1"
)
f_G_mat_exon <- sapply(
    l_G_mat_exon,
    filter_transfrags,
    simplify = FALSE,
    USE.NAMES = TRUE,
    y = "G1"
)
f_G_mat_CDS <- sapply(
    l_G_mat_CDS,
    filter_transfrags,
    simplify = FALSE,
    USE.NAMES = TRUE,
    y = "G1"
)
f_G_mat_introns_filtered <- sapply(
    l_G_mat_introns_filtered,
    filter_transfrags,
    simplify = FALSE,
    USE.NAMES = TRUE,
    y = "G1"
)

#  Q ----------------------------------
f_Q_mat_locus <- sapply(
    l_Q_mat_locus,
    filter_transfrags,
    simplify = FALSE,
    USE.NAMES = TRUE,
    y = "Q"
)
f_Q_mat_mRNA <- sapply(
    l_Q_mat_mRNA,
    filter_transfrags,
    simplify = FALSE,
    USE.NAMES = TRUE,
    y = "Q"
)
f_Q_mat_exon <- sapply(
    l_Q_mat_exon,
    filter_transfrags,
    simplify = FALSE,
    USE.NAMES = TRUE,
    y = "Q"
)
f_Q_mat_CDS <- sapply(
    l_Q_mat_CDS,
    filter_transfrags,
    simplify = FALSE,
    USE.NAMES = TRUE,
    y = "Q"
)
f_Q_mat_introns_filtered <- sapply(
    l_Q_mat_introns_filtered,
    filter_transfrags,
    simplify = FALSE,
    USE.NAMES = TRUE,
    y = "Q"
)
```
</details>
<br />

### Subset `gtf` `data.frame` to retain percentile-filtered `gene_id`s
#### Code
<details>
<summary><i>Code: Subset `gtf` `data.frame` to retain percentile-filtered `gene_id`s</i></summary>

```{r}
#!/usr/bin/env Rscript

#  Initialize functions ---------------
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) {
    # ...
    # 
    # :param x: <list>
    # :return: ...
    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)
}


subset_dataframes <- function(x, y, z) {
    # ...
    # 
    # :param x: list of *.gff3 files <list>
    # :param y: list of putative transcripts categorized by percentile <list>
    # :param z: feature type to examine; must be "locus", "mRNA", "exon", "CDS", or "gene" <chr>
    
    #  Tests
    # x <- l_Q_gtf
    # y <- f_Q_mat_CDS
    # z <- "CDS"
    # 
    # x <- l_G_gtf_introns_filtered
    # y <- f_G_mat_introns_filtered
    # z <- "gene"
    # 
    # x <- l_Q_gtf
    # y <- f_Q_mat_exon
    # z <- "exon"
    
    if((base::isFALSE(stringr::str_detect(z, "locus|mRNA|exon|CDS|gene")))) {
        stop("z must be \"locus\", \"mRNA\", \"exon\", \"CDS\", or \"gene\"")
    }
    
    df <- list()
    mat <- list()
    out <- list()
    
    for(i in names(x)) {
        #  Tests
        # i <- "k1"
        # i <- "k4"
        
        mat[[i]] <- y[[i]]
        
        if(base::isTRUE(z == "locus")) {
            df[[i]] <- x[[i]][x[[i]]["type"] == "locus", ]
        } else if(base::isTRUE(z == "mRNA")) {
            df[[i]] <- x[[i]][x[[i]]["type"] == "mRNA", ]
        } else if(base::isTRUE(z == "exon")) {
            df[[i]] <- x[[i]][x[[i]]["type"] == "exon", ]
        } else if(base::isTRUE(z == "CDS")) {
            df[[i]] <- x[[i]][x[[i]]["type"] == "CDS", ]
        } else if(base::isTRUE(z == "gene")) {
            df[[i]] <- x[[i]][x[[i]]["type"] == "gene", ]
        } 
        
        if(base::isTRUE(stringr::str_detect(z, "locus|mRNA|gene"))) {
            out[[i]]$below_05 <- df[[i]][df[[i]]$ID %in% mat[[i]]$below_05, ]
            out[[i]]$below_10 <- df[[i]][df[[i]]$ID %in% mat[[i]]$below_10, ]
            out[[i]]$below_15 <- df[[i]][df[[i]]$ID %in% mat[[i]]$below_15, ]
            out[[i]]$below_20 <- df[[i]][df[[i]]$ID %in% mat[[i]]$below_20, ]
            out[[i]]$below_25 <- df[[i]][df[[i]]$ID %in% mat[[i]]$below_25, ]
            out[[i]]$below_30 <- df[[i]][df[[i]]$ID %in% mat[[i]]$below_30, ]
            out[[i]]$below_35 <- df[[i]][df[[i]]$ID %in% mat[[i]]$below_35, ]
            out[[i]]$below_40 <- df[[i]][df[[i]]$ID %in% mat[[i]]$below_40, ]
            out[[i]]$below_45 <- df[[i]][df[[i]]$ID %in% mat[[i]]$below_45, ]
            out[[i]]$below_50 <- df[[i]][df[[i]]$ID %in% mat[[i]]$below_50, ]
            out[[i]]$below_55 <- df[[i]][df[[i]]$ID %in% mat[[i]]$below_55, ]
            out[[i]]$below_60 <- df[[i]][df[[i]]$ID %in% mat[[i]]$below_60, ]
            out[[i]]$below_65 <- df[[i]][df[[i]]$ID %in% mat[[i]]$below_65, ]
            out[[i]]$below_70 <- df[[i]][df[[i]]$ID %in% mat[[i]]$below_70, ]
            out[[i]]$below_75 <- df[[i]][df[[i]]$ID %in% mat[[i]]$below_75, ]
            out[[i]]$below_80 <- df[[i]][df[[i]]$ID %in% mat[[i]]$below_80, ]
            out[[i]]$below_85 <- df[[i]][df[[i]]$ID %in% mat[[i]]$below_85, ]
            out[[i]]$below_90 <- df[[i]][df[[i]]$ID %in% mat[[i]]$below_90, ]
            out[[i]]$below_95 <- df[[i]][df[[i]]$ID %in% mat[[i]]$below_95, ]
        } else if(base::isTRUE(stringr::str_detect(z, "exon|CDS"))) {
            out[[i]]$below_05 <- df[[i]][
                process_list_column(df[[i]]$Parent) %in% mat[[i]]$below_05, 
            ]
            out[[i]]$below_10 <- df[[i]][
                process_list_column(df[[i]]$Parent) %in% mat[[i]]$below_10, 
            ]
            out[[i]]$below_15 <- df[[i]][
                process_list_column(df[[i]]$Parent) %in% mat[[i]]$below_15, 
            ]
            out[[i]]$below_20 <- df[[i]][
                process_list_column(df[[i]]$Parent) %in% mat[[i]]$below_20, 
            ]
            out[[i]]$below_25 <- df[[i]][
                process_list_column(df[[i]]$Parent) %in% mat[[i]]$below_25, 
            ]
            out[[i]]$below_30 <- df[[i]][
                process_list_column(df[[i]]$Parent) %in% mat[[i]]$below_30, 
            ]
            out[[i]]$below_35 <- df[[i]][
                process_list_column(df[[i]]$Parent) %in% mat[[i]]$below_35, 
            ]
            out[[i]]$below_40 <- df[[i]][
                process_list_column(df[[i]]$Parent) %in% mat[[i]]$below_40, 
            ]
            out[[i]]$below_45 <- df[[i]][
                process_list_column(df[[i]]$Parent) %in% mat[[i]]$below_45, 
            ]
            out[[i]]$below_50 <- df[[i]][
                process_list_column(df[[i]]$Parent) %in% mat[[i]]$below_50, 
            ]
            out[[i]]$below_55 <- df[[i]][
                process_list_column(df[[i]]$Parent) %in% mat[[i]]$below_55, 
            ]
            out[[i]]$below_60 <- df[[i]][
                process_list_column(df[[i]]$Parent) %in% mat[[i]]$below_60, 
            ]
            out[[i]]$below_65 <- df[[i]][
                process_list_column(df[[i]]$Parent) %in% mat[[i]]$below_65, 
            ]
            out[[i]]$below_70 <- df[[i]][
                process_list_column(df[[i]]$Parent) %in% mat[[i]]$below_70, 
            ]
            out[[i]]$below_75 <- df[[i]][
                process_list_column(df[[i]]$Parent) %in% mat[[i]]$below_75, 
            ]
            out[[i]]$below_80 <- df[[i]][
                process_list_column(df[[i]]$Parent) %in% mat[[i]]$below_80, 
            ]
            out[[i]]$below_85 <- df[[i]][
                process_list_column(df[[i]]$Parent) %in% mat[[i]]$below_85, 
            ]
            out[[i]]$below_90 <- df[[i]][
                process_list_column(df[[i]]$Parent) %in% mat[[i]]$below_90, 
            ]
            out[[i]]$below_95 <- df[[i]][
                process_list_column(df[[i]]$Parent) %in% mat[[i]]$below_95, 
            ]
        }
    }
    
    return(out)
}


#  G1 ---------------------------------
out_G_gtf_locus <- subset_dataframes(l_G_gtf, f_G_mat_locus, "locus")
out_G_gtf_mRNA <- subset_dataframes(l_G_gtf, f_G_mat_mRNA, "mRNA")
out_G_gtf_exon <- subset_dataframes(l_G_gtf, f_G_mat_exon, "exon")
out_G_gtf_CDS <- subset_dataframes(l_G_gtf, f_G_mat_CDS, "CDS")
out_G_gtf_introns_filtered <- subset_dataframes(
    l_G_gtf_introns_filtered, f_G_mat_introns_filtered, "gene"
)


#  Q ----------------------------------
out_Q_gtf_locus <- subset_dataframes(l_Q_gtf, f_Q_mat_locus, "locus")
out_Q_gtf_mRNA <- subset_dataframes(l_Q_gtf, f_Q_mat_mRNA, "mRNA")
out_Q_gtf_exon <- subset_dataframes(l_Q_gtf, f_Q_mat_exon, "exon")
out_Q_gtf_CDS <- subset_dataframes(l_Q_gtf, f_Q_mat_CDS, "CDS")
out_Q_gtf_introns_filtered <- subset_dataframes(
    l_Q_gtf_introns_filtered, f_Q_mat_introns_filtered, "gene"
)
```
</details>
<br />

### Write out processed `gtf`s
#### Code
<details>
<summary><i>Code: Write out processed `gtf`s</i></summary>

```{r}
#!/usr/bin/env Rscript

#  Initialize functions -------------------------------------------------------
write_gtf <- function(x, y) {
    # ...
    # 
    # :param x: tibble <tibble, data.frame>
    # :param y: outfile <chr>
    # :return: NA
    readr::write_tsv(
        x,
        y,
        col_names = FALSE,
        quote = "none",
        escape = "none"
    )
}


#  Create outdirectories ------------------------------------------------------
#  G1 ---------------------------------
path_G_locus <- "outfiles_gtf-gff3/Trinity-GG/G_N/filtered/locus"
if(base::isFALSE(dir.exists(path_G_locus))) dir.create(path_G_locus, recursive = TRUE)

path_G_mRNA <- "outfiles_gtf-gff3/Trinity-GG/G_N/filtered/mRNA"
if(base::isFALSE(dir.exists(path_G_mRNA))) dir.create(path_G_mRNA, recursive = TRUE)

path_G_exon <- "outfiles_gtf-gff3/Trinity-GG/G_N/filtered/exon"
if(base::isFALSE(dir.exists(path_G_exon))) dir.create(path_G_exon, recursive = TRUE)

path_G_CDS <- "outfiles_gtf-gff3/Trinity-GG/G_N/filtered/CDS"
if(base::isFALSE(dir.exists(path_G_CDS))) dir.create(path_G_CDS, recursive = TRUE)

path_G_introns_filtered <- "outfiles_gtf-gff3/Trinity-GG/G_N/filtered/introns_filtered"
if(base::isFALSE(dir.exists(path_G_introns_filtered))) dir.create(path_G_introns_filtered, recursive = TRUE)

#  Q ----------------------------------
path_Q_locus <- "outfiles_gtf-gff3/Trinity-GG/Q_N/filtered/locus"
if(base::isFALSE(dir.exists(path_Q_locus))) dir.create(path_Q_locus, recursive = TRUE)

path_Q_mRNA <- "outfiles_gtf-gff3/Trinity-GG/Q_N/filtered/mRNA"
if(base::isFALSE(dir.exists(path_Q_mRNA))) dir.create(path_Q_mRNA, recursive = TRUE)

path_Q_exon <- "outfiles_gtf-gff3/Trinity-GG/Q_N/filtered/exon"
if(base::isFALSE(dir.exists(path_Q_exon))) dir.create(path_Q_exon, recursive = TRUE)

path_Q_CDS <- "outfiles_gtf-gff3/Trinity-GG/Q_N/filtered/CDS"
if(base::isFALSE(dir.exists(path_Q_CDS))) dir.create(path_Q_CDS, recursive = TRUE)

path_Q_introns_filtered <- "outfiles_gtf-gff3/Trinity-GG/Q_N/filtered/introns_filtered"
if(base::isFALSE(dir.exists(path_Q_introns_filtered))) dir.create(path_Q_introns_filtered, recursive = TRUE)


#  Write out gtfs using bespoke approaches/code -------------------------------
for(i in names(out_G_gtf_locus)) {
    for(j in c(paste0(
        "below_",
        c(
            "05", "10", "15", "20", "25", "30", "35", "40", "45", "50",
            "55", "60", "65", "70", "75", "80", "85", "90", "95"
        )
    ))) {
        for(k in c(
            "path_G_locus", "path_G_mRNA", "path_G_exon", "path_G_CDS",
            "path_G_introns_filtered",
            "path_Q_locus", "path_Q_mRNA", "path_Q_exon", "path_Q_CDS",
            "path_Q_introns_filtered"
        )) {
            # i <- "k4"
            # j <- "below_25"
            # k <- "path_G_introns_filtered"
            if(base::isTRUE(unlist(stringr::str_split(k, "_"))[2] == "G")) {
                prefix_file <- "/G1_mkc-"
                prefix_message <- "Saving G1"
            } else if(base::isTRUE(unlist(stringr::str_split(k, "_"))[2] == "Q")) {
                prefix_file <- "/Q_mkc-"
                prefix_message <- "Saving Q"
            }
            outfile <- paste0(
                eval(parse(text = k)), prefix_file, gsub("k", "", i), "_",
                gsub("below_", "gte-pctl-", j), ".gtf"
            )
            
            cat(prefix_message, i, j, "in", k, "\n")
            cat(paste0(outfile, "\n"))
            cat("\n\n")
            
            if(base::isTRUE(unlist(stringr::str_split(k, "_"))[3] == "locus")) {
                if(base::isTRUE(unlist(stringr::str_split(k, "_"))[2] == "G")) {
                    gtf <- out_G_gtf_locus[[i]][[j]]
                }
                if(base::isTRUE(unlist(stringr::str_split(k, "_"))[2] == "Q")) {
                    gtf <- out_Q_gtf_locus[[i]][[j]]
                }
                gtf <- gtf %>%
                    tibble::as_tibble() %>%
                    dplyr::select(-c(
                        "phase", "width", "geneID", "Name", "coverage",
                        "identity", "matches", "mismatches", "indels",
                        "unknowns", "Target", "CDS_Target", "locus", "Parent",
                        "partialness", "InFrameStop"
                    )) %>%
                    dplyr::arrange(seqnames, start) %>%
                    dplyr::rename(
                        "feature" = "type",
                        "id" = "ID",
                        "gene_id" = "geneIDs",
                        "transcript_id" = "transcripts"
                    ) %>%
                    dplyr::mutate(
                        source = paste0(
                            "sacCer3; trinity; ",
                            source,
                            ", gffread; Greenlaw et al."
                        ),
                        score = ".",
                        frame = ".",
                        gene_id = process_list_column(gene_id),
                        transcript_id = process_list_column(transcript_id)
                    ) %>%
                    dplyr::relocate(c(source, feature), .after = seqnames) %>%
                    dplyr::relocate(c(score, strand, frame), .after = end) %>%
                    dplyr::mutate(
                        attribute = paste(
                            paste0("id \"", id, "\""),
                            paste0("gene_id \"", gene_id, "\""),
                            paste0("transcript_id \"", transcript_id, "\""),
                            paste0("locus_source \"GffRead v0.12.7\""),
                            paste0("gene_source \"Trinity v2.14\""),
                            paste0("transcript_source \"Trinity v2.14\""),
                            sep = "; "
                        )
                    ) %>%
                    dplyr::select(-c(id, gene_id, transcript_id))
            } else if(base::isTRUE(unlist(stringr::str_split(k, "_"))[3] == "mRNA")) {
                if(base::isTRUE(unlist(stringr::str_split(k, "_"))[2] == "G")) {
                    gtf <- out_G_gtf_mRNA[[i]][[j]]
                }
                if(base::isTRUE(unlist(stringr::str_split(k, "_"))[2] == "Q")) {
                    gtf <- out_Q_gtf_mRNA[[i]][[j]]
                }
                gtf <- gtf %>%
                    tibble::as_tibble() %>%
                    dplyr::select(-c(
                        "phase", "width", "geneIDs", "transcripts", "Parent"
                    )) %>%
                    dplyr::arrange(seqnames, start) %>%
                    dplyr::rename(
                        "feature" = "type",
                        "id" = "ID",
                        "gene_id" = "geneID"
                    ) %>%
                    dplyr::mutate(
                        source = paste0(source, "; trinity; gffread; Greenlaw et al."),
                        score = ".",
                        frame = "."
                    ) %>%
                    dplyr::relocate(c(source, feature), .after = seqnames) %>%
                    dplyr::relocate(c(score, strand, frame), .after = end) %>%
                    dplyr::mutate(
                        attribute = paste(
                            paste0("id \"", id, "\""),
                            paste0("name \"", Name, "\""),
                            paste0("gene_id \"", gene_id, "\""),
                            paste0("locus_id \"", locus, "\""),
                            paste0("mRNA_source \"Trinity v2.14\""),
                            paste0("gene_source \"Trinity v2.14\""),
                            paste0("locus_source \"GffRead v0.12.7\""),
                            paste0("coverage \"", coverage, "\""),
                            paste0("identity \"", identity, "\""),
                            paste0("matches \"", matches, "\""),
                            paste0("mismatches \"", mismatches, "\""),
                            paste0("indels \"", indels, "\""),
                            paste0("unknowns \"", unknowns, "\""),
                            paste0("target \"", Target, "\""),
                            paste0("CDS_target \"", CDS_Target, "\""),
                            paste0("partialness \"", partialness, "\""),
                            paste0("in_frame_stop \"", InFrameStop, "\""),
                            sep = "; "
                        )
                    ) %>%
                    dplyr::select(-c(
                        id, Name, gene_id, locus, coverage, identity, matches,
                        mismatches, indels, unknowns, Target, CDS_Target,
                        partialness, InFrameStop
                    ))
            } else if(base::isTRUE(unlist(stringr::str_split(k, "_"))[3] == "exon")) {
                if(base::isTRUE(unlist(stringr::str_split(k, "_"))[2] == "G")) {
                    gtf <- out_G_gtf_exon[[i]][[j]]
                }
                if(base::isTRUE(unlist(stringr::str_split(k, "_"))[2] == "Q")) {
                    gtf <- out_Q_gtf_exon[[i]][[j]]
                }
                gtf <- gtf %>%
                    tibble::as_tibble() %>%
                    dplyr::select(-c(
                        "width", "phase", "ID", "geneIDs", "transcripts",
                        "geneID", "Name", "coverage", "identity", "matches",
                        "mismatches", "indels", "unknowns", "CDS_Target",
                        "partialness", "locus", "InFrameStop"
                    )) %>%
                    dplyr::arrange(seqnames, start) %>%
                    dplyr::rename("id" = "Parent", "feature" = "type") %>%
                    dplyr::mutate(
                        source = paste0(source, "; trinity; gffread; Greenlaw et al."),
                        frame = ".",
                        id = process_list_column(id)
                    ) %>%
                    dplyr::relocate(c("source", "feature"), .after = "seqnames") %>%
                    dplyr::relocate("frame", .after = "score") %>%
                    dplyr::mutate(
                        attribute = paste(
                            paste0("id \"", id, "\""),
                            paste0("source_id \"Trinity v2.14\""),
                            paste0("source_additional \"GffRead v0.12.7\""),
                            paste0("target \"", Target, "\""),
                            sep = "; "
                        )
                    ) %>%
                    dplyr::select(-c("Target", "id"))
            } else if(base::isTRUE(unlist(stringr::str_split(k, "_"))[3] == "CDS")) {
                if(base::isTRUE(unlist(stringr::str_split(k, "_"))[2] == "G")) {
                    gtf <- out_G_gtf_CDS[[i]][[j]]
                }
                if(base::isTRUE(unlist(stringr::str_split(k, "_"))[2] == "Q")) {
                    gtf <- out_Q_gtf_CDS[[i]][[j]]
                }
                gtf <- gtf %>%
                    tibble::as_tibble() %>%
                    dplyr::select(-c(
                        "width", "ID", "geneIDs", "transcripts", "geneID",
                        "Name", "coverage", "identity", "matches",
                        "mismatches", "indels", "unknowns", "CDS_Target",
                        "partialness", "locus", "InFrameStop"
                    )) %>%
                    dplyr::arrange(seqnames, start) %>%
                    dplyr::rename("id" = "Parent", "feature" = "type") %>%
                    dplyr::mutate(
                        source = paste0(
                            source,
                            "; trinity; gffread; Greenlaw et al."
                        ),
                        score = ".",
                        frame = ".",
                        id = process_list_column(id)
                    ) %>%
                    dplyr::relocate(c("source", "feature"), .after = "seqnames") %>%
                    dplyr::relocate(c("strand", "frame"), .after = "score") %>%
                    dplyr::mutate(
                        attribute = paste(
                            paste0("id \"", id, "\""),
                            paste0("source_id \"Trinity v2.14\""),
                            paste0("source_additional \"GffRead v0.12.7\""),
                            paste0("phase \"", phase, "\""),
                            paste0("target \"", Target, "\""),
                            sep = "; "
                        )
                    ) %>%
                    dplyr::select(-c("phase", "Target", "id"))
            } else if(base::isTRUE(unlist(stringr::str_split(k, "_"))[3] == "introns")) {
                if(base::isTRUE(unlist(stringr::str_split(k, "_"))[2] == "G")) {
                    gtf <- out_G_gtf_introns_filtered[[i]][[j]]
                }
                if(base::isTRUE(unlist(stringr::str_split(k, "_"))[2] == "Q")) {
                    gtf <- out_Q_gtf_introns_filtered[[i]][[j]]
                }
                gtf <- gtf %>%
                    tibble::as_tibble() %>%
                    dplyr::select(-c("width", "phase", "Parent")) %>%
                    dplyr::arrange(seqnames, start) %>%
                    dplyr::rename("feature" = "type", "id" = "ID") %>%
                    dplyr::mutate(
                        source = paste0(source, "; trinity; gffread; Greenlaw et al."),
                        score = ".",
                        frame = "."
                    ) %>%
                    dplyr::relocate(c(source, feature), .after = seqnames) %>%
                    dplyr::relocate(c(score, strand, frame), .after = end) %>%
                    dplyr::mutate(
                        attribute = paste(
                            paste0("id \"", id, "\""),
                            paste0("source_id \"Trinity v2.14\""),
                            paste0("source_additional \"GffRead v0.12.7\""),
                            sep = "; "
                        )
                    ) %>%
                    dplyr::select(-id)
            }
            
            write_gtf(gtf, outfile)
        }
    }
}

# rm(i, j, k, outfile, gtf)
```
<br />

```{r, eval=FALSE}
#  Write out gtfs using rtracklayer::export() ---------------------------------
#DONOTRUN #PREVIOUS
for(i in names(out_G_gtf_locus)) {
    for(j in c(paste0(
        "below_",
        c(
            "05", "10", "15", "20", "25", "30", "35", "40", "45", "50",
            "55", "60", "65", "70", "75", "80", "85", "90", "95"
        )
    ))) {
        for(k in c(
            path_G_locus, path_G_mRNA, path_G_exon, path_G_CDS,
            path_G_introns_filtered,
            path_Q_locus, path_Q_mRNA, path_Q_exon, path_Q_CDS,
            path_Q_introns_filtered
        )) {
            # i <- "k2"
            # j <- "below_25"
            # k <- path_Q_locus
            if(base::isTRUE(k == path_G_locus)) {
                outfile <- paste0(
                    path_G_locus, "/G1_mkc-", gsub("k", "", i), "_",
                    gsub("below_", "gte-pctl-", j), ".gtf"
                )
                
                cat("Saving G1", i, j, "in", k, "\n")
                cat(paste0(outfile, "\n"))
                cat("\n\n")
                
                out_G_conv <- GenomicRanges::makeGRangesFromDataFrame(
                    out_G_gtf_locus[[i]][[j]], keep.extra.columns = TRUE
                )

                rtracklayer::export(
                    out_G_gtf_locus[[i]][[j]],
                    paste0(
                        path_G_locus, "/G1_mkc-", gsub("k", "", i), "_",
                        gsub("below_", "gte-pctl-", j), ".gtf"
                    )
                )
            } else if(base::isTRUE(k == path_G_mRNA)) {
                outfile <- paste0(
                    path_G_mRNA, "/G1_mkc-", gsub("k", "", i), "_",
                    gsub("below_", "gte-pctl-", j), ".gtf"
                )
                
                cat("Saving G1", i, j, "in", k, "\n")
                cat(paste0(outfile, "\n"))
                cat("\n\n")
                
                out_G_conv <- GenomicRanges::makeGRangesFromDataFrame(
                    out_G_gtf_mRNA[[i]][[j]], keep.extra.columns = TRUE
                )

                rtracklayer::export(
                    out_G_gtf_mRNA[[i]][[j]],
                    paste0(
                        path_G_mRNA, "/G1_mkc-", gsub("k", "", i), "_",
                        gsub("below_", "gte-pctl-", j), ".gtf"
                    )
                )
            } else if(base::isTRUE(k == path_G_exon)) {
                outfile <- paste0(
                    path_G_exon, "/G1_mkc-", gsub("k", "", i), "_",
                    gsub("below_", "gte-pctl-", j), ".gtf"
                )
                
                cat("Saving G1", i, j, "in", k, "\n")
                cat(paste0(outfile, "\n"))
                cat("\n\n")
                
                out_G_conv <- GenomicRanges::makeGRangesFromDataFrame(
                    out_G_gtf_exon[[i]][[j]], keep.extra.columns = TRUE
                )

                rtracklayer::export(
                    out_G_gtf_exon[[i]][[j]],
                    paste0(
                        path_G_exon, "/G1_mkc-", gsub("k", "", i), "_",
                        gsub("below_", "gte-pctl-", j), ".gtf"
                    )
                )
            } else if(base::isTRUE(k == path_G_CDS)) {
                outfile <- paste0(
                    path_G_CDS, "/G1_mkc-", gsub("k", "", i), "_",
                    gsub("below_", "gte-pctl-", j), ".gtf"
                )
                
                cat("Saving G1", i, j, "in", k, "\n")
                cat(paste0(outfile, "\n"))
                cat("\n\n")
                
                out_G_conv <- GenomicRanges::makeGRangesFromDataFrame(
                    out_G_gtf_CDS[[i]][[j]], keep.extra.columns = TRUE
                )

                rtracklayer::export(
                    out_G_gtf_CDS[[i]][[j]],
                    paste0(
                        path_G_CDS, "/G1_mkc-", gsub("k", "", i), "_",
                        gsub("below_", "gte-pctl-", j), ".gtf"
                    )
                )
            } else if(base::isTRUE(k == path_G_introns_filtered)) {
                outfile <- paste0(
                    path_G_introns_filtered, "/G1_mkc-", gsub("k", "", i), "_",
                    gsub("below_", "gte-pctl-", j), ".gtf"
                )
                
                cat("Saving G1", i, j, "in", k, "\n")
                cat(paste0(outfile, "\n"))
                cat("\n\n")
                
                out_G_conv <- GenomicRanges::makeGRangesFromDataFrame(
                    out_G_gtf_introns_filtered[[i]][[j]], keep.extra.columns = TRUE
                )

                rtracklayer::export(
                    out_G_gtf_introns_filtered[[i]][[j]],
                    paste0(
                        path_G_introns_filtered, "/G1_mkc-", gsub("k", "", i), "_",
                        gsub("below_", "gte-pctl-", j), ".gtf"
                    )
                )
            } else if(base::isTRUE(k == path_Q_locus)) {
                outfile <- paste0(
                    path_Q_locus, "/Q_mkc-", gsub("k", "", i), "_",
                    gsub("below_", "gte-pctl-", j), ".gtf"
                )
                
                cat("Saving Q", i, j, "in", k, "\n")
                cat(paste0(outfile, "\n"))
                cat("\n\n")
                
                out_Q_conv <- GenomicRanges::makeGRangesFromDataFrame(
                    out_Q_gtf_locus[[i]][[j]], keep.extra.columns = TRUE
                )

                rtracklayer::export(
                    out_Q_gtf_locus[[i]][[j]],
                    paste0(
                        path_Q_locus, "/Q_mkc-", gsub("k", "", i), "_",
                        gsub("below_", "gte-pctl-", j), ".gtf"
                    )
                )
            } else if(base::isTRUE(k == path_Q_mRNA)) {
                outfile <- paste0(
                    path_Q_mRNA, "/Q_mkc-", gsub("k", "", i), "_",
                    gsub("below_", "gte-pctl-", j), ".gtf"
                )
                
                cat("Saving Q", i, j, "in", k, "\n")
                cat(paste0(outfile, "\n"))
                cat("\n\n")
                
                out_Q_conv <- GenomicRanges::makeGRangesFromDataFrame(
                    out_Q_gtf_mRNA[[i]][[j]], keep.extra.columns = TRUE
                )

                rtracklayer::export(
                    out_Q_gtf_mRNA[[i]][[j]],
                    paste0(
                        path_Q_mRNA, "/Q_mkc-", gsub("k", "", i), "_",
                        gsub("below_", "gte-pctl-", j), ".gtf"
                    )
                )
            } else if(base::isTRUE(k == path_Q_exon)) {
                outfile <- paste0(
                    path_Q_exon, "/Q_mkc-", gsub("k", "", i), "_",
                    gsub("below_", "gte-pctl-", j), ".gtf"
                )
                
                cat("Saving Q", i, j, "in", k, "\n")
                cat(paste0(outfile, "\n"))
                cat("\n\n")
                
                out_Q_conv <- GenomicRanges::makeGRangesFromDataFrame(
                    out_Q_gtf_exon[[i]][[j]], keep.extra.columns = TRUE
                )

                rtracklayer::export(
                    out_Q_gtf_exon[[i]][[j]],
                    paste0(
                        path_Q_exon, "/Q_mkc-", gsub("k", "", i), "_",
                        gsub("below_", "gte-pctl-", j), ".gtf"
                    )
                )
            } else if(base::isTRUE(k == path_Q_CDS)) {
                outfile <- paste0(
                    path_Q_CDS, "/Q_mkc-", gsub("k", "", i), "_",
                    gsub("below_", "gte-pctl-", j), ".gtf"
                )
                
                cat("Saving Q", i, j, "in", k, "\n")
                cat(paste0(outfile, "\n"))
                cat("\n\n")
                
                out_Q_conv <- GenomicRanges::makeGRangesFromDataFrame(
                    out_Q_gtf_CDS[[i]][[j]], keep.extra.columns = TRUE
                )

                rtracklayer::export(
                    out_Q_gtf_CDS[[i]][[j]],
                    paste0(
                        path_Q_CDS, "/Q_mkc-", gsub("k", "", i), "_",
                        gsub("below_", "gte-pctl-", j), ".gtf"
                    )
                )
            } else if(base::isTRUE(k == path_Q_introns_filtered)) {
                outfile <- paste0(
                    path_Q_introns_filtered, "/Q_mkc-", gsub("k", "", i), "_",
                    gsub("below_", "gte-pctl-", j), ".gtf"
                )
                
                cat("Saving Q", i, j, "in", k, "\n")
                cat(paste0(outfile, "\n"))
                cat("\n\n")
                
                out_Q_conv <- GenomicRanges::makeGRangesFromDataFrame(
                    out_Q_gtf_introns_filtered[[i]][[j]], keep.extra.columns = TRUE
                )

                rtracklayer::export(
                    out_Q_gtf_introns_filtered[[i]][[j]],
                    paste0(
                        path_Q_introns_filtered, "/Q_mkc-", gsub("k", "", i), "_",
                        gsub("below_", "gte-pctl-", j), ".gtf"
                    )
                )
            }
        }
    }
}
rm(i, j, k, outfile, out_G_conv, out_Q_conv)
```
</details>
<br />