GRCh38_STEP_1_annotation_forging.Rmd

---
title: "wind: wORKFLOW FOR PiRNAs AnD BEYONd"
subtitle: "Computational workflow for the creation of Gene transfer format file with small RNA sequences, GRCh38"
author: "Constantinos Yeles (Konstantinos Geles)"
date: "`r format(Sys.Date(), '%a %b %d %Y')`"
output:
  html_document:
    toc: yes
    toc_depth: 3
    theme: paper 
  pdf_document:
    toc: yes
    toc_depth: 3
editor_options:
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, eval = FALSE)
```

# Introduction
With the intent to annotate and quantify small RNA sequence data (and in 
particular piRNA) derived from Next-Generation Sequencing, we have developed wind. 
For the generation of annotation files and results, widely used tools of alignment, annotation, quantification and differential expression algorithms have been used. 
Although the workflow is focused particularly on piRNAs 
(as it is our main subject of research) with slight modifications can be applied
to all small RNA categories of interest.

To make it more versatile and reproducible, we adopted the _[containerization approach](https://www.docker.com/resources/what-container)_ as the software 
deployment is fast, efficient, and potentially bug-free. It can be used in 
various operating systems with only requirements the installation of the docker
engine and some minimum requirements of processing power and RAM to 
run the most memory demanding tools.

# Materials and Methods
The workflow has been primarily carried out on a Linux server, 
but it can be used easily on a Windows or Mac OS machine 
as long as changes would be done for particular functions/operations.

The workflow utilizes _[Bash](https://www.gnu.org/software/bash/)_ and 
_[R](https://www.r-project.org/)_ scripting for various operations.
For the application of the workflow, the following tools have been used:

*  _[Rstudio](https://rstudio.com/)_ for R scripting,

*  _[STAR](https://www.ncbi.nlm.nih.gov/pubmed/23104886)_ for alignment,

*  _[Samtools](https://www.htslib.org/)_ for various modifications and extraction
of reads from resulted aligned files,

*  _[FastQC](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)_ for quality control,

*  _[Cutadapt](https://journal.embnet.org/index.php/embnetjournal/article/view/200)_ for adapter trimming,

*  _[bedtools](https://bedtools.readthedocs.io/en/latest/index.html)_ for bam to bed manipulation, 

*  _[Salmon](https://www.nature.com/articles/nmeth.4197/)_ for transcript-level quantification,

*  _[featureCounts](https://academic.oup.com/bioinformatics/article/30/7/923/232889)_ for transcript-level quantification. 

Databases that have been used:

*  _[piRNABank](http://pirnabank.ibab.ac.in/)_ for piRNA sequences,

*  _[RNAcentral](https://rnacentral.org/)_ for smallRNA sequences.

# Workflow
## 1. Acquisition and Preprocessing of the small non-coding RNA (ncRNA) sequences
### i. Downloading the files for the generation of a Gene transfer format(gtf) 

To begin with, Human piRNA sequences are downloaded from piRNABank to enrich 
in piRNA sequences the gtf file, and we get the small RNA genome 
coordinates(bed files) from RNAcentral.

```{bash download the Databases}
# start 1st the docker container
docker run --rm -ti -v $(pwd):/home/my_data  congelos/sncrna_workflow

# all the files and folders for the workflow are created in the working directory 
#create the folder that will have the genome and smallRNA sequences information
mkdir -p my_data/human_data/GRCh38 my_data/human_data/piRNABank my_data/human_data/RNACentral

# download the piRNAbank sequences
wget http://pirnabank.ibab.ac.in/downloads/all/human_all.zip -O my_data/human_data/piRNABank/piRNA_human_all.zip

unzip -d my_data/human_data/piRNABank  my_data/human_data/piRNABank/piRNA_human_all.zip && rm my_data/human_data/piRNABank/piRNA_human_all.zip

# download the RNAcentral genomic coordinates
wget ftp://ftp.ebi.ac.uk/pub/databases/RNAcentral/releases/15.0/genome_coordinates/bed/homo_sapiens.GRCh38.bed.gz\
genome_coordinates/bed/homo_sapiens.GRCh38.bed.gz -O my_data/human_data/RNACentral/homo_sapiens.GRCh38.bed.gz

# get the GRCh38 fasta for STAR index 
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_34/GRCh38.primary_assembly.genome.fa.gz -O my_data/human_data/GRCh38/GRCh38.primary_assembly.genome.fa.gz

pigz -d my_data/human_data/GRCh38/GRCh38.primary_assembly.genome.fa.gz 

# get the GRCh38 annotation in order to exclude sequences of piRNA
# that are inside other sequences see 2.ix. chunk
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_34/gencode.v34.primary_assembly.annotation.gtf.gz -O my_data/human_data/GRCh38/gencode.v34.primary_assembly.annotation.gtf.gz
```
### ii. Preprocessing of the piRNAbank file

The fasta file from piRNAbank has U character instead of T in the sequences, 
so we can make an adjustment using [sed](https://www.gnu.org/software/sed/)
```{bash sed fasta}
sed 's/U/T/g' my_data/human_data/piRNABank/human_pir.txt > my_data/human_data/piRNABank/pirnaBank_human.fasta
# exit docker container
exit
```
### iii. Removing of the duplicated sequences from the piRNAbank file

We deploy the docker container with Rstudio server
```{bash Run docker bioc I}
docker run --rm -v $(pwd):/home/0 -p 8787:8787 -e PASSWORD=12345 \
-e USER=$UID congelos/rocker_tidyverse_plus_de_pckages
```
We prefer to work on Rstudio to perform everything on R otherwise R on 
bash could be used directly. We load the libraries that would assist us
in the creation of the piRNAbank fasta file.

```{r load the libraries I}
suppressPackageStartupMessages({
  library('stringr')
  library('plyranges')
  library('magrittr')
})
```

In the piRNABank fasta duplicated sequences exist and need to be removed
```{r remove duplicates}
# import the fasta  ----
pirnaB_hg19 <- Biostrings::readDNAStringSet("human_data/piRNABank/pirnaBank_human.fasta")
pirnaB_hg19 %>% length() ## should be: 667,944 sequences

# remove duplicated sequences---- 
pirnaB_hg19 <- pirnaB_hg19[!duplicated(pirnaB_hg19)]
pirnaB_hg19 %>% length() ##  23.439 sequences

# clean the names----
names(pirnaB_hg19) <- names(pirnaB_hg19) %>% 
  str_remove("\\|H.+") %>% 
  str_replace("\\|gb\\|","_")

# write the fasta ----
pirnaB_hg19 %>% 
  Biostrings::writeXStringSet("human_data/piRNABank/piRNAbank_hg19_rem_duplic.fa")
```
exit the docker container

```{bash exit bioc}
# exit docker container
exit
```
### iv. Align piRNA sequences to human genome

We proceed with the alignment of piRNA sequences to the human genome 
utilizing STAR aligner and export them in fasta format

```{bash STAR_SAM}
docker run --rm -ti -v $(pwd):/home/my_data  congelos/sncrna_workflow

# create index
STAR --runMode genomeGenerate --genomeDir my_data/human_data/GRCh38_v34 --genomeFastaFiles my_data/human_data/GRCh38/GRCh38.primary_assembly.genome.fa --runThreadN 8 && cp Log.out my_data/human_data/GRCh38_v34/Log.out
 
# align the piRNABank sequences
STAR --genomeDir my_data/human_data/GRCh38_v34/ --genomeLoad LoadAndKeep \
--readFilesIn "my_data/human_data/piRNABank/piRNAbank_hg19_rem_duplic.fa"  \
--runThreadN 4 --alignIntronMax 0 --outSAMattributes NH HI NM MD \
--outFilterMultimapNmax 100 --outReadsUnmapped Fastx --outFilterMismatchNmax 0 \
--outFilterMatchNmin 16  --outSAMtype BAM SortedByCoordinate \
--limitBAMsortRAM 30000000000 \
--outFileNamePrefix "my_data/human_data/piRNABank/aligned/piBnk_GRCh38_v34_"

# BAM to fasta format
samtools fasta -F 4 -@ 4 \
my_data/human_data/piRNABank/aligned/piBnk_GRCh38_v34_Aligned.sortedByCoord.out.bam > \
my_data/human_data/piRNABank/piBnk_GRCh38_v34_fin.fasta

# BAM to bed format
bedtools bamtobed < my_data/human_data/piRNABank/aligned/piBnk_GRCh38_v34_Aligned.sortedByCoord.out.bam > my_data/human_data/piRNABank/piBnk_GRCh38_v34_fin.bed

exit
```
## 2. Join of piRNABank sequences and RNAcentral ncRNA sequences 
### Run docker

```{bash Run docker bioc II}
docker run --rm -v $(pwd):/home/0 -p 8787:8787 -e PASSWORD=12345 \
-e USER=$UID congelos/rocker_tidyverse_plus_de_pckages
```

In order to minimize issues with the paths of  folders and files
we use wherever possible the package [here](https://github.com/jennybc/here_here#readme)

### i. Load libraries

```{r load the libraries II}
suppressPackageStartupMessages({
  library('tidyverse') 
  library('data.table')
  library('plyranges')
  library("BSgenome.Hsapiens.UCSC.hg38")
  library("here")
})
```
### ii. RNAcentral. import RNAcentral file 

```{r import RNAcentral}
sRNA <- here("human_data", "RNACentral", "homo_sapiens.GRCh38.bed.gz") %>% 
  read_bed() %>% 
  select("sRNA_id" = name, "gene_type" = NA.1, "source" = NA.2) %>% 
  mutate(type = "exon")

sInfo <- Seqinfo(genome="hg38")
seqlevels(sInfo) <- seqlevels(sRNA)
seqinfo(sRNA) <- sInfo
```
#### Exploring the lengths of the RNAcentral sequences

```{r summaries of length RNAcentral seq}
sRNA %>% 
  as_tibble() %>% 
  group_by(gene_type, sRNA_id) %>% 
  summarise(count_less = sum(width < 100), count_more = sum(width >= 100)) %>%
  group_by(gene_type) %>% 
  summarise(longer_or_equal_than_100_seqs = sum(count_more>0), 
            shorter_than_100_seqs= sum(count_less>0),
            longer_or_equal_than_100_GRs = sum(count_more), 
            shorter_than_100_GRs= sum(count_less))
```
### iii. RNAcentral. filtering for sequences smaller than 100 nucleotides

```{r filter 100bp}
tr_sRNA <- sRNA %>% 
  as_tibble() %>% # [602,197] genomic ranges (GR) / [446,265] sRNA_ids
  filter(width < 100) %>% # [162,958] GR / [44,556] sRNA_ids
  mutate(sRNA_id = str_remove(sRNA_id,"_9606")) %>% 
  as_granges() %>% 
  # keep info about the standard chromosomes
  keepStandardChromosomes(pruning.mode = "coarse") %>% # [162,958] -> [160,980] GRs / [44,556] -> [44,529] sRNA_ids
  # remove the duplicated GR entries from RNAcentral
  as_tibble() %>% 
  unite(col = "seq_s",seqnames:strand, sep = "_") %>% 
  distinct(seq_s, .keep_all = TRUE) %>% # remove dupl GRs [160,980] -> [153,043] GRs / [44,529] -> [41,496] sRNA_ids
  separate(col = seq_s,
           into = c("seqnames","start","end","width","strand"),
           sep = "_") %>% 
  mutate(start = as.numeric(start),
    end = as.numeric(end),
    width = as.numeric(width)) %>% 
  as_granges()

biotypes <- tr_sRNA %>% 
  as_tibble() %>% 
  select(sRNA_id,gene_type) %>% 
  distinct(sRNA_id, .keep_all = T)
```
### iv. RNAcentral. keep sequence information

```{r keep seq info}
transcripts_human <- Views(BSgenome.Hsapiens.UCSC.hg38, tr_sRNA)

# search for duplicated sequences ----
fasta_tr_hs <- DNAStringSet(transcripts_human)
names(fasta_tr_hs) <- mcols(transcripts_human)$sRNA_id
fasta_tr_hs <- fasta_tr_hs[sort(fasta_tr_hs@ranges@NAMES)]

fasta_tr_hs_tbl <- fasta_tr_hs %>% 
  as.character() %>% 
  enframe(name = "tr_hg38", value = "hg38") %>% 
  left_join(biotypes, by = c("tr_hg38" = "sRNA_id"))

fasta_tr_hs_tbl %>% 
  distinct(tr_hg38, hg38, gene_type, .keep_all = TRUE) %>% 
  filter(duplicated(hg38)) %>%
  count(hg38,sort = TRUE) # 27 duplicated sequences

### duplicates between sequences and genomic locations(GRs)----
# make a tibble with all GR, seq and ids
transcripts_GR <- transcripts_human %>% 
    as_granges() %>% 
    as_tibble() %>% 
    unite(col = "seq_RCent", seqnames:strand, sep = "_") 

# find unique pairs of sequences and GR
uniq_seq <- transcripts_GR  %>% 
  distinct(dna, .keep_all = TRUE) %>% 
  arrange(dna) %>% 
  mutate(sRNA_id2 = str_c(sRNA_id,"_GR_",seq_RCent)) %>% 
  select(dna,sRNA_id2)

transcripts_GR <- transcripts_GR %>%  # [153.043] GRs / [41,496] sRNA_id -> [43,575] sRNA_id2 
  left_join(uniq_seq)

```
### v. piRNABank. import the piRNA sequences aligned to genome 

```{r import piRNABank}
piRNAbank_hg38 <- here("human_data", "piRNABank", "piBnk_GRCh38_v34_fin.fasta") %>% 
  Biostrings::readDNAStringSet()

piRNAbank_hg38_tib <- piRNAbank_hg38 %>% # 23,120 sequences
  as.character() %>% 
  enframe(value = "seq_piBn") %>% 
  mutate(sRNA_type ="piRNA",
    bpairs_piR = str_length(seq_piBn)) %>% 
  arrange(desc(bpairs_piR))
```
### vi. piRNAbank. make Genomic Ranges and remove duplicates from hg38 

```{r import piRNABank GRs}
# piRNABank. import the Genomic Ranges and filter them -----
  piRNAbank_hg38_ranges <- here("human_data","piRNABank", "piBnk_GRCh38_v34_fin.bed") %>%
  read_bed() %>% 
  as_tibble() %>% 
  arrange(desc(width)) %>% # [46,552] GRs / 23,120 sequences
  filter(width < 100) %>% #  [46,552] -> [46,503] GRs / 23,120
  as_granges() %>% 
  keepStandardChromosomes(pruning.mode = "coarse") # [46,503] -> [45,818] GRs / 23,120 seq

transcripts_pi_hg38 <- Views(BSgenome.Hsapiens.UCSC.hg38, piRNAbank_hg38_ranges) %>%
  as_granges() %>%
  keepStandardChromosomes(pruning.mode = "coarse") %>% 
  as_tibble()

# we need to apply a second width filter at 69 
# as we know that piRNAs are ~32 base pairs
# but if we search the length of GRs there are some with width >= 69, 8 in particular
transcripts_pi_hg38 <- transcripts_pi_hg38 %>% filter(width < 69)  # [45,818] -> [45,810] GRs / 23,120 seq
transcripts_pi_hg38 %>% count(name) %>% nrow #> 23120 piRNAs from piRNABANK

# checking sequences of alignments with lower length
# here we explore the GRs that have one less at least base than the 
# piRNA sequence we actually have from the principal piRNABank fasta
sequen_pi_false <- transcripts_pi_hg38 %>% 
  as_tibble() %>% 
  left_join(piRNAbank_hg38_tib) %>% 
  arrange(desc(width)) %>% 
  mutate(sequences_true = (dna == seq_piBn)) %>% 
  filter(sequences_true == FALSE) %>% 
  unite(col = "seq_s",seqnames:strand, sep = "_") 

# piRNABank. removing duplicated GR ----
piRNAbank_hg38_ranges %>% 
  as_tibble() %>%
  unite(col = "seq_s", seqnames:strand, sep = "_") %>% 
  count(seq_s, sort = TRUE) %>% 
  filter(n > 1) %>% # 38 sequences have the same genomic range with at least another 
  .$seq_s %>% 
  map(~sequen_pi_false %>% 
  filter(seq_s == .x)) %>% 
  bind_rows()

transcripts_pi_hg38_clean <- transcripts_pi_hg38 %>% 
  as_tibble() %>% 
  left_join(piRNAbank_hg38_tib) %>% 
  arrange(desc(width)) %>% 
  mutate(sequences_true = (dna == seq_piBn)) %>% 
  filter(sequences_true == TRUE) %>% # [45,810] -> [44,557] GRs/ [23,120] -> [23,116] sequences
  select(-c(score, seq_piBn, bpairs_piR, sequences_true)) %>% 
  unite(col = "seq_piBNK", seqnames:strand, sep = "_")

transcripts_pi_hg38_clean %>% count(name) %>% nrow #>  23,116 piRNAs final piRNABANK
```
### vii. RNAcentral. + piRNABank. make annotation tibble
create a tibble with the information of RNAcentral, piRNAbank sequences and IDs

```{r annotation tibble}
hg38_piBAnk_RCent <- transcripts_GR %>% 
  left_join(piRNAbank_hg38_tib, by = c("dna" = "seq_piBn")) 

# check gene_types
hg38_piBAnk_RCent %>% 
  filter(is.na(name)) %>% 
  count(gene_type) 

hg38_piBAnk_RCent %>% 
  filter(!is.na(name)) %>% 
  count(gene_type,sRNA_type) 

concated_hg38_piBAnk <- hg38_piBAnk_RCent  %>% 
  mutate(
    seq_id = case_when(
      is.na(gene_type) ~ name,
      gene_type != "piRNA" ~ sRNA_id2,
      is.na(sRNA_type) ~ sRNA_id2,
      TRUE ~ name
    )
  ) %>% 
  select(seq_RCent, seq_id, gene_type,sRNA_id,sRNA_id2,  everything())

# sanity checks 1 ----
## check for the NA values, should be only true
(concated_hg38_piBAnk %>% 
    filter(is.na(name)) %>% .$sRNA_id2 == 
    concated_hg38_piBAnk %>% 
    filter(is.na(name)) %>% .$seq_id
  ) %>% table

## check for the miRNA values, should be only true
 (concated_hg38_piBAnk %>% filter(gene_type == "miRNA") %>% .$sRNA_id2 == 
     concated_hg38_piBAnk %>% filter(gene_type == "miRNA") %>% .$seq_id
   ) %>% table

## function for all gene_types
fun_unm <- function(x){
  (concated_hg38_piBAnk %>% 
      filter(gene_type == x) %>% 
      .$sRNA_id2 == 
      concated_hg38_piBAnk %>% 
      filter(gene_type == x) %>% 
      .$seq_id
    ) %>% table
}
## check for all gene_types, should be only true except piRNAs
concated_hg38_piBAnk %>% 
  count(gene_type) %>% 
  .$gene_type %>% set_names(.) %>%  
  map(~fun_unm(.x)) %>% 
  bind_rows(.id = "ID")

## check for the piRNA values
is.na(concated_hg38_piBAnk$seq_id) %>% table

## check for duplicates
concated_hg38_piBAnk %>% 
  filter(duplicated(seq_id)) %>% 
  arrange(name)

concated_hg38_piBAnk %>% 
  filter(duplicated(name),!is.na(name)) %>% 
  arrange(name)

concated_hg38_piBAnk %>% 
  filter(duplicated(sRNA_id2),!is.na(sRNA_id2)) %>% 
  arrange(name)

dupl_seqs <- concated_hg38_piBAnk %>% 
  filter(duplicated(dna)) %>% 
  arrange(name)

fasta_tr_hs_tbl %>% 
  filter(hg38 %in% dupl_seqs$dna)

# sanity checks 2 ----
concated_hg38_piBAnk %>% 
  filter(!duplicated(dna)) %>% 
  select(sRNA_id,seq_id) %>% 
  arrange(sRNA_id)

tr_test <- transcripts_human %>% as_granges()

concated_hg38_piBAnk %>% filter(sRNA_id == "URS0000000096")
transcripts_pi_hg38_clean %>%  filter(name == "hsa_piR_009796_DQ583192")

tr_test %>% filter(sRNA_id == "URS0000000096")

## check for the duplicated sequences
concated_hg38_piBAnk %>% 
  filter(duplicated(dna)) %>% 
  select(seq_RCent,seq_id)

concated_hg38_piBAnk %>% filter(sRNA_id == "URS00001B5714")
transcripts_pi_hg38_clean %>%  filter(name == "hsa_piR_011289_DQ585240")
tr_test %>% filter(sRNA_id == "URS00001B5714")
```
### viii. RNAcentral + piRNABank. generation of GRs

```{r GRs generation}
concated_hg38_piBAnk  # dataframe with combined sequences piRNAbank+RNAcentral
transcripts_pi_hg38_clean # dataframe that has filtered alignments from piRNAbank 

# here we change the small-RNA category of a piRNA from piRNABank if it is found 
# in RNACentral with another small-RNA category
c_piBNK_RCent <- concated_hg38_piBAnk %>%
  full_join(transcripts_pi_hg38_clean, 
            by = c("dna", "name", "sRNA_type", "seq_RCent" =  "seq_piBNK")) %>% 
  select(seq_RCent, sRNA_id, name, seq_id, gene_type, sRNA_type, everything()) %>%
  mutate(
    source = 
      case_when(
        is.na(source) ~ "piRNA_BANK",
        !is.na(sRNA_type) ~ str_c("piRNA_BANK,",source),
        TRUE ~ source),
    gene_type = 
      case_when(
        is.na(gene_type) ~ sRNA_type,
        TRUE ~ gene_type),
    seq_id =
      case_when(
        is.na(seq_id) ~ name,
        TRUE ~ seq_id),
    type =
      case_when(
        is.na(type) ~ "exon",
        TRUE ~ type)
  ) 
      

names(c_piBNK_RCent) 
c_piBNK_RCent %>% count(sRNA_id, sort = T)

c_piBNK_RCent %>% count(name, sort = T)

c_piBNK_RCent %>% count(seq_id, sort = T)

c_piBNK_RCent %>% count(gene_type, sort = T)# GR per gene_type

c_piBNK_RCent %>% count(source, sort = T)

c_piBNK_RCent %>% count(type, sort = T)

c_piBNK_RCent %>% count(dna, sort = T) 

c_piBNK_RCent %>% count(sRNA_id2, sort = T)

c_piBNK_RCent %>% count(bpairs_piR, sort = T)

# combined Genomic ranges -----
c_piBNK_RCent_GR <- c_piBNK_RCent %>% 
  select(-name, -sRNA_type, -bpairs_piR) %>% 
  dplyr::rename(  "seq_RNA" = dna ) %>% 
  separate(col = seq_RCent,into = c("seqnames",
    "start","end","width","strand"),sep = "_") %>% 
  mutate(start = as.numeric(start),
    end = as.numeric(end),
    width = as.numeric(width)) %>% 
  as_granges
```
### ix.  Filtering piRNAs sequences inside other Genes

Studies of Tosar et. al[1](https://pubmed.ncbi.nlm.nih.gov/33376191/) [2](https://pubmed.ncbi.nlm.nih.gov/30271890/) have demonstrated that some piRNAs
sequences can been found inside other small non-coding RNAs or genes.
These "piRNA" most probably are mRNA or ncRNA fragments that have been misstyped 
as such.
For this reason we will exclude the Genomic ranges of piRNAs that are inside 
other protein coding genes. Regarding piRNAs inside  regions of small noncoding
RNAs we can follow up after the final downstream analysis has been done.

Always consider to experimentally validate the sequences of interest
with wet-lab techniques (after the smallRNA sequencing) that could prove that
these sequences are actual PIWI-interacting small non-coding RNAs.

We first evaluate which Biotypes should use from
gencode annotation.
A cut-off of at least 1 base of overlap with the piRNA sequence in order
to be as strict as possible is used

We filter the GRs that are found inside protein coding of CDS and exon tags

```{r filtering protein coding GRs}
gene_annot <- here("human_data", "GRCh38", "gencode.v34.primary_assembly.annotation.gtf.gz") %>% 
  read_gff2()

# how many GRs for each type are in the gencode annotation:
gene_annot %>% 
  as_tibble %>% 
  count(type, sort = T) # exon:1379143, CDS:764194
# how many GRs for each type per gene_type are in the gencode annotation:
gene_annot %>% 
  as_tibble %>% 
  count(type, gene_type, sort = T) # exon:protein_coding= 1149694,
                                   # CDS:protein_coding=762887

#filter gene_annot for only exon and CDS
gene_annot_exon_CDS <- gene_annot %>% filter(type %in% c("exon", "CDS"))

#filter exon and CDS for protein_coding and immunoglobulin genes
gene_annot_exon_CDS_prot_cod <- gene_annot_exon_CDS %>% 
  filter(type %in% c("exon", "CDS"), 
         str_detect(gene_type, "protein_coding|IG_|TR_"))

# we will consider overlapping piRNAs with exon and CDS of protein_coding first
c_piRNAs_GR_prot_cod <- c_piBNK_RCent_GR %>% # GR = 155.143
  filter(gene_type == "piRNA") %>% # GR = 139.196
  find_overlaps_directed(gene_annot_exon_CDS_prot_cod, maxgap = -1L, 
                         minoverlap = 1L, suffix = c("_piRNAs", "_genes")) %>% 
  select(seq_id, starts_with(c("gene_","typ"))) %>% 
  as_tibble() %>% 
  unite(col = "seq_piBNK", seqnames:strand, sep = "_") %>% 
  distinct(seq_piBNK, .keep_all = TRUE) %>% 
  write_tsv(file.path("human_data", "excluded_sncRNA_piRNBnk_RNACent_GRCh38_v34.txt"))# GR = 5.594 to be excluded

# exclude them from combined ranges
c_piBNK_RCent_GR_filtered <- c_piBNK_RCent_GR %>% 
  as_tibble() %>% 
  unite(col = "seq_piBNK", seqnames:strand, sep = "_") %>% 
  filter(!seq_piBNK %in% c_piRNAs_GR_prot_cod$seq_piBNK) %>%  # 155.143 -> 149.549
  separate(col = seq_piBNK,
           into = c("seqnames", "start","end","width","strand"), 
           sep = "_", convert = TRUE) %>% 
  as_granges()
```
Then we will add some information regarding piRNAs that overlap with 
small non coding RNAs and other pseudogenes

```{r piRNA overlapping GENCODE_pseudogenes}
# piRNA overlap with sequences of the 25 types(containing pseudogenes and smallRNAs) from GENCODE:
gene_annot_exon_CDS %>% 
  filter(!gene_type %in% gene_annot_exon_CDS_prot_cod$gene_type) %>% 
  as_tibble %>% 
  count(gene_type, sort = T) %>%
  .$gene_type

# keep only the sequences of the 25 types
gene_annot_pseud_n_others <- gene_annot_exon_CDS %>% 
  filter(!gene_type %in% gene_annot_exon_CDS_prot_cod$gene_type) 

# make a dataframe with the overlapping GRs of piRNAs
piRNAs_GR_filt <- c_piBNK_RCent_GR_filtered %>% 
  filter(gene_type == "piRNA") %>% # GR = 133.602
  find_overlaps_directed(gene_annot_pseud_n_others, maxgap = -1L, 
                         minoverlap = 1L, suffix = c("_piRNAs", "_genes")) %>% 
  as_tibble() %>% 
  unite(col = "seq_piBNK", seqnames:strand, sep = "_") %>% 
  distinct(seq_piBNK, .keep_all = TRUE) %>% # GR = 12.897 to be recategorized
  select(seq_piBNK, gene_name, gene_id, transcript_id, gene_type_genes,tag) %>% 
  unite(col = "GENCODE_annot", gene_name:tag, sep = "_GNC_")
```
Using the information from GENECODE
we can add it to the final gtf as extra information

```{r miscellaneous piRNA merge}
# complete GRs with misc_piRNAs
complete_piBnk_RCent <- c_piBNK_RCent_GR_filtered %>% 
  as_tibble() %>% 
  unite(col = "seq_piBNK", seqnames:strand, sep = "_") %>% 
  left_join(piRNAs_GR_filt) %>% 
  mutate(across(.cols = GENCODE_annot,
                 .fns = ~if_else(is.na(.x),"no_overlap",.x))
           ) %>% 
  select(seq_piBNK, seq_id, GENCODE_annot, everything()) %>% 
  separate(col = seq_piBNK,
           into = c("seqnames", "start","end","width","strand"), 
           sep = "_", convert = TRUE) %>% 
  as_granges()
```
Last check for the sequences of all smallRNAs

```{r  objects to export}
# test the sequences-----
piRNAbank_rCentral_seqs <- Views(BSgenome.Hsapiens.UCSC.hg38, complete_piBnk_RCent)

piRNAbank_rCentral_seqs %>%
  as_granges() %>% 
  as_tibble() %>% 
  mutate( is_it_TR = (seq_RNA == dna)) %>% 
  filter(is_it_TR == FALSE)# should be 0

# final objects to export-----
piRNAbank_rCentral_fasta <- DNAStringSet(piRNAbank_rCentral_seqs)
names(piRNAbank_rCentral_fasta) <- mcols(piRNAbank_rCentral_seqs)$seq_id  

piRNAbank_rCentral_fasta <- piRNAbank_rCentral_fasta[!duplicated(piRNAbank_rCentral_fasta)]
```
## 3. Save the results to fasta and gtf format

```{r export and save res}
piRNAbank_rCentral_fasta %>% 
  Biostrings::writeXStringSet(file.path("human_data", "sncRNA_piRNBnk_RNACent_GRCh38_v34.fa"))

gtf_piB_RCentr <- complete_piBnk_RCent %>% 
  as_tibble() %>% 
  dplyr::rename("gene_id" = seq_id) %>% 
  as_granges()

sInfo <- Seqinfo(genome="hg38")
seqlevels(sInfo) <- seqlevels(gtf_piB_RCentr)
seqinfo(gtf_piB_RCentr)<- sInfo

gtf_piB_RCentr %>% 
  write_gff2(file.path("human_data", "sncRNA_piRNBnk_RNACent_GRCh38_v34.gtf"))
```

Until now, we have prepared the files for annotation and quantification (GTF, FASTA) 
of smallRNA sequencing samples of human transcriptome. 
Afterwards, the steps in the pre-processing of the samples, alignment, 
quantification and calculation of transcript abundances should be followed.

## 5. Make the directories for each dataset
We will create folders for each dataset we will analyse.

```{r}
data_sets <- c("Testis_COLO205", "TCGA_BRCA", "GSE68246")
data_sets %>% 
  file.path("Datasets_analysis",.) %>% 
  map(~dir.create(., recursive = T))
```
## 6. Optional* Make the indexes for STAR and salmon with spike-in sequences
### i. make the gtf merged file(spike-ins and smallRNAs)

```{r adding spikes-ins in genome and gtf}
#make the folder in which we will put the indexes
index_folders <- file.path("human_data","indexes",c("GRCh38_v34_spike_ins", "GRCh38_v34_public")) %>% 
  str_c(rep(c("_STAR","_salmon"),each=2)) 

index_folders %>% 
  map(~dir.create(., recursive = TRUE))

# load the spike in sequences and make it as Genomic ranges
spike_ins <- read_tsv("spike-ins.txt", col_names = c("names", "seq_RNA")) %>% 
  dplyr::rename("seqnames" = names) %>% 
  mutate(
  start = 1,
  end = str_length(seq_RNA),
  width = str_length(seq_RNA),
  strand = "+",
  gene_id = seqnames,
  type = "exon",
  source = "spike_in",
  gene_type = "spike_in"
)

# make the merged gtf with spike-in and smallRNA sequences
spike_ins_GR <- spike_ins %>% 
  as_granges()
genome(spike_ins_GR) <- "hg38"
 
piB_RCentr_spike_ins <- spike_ins_GR %>% 
  bind_ranges(gtf_piB_RCentr)

piB_RCentr_spike_ins %>% 
  write_gff2(file.path("human_data","sncRNA_spike_ins_piRNBnk_RNACent_GRCh38_v34.gtf"))

# make a fasta file of the sequences and then use fcat to concat the files on bash (way faster)
spikes_Fasta <- spike_ins$seq_RNA %>%
  DNAStringSet(start = rep(1,nrow(spike_ins)), end = str_length(spike_ins$seq_RNA)) %>% 
  setNames(spike_ins$seqnames) 

spikes_Fasta %>% 
  Biostrings::writeXStringSet(file.path("human_data","spike_ins.fasta"))
```
### ii. fcat bash use, make the fasta genome transcriptome
we will use fcat to make the fasta files, one of the genome and spike-ins and one of
the transcriptome and the spike-ins. Both files will be used to generate the
indexes of salmon and STAR

```{bash fcat}
# for now it works on linux only 
# cargo install fcat

# transriptome fasta
fcat my_data/human_data/spike_ins.fasta my_data/human_data/sncRNA_piRNBnk_RNACent_GRCh38_v34.fa  > my_data/human_data/transcriptome_sncRNA_spike_ins_piRNBnk_RNACent_GRCh38_v34.fasta

# genome fasta
fcat my_data/human_data/spike_ins.fasta my_data/human_data/GRCh38/GRCh38.primary_assembly.genome.fa > my_data/human_data/genome_sncRNA_spike_ins_piRNBnk_RNACent_GRCh38_v34.fasta
```
### iii. make a txt file for future use about the generation of histograms

```{r file for selecting piRNA for histograms}
# file for histograms of piRNAs
piB_RCentr_spike_ins %>% 
  as_tibble() %>%
  filter(gene_type %in% c("amb_piRNA", "misc_piRNA", "piRNA","spike_in")) %>% 
  plyranges::select(gene_id, seq_RNA, gene_type) %>%
  distinct(seq_RNA , .keep_all = TRUE) %>% 
  write_tsv(file.path("human_data","all_piRNAs_spike_ins_for_hists.txt"))
```
### iv. Indexes for STAR and Salmon with Spike-ins

```{bash spike-ins indexes}
docker run --rm -ti -v $(pwd):/home/my_data  congelos/sncrna_workflow

# make the STAR index
STAR --runMode genomeGenerate \
--genomeDir my_data/human_data/indexes/GRCh38_v34_spike_ins_STAR \
--genomeFastaFiles my_data/human_data/genome_sncRNA_spike_ins_piRNBnk_RNACent_GRCh38_v34.fasta --runThreadN 8 && cp Log.out my_data/human_data/indexes/GRCh38_v34_spike_ins_STAR/Log.out 

# create the STAR index without  spike-ins
STAR --runMode genomeGenerate \
--genomeDir my_data/human_data/indexes/GRCh38_v34_public_STAR \
--genomeFastaFiles my_data/human_data/GRCh38/GRCh38.primary_assembly.genome.fa --runThreadN 8 && cp Log.out my_data/human_data/indexes/GRCh38_v34_public_STAR/Log.out

# following the instructions for salmon decoy aware indexing
# https://combine-lab.github.io/alevin-tutorial/2019/selective-alignment/
grep "^>" my_data/human_data/GRCh38/GRCh38.primary_assembly.genome.fa |cut -d " " -f 1 > my_data/human_data/decoys_GRCh38.txt

sed -i.bak -e 's/>//g' my_data/human_data/decoys_GRCh38.txt

#concat 
fcat my_data/human_data/transcriptome_sncRNA_spike_ins_piRNBnk_RNACent_GRCh38_v34.fasta my_data/human_data/genome_sncRNA_spike_ins_piRNBnk_RNACent_GRCh38_v34.fasta > my_data/human_data/gentrome_sncRNA_spike_ins_piRNBnk_RNACent_GRCh38_v34.fasta

pigz --best -p 10 my_data/human_data/gentrome_sncRNA_spike_ins_piRNBnk_RNACent_GRCh38_v34.fasta 

# the same without the spike-ins
fcat my_data/human_data/sncRNA_piRNBnk_RNACent_GRCh38_v34.fa my_data/human_data/GRCh38/GRCh38.primary_assembly.genome.fa > my_data/human_data/gentrome_sncRNA_piRNBnk_RNACent_GRCh38_v34.fasta

pigz --best -p 10 my_data/human_data/gentrome_sncRNA_piRNBnk_RNACent_GRCh38_v34.fasta 

mkdir  my_data/human_data/GRCh38_v34_salmon
exit
# run the docker
docker run --rm -it -v $(pwd):/home/my_data combinelab/salmon
# create the index with spike-ins
## 
salmon index -t my_data/human_data/gentrome_sncRNA_spike_ins_piRNBnk_RNACent_GRCh38_v34.fasta.gz  \
-d my_data/human_data/decoys_GRCh38.txt \
-imy_data/human_data/indexes/GRCh38_v34_spike_ins_salmon \
-k 15 -p 14

# the same without the spike-ins
salmon index -t my_data/human_data/gentrome_sncRNA_piRNBnk_RNACent_GRCh38_v34.fasta.gz \
-d my_data/human_data/decoys_GRCh38.txt \
-i my_data/human_data/indexes/GRCh38_v34_public_salmon -k 15 -p 14

exit
```