GRCm38_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, GRCm38"
author: "Constantinos Yeles (Konstantinos Geles)"
date: "`r format(Sys.Date(), '%a %b %d %Y')`"
output:
  pdf_document:
    toc: yes
    toc_depth: 3
  html_document:
    toc: yes
    toc_depth: 3
    theme: paper 
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 have 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 have been done to 
appropriate 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) 

Mouse piRNA sequences were downloaded from piRNABank to enrich in piRNA sequences the gtf file, and small-RNA genome coordinates (bed files) from RNACentral have been acquired
```{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/mouse_data/GRCm38 my_data/mouse_data/piRNABank my_data/mouse_data/RNACentral

# download the piRNAbank sequences
wget http://pirnabank.ibab.ac.in/downloads/all/mouse_all.zip -O my_data/mouse_data/piRNABank/piRNA_mouse_all.zip

unzip -d my_data/mouse_data/piRNABank  my_data/mouse_data/piRNABank/piRNA_mouse_all.zip && rm my_data/mouse_data/piRNABank/piRNA_mouse_all.zip

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

# download the RNAcental id mappings
wget http://ftp.ebi.ac.uk/pub/databases/RNAcentral/releases/15.0/id_mapping/id_mapping.tsv.gz \
 -O my_data/mouse_data/RNACentral/mus_musculus.GRCm38.id_mapping.tsv.gz

# get the GRCm38 fasta for STAR index 
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M25/GRCm38.primary_assembly.genome.fa.gz \
 -O my_data/mouse_data/GRCm38/GRCm38.primary_assembly.genome.fa.gz

pigz -d my_data/mouse_data/GRCm38/GRCm38.primary_assembly.genome.fa.gz

# get the GRCm38 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_mouse/release_M25/gencode.vM25.primary_assembly.annotation.gtf.gz \
 -O my_data/mouse_data/GRCm38/gencode.vM25.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/mouse_data/piRNABank/mouse_pir.txt > my_data/mouse_data/piRNABank/pirnaBank_mouse.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')
})
```
### iv. Remove duplicated sequences

In the piRNABank fasta duplicated sequences exist and need to be removed
```{r remove duplicates}
pirnaB_mm8 <- Biostrings::readDNAStringSet(file.path("mouse_data", 
                                                     "piRNABank",
                                                     "pirnaBank_mouse.fasta"))
pirnaB_mm8 %>% length() ## >[1] 1399813

# remove duplicate sequences-----
pirnaB_mm8 <- pirnaB_mm8[!duplicated(pirnaB_mm8)]
pirnaB_mm8 %>% length() ## >[1] 39986

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

# write the fasta ----
pirnaB_mm8 %>% 
  Biostrings::writeXStringSet(file.path("mouse_data", "piRNABank" ,"pirnaB_mm8_removed_duplicates.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/mouse_data/GRCm38 --genomeFastaFiles my_data/mouse_data/GRCm38/GRCm38.primary_assembly.genome.fa --runThreadN 6

mkdir my_data/mouse_data/piRNABank_mouse_mm10

# align the piRNABank sequences
STAR --genomeDir my_data/mouse_data/GRCm38_v34  \
--readFilesIn "my_data/mouse_data/piRNABank/pirnaB_mm8_removed_duplicates.fa"  \
--runThreadN 8 --alignIntronMax 0 --outSAMattributes NH HI NM MD \
--outFilterMultimapNmax 100 --outReadsUnmapped Fastx --outFilterMismatchNmax 0 \
--outFilterMatchNmin 16 --outSAMtype BAM SortedByCoordinate \
--outFileNamePrefix "my_data/mouse_data/piRNABank/aligned/piBnk_GRCm38_v34_"

# BAM to fasta format
samtools fasta -F 4 -@ 4 \
my_data/mouse_data/piRNABank/aligned/piBnk_GRCm38_v34_Aligned.sortedByCoord.out.bam > my_data/mouse_data/piRNABank/piBnk_GRCm38_v34_fin.fasta


# BAM to bed format
bedtools bamtobed < my_data/mouse_data/piRNABank/aligned/piBnk_GRCm38_v34_Aligned.sortedByCoord.out.bam > my_data/mouse_data/piRNABank/piBnk_GRCm38_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.Mmusculus.UCSC.mm10")
  library("here")
})
```
### ii. RNAcentral. import RNAcentral file 

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

sInfo <- Seqinfo(genome="mm10")
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 bps

```{r filter 100bp}
tr_sRNA <- sRNA %>% 
  as_tibble() %>% # [1,065,469] genomic ranges(GR) / [180,780] sRNA_ids
  filter(width < 100) %>% # [962,711] GRs / [94,653] sRNA_ids
  mutate(sRNA_id = str_remove(sRNA_id,"_10090")) %>% 
  as_granges() %>% 
  # keep info about the standard chromosomes
  keepStandardChromosomes(pruning.mode = "coarse") %>% # [962,711] -> [947,408] GRs / [94,653] -> [94,644] sRNA_ids
  # remove the duplicated entries from RNAcentral
  as_tibble() %>% 
  unite(col = "seq_s",seqnames:strand, sep = "_") %>% 
  distinct(seq_s, .keep_all = TRUE) %>% # remove dupl GRs [947,408] -> [929,168] GRs / [94,644] -> [90,056] 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_mouse <- Views(BSgenome.Mmusculus.UCSC.mm10, tr_sRNA)

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

fasta_tr_ms_tbl <- fasta_tr_ms  %>% 
  as.character() %>% 
  enframe(name = "tr_mm10" ,value = "mm10") %>% 
  left_join(biotypes, by = c("tr_mm10" = "sRNA_id"))

fasta_tr_ms_tbl %>% 
  distinct(tr_mm10, mm10, gene_type, .keep_all = TRUE) %>% 
  filter(duplicated(mm10)) %>%
  count(mm10,sort = TRUE) # 48 duplicates

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

## find unique pairs of seq 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 %>%  # [929,168] GRs / [90,056] sRNA_id -> [97,367] sRNA_id2 
  left_join(uniq_seq)
```
### v. piRNABank. import the piRNA sequences aligned to genome 

```{r import piRNABank}
piRNAbank_mm10 <- here("mouse_data","piRNABank", "piBnk_GRCm38_v34_fin.fasta") %>% 
  Biostrings::readDNAStringSet()

piRNAbank_mm10_tib <- piRNAbank_mm10 %>% # 55,548 GRs
  as.character() %>% 
  enframe(value = "seq_piBn") %>% 
  mutate(sRNA_type ="piRNA",
    bpairs_piR = str_length(seq_piBn)) %>% 
  distinct(name, seq_piBn, .keep_all = TRUE) %>% 
  arrange(desc(bpairs_piR))
```
### vi. piRNAbank. make Genomic Ranges and remove duplicates from GRCm38 

```{r import piRNABank GRanges}
# piRNABank. import the Genomic Ranges and filter them -----
  piRNAbank_mm10_ranges <- here("mouse_data", "piRNABank", "piBnk_GRCm38_v34_fin.bed") %>% 
  read_bed() %>% 
  as_tibble() %>% 
  arrange(desc(width)) %>% # [55,905] GRs / [39,380] sequences
  filter(width < 100) %>% # [55,905] -> [55,903] GRs / [39,380] sequences
  as_granges() %>% 
  keepStandardChromosomes(pruning.mode = "coarse") # [55,903] -> [54,677] GRs / [39,380]seq

transcripts_pi_mm10 <- Views(BSgenome.Mmusculus.UCSC.mm10, piRNAbank_mm10_ranges) %>%
  as_granges() %>%
  keepStandardChromosomes(pruning.mode = "coarse") %>% 
  as_tibble()

# we need to apply a second width filter at 68 (so as to have sequences of pirnas with max 38)
# as we know that piRNAs are ~32 base pairs
transcripts_pi_mm10 <- transcripts_pi_mm10 %>% filter(width < 68) # [54,677] -> [54,676] GRs / [39,380] names, [39,691] sequences
transcripts_pi_mm10 %>% count(name) %>% nrow #> 39380 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
# or a missmatch
sequen_pi_false <- transcripts_pi_mm10 %>% 
  as_tibble() %>% 
  left_join(piRNAbank_mm10_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_mm10_ranges %>% 
  as_tibble() %>%
  unite(col = "seq_s",seqnames:strand, sep = "_") %>% 
  count(seq_s) %>% 
  filter(n > 1) %>% 
  .$seq_s %>% 
  map(~sequen_pi_false %>% 
  filter(seq_s == .x)) %>% 
  bind_rows()

transcripts_pi_mm10_clean <- transcripts_pi_mm10 %>% 
  as_tibble() %>% 
  left_join(piRNAbank_mm10_tib) %>% 
  arrange(desc(width)) %>% 
  mutate(sequences_true = (dna == seq_piBn)) %>% 
  filter(sequences_true == TRUE) %>% # [54,676] -> [54,014] GRs / [39,380] -> [39,380] sequences 
  select(-score, -seq_piBn, -bpairs_piR, -sequences_true) %>% 
  unite(col = "seq_piBNK",seqnames:strand, sep = "_")

transcripts_pi_mm10_clean %>% count(name) %>% nrow #> 39380 piRNAs final piRNABANK
#> 54,014 Genomic ranges

```
### vii. RNAcentral. + piRNABank. make annotation tibble
create a tibble with that information of RNAcentral and piRNAbank sequences and IDs
```{r annotation tibble}
mm10_piBAnk_RCent <- transcripts_GR %>% 
  left_join(piRNAbank_mm10_tib, by = c("dna" = "seq_piBn"))

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

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

mm10_piBAnk_RCent%>%
  filter(gene_type != "piRNA",sRNA_type == "piRNA") %>% 
  count(gene_type)


# in case that in piRNAbank a pirna is the same but with different type in RNAcentral
# we will keep the gene type of RNAcentral 

concated_mm10_piBAnk <- mm10_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
    )
)


concated_mm10_piBAnk %>%
  filter(!is.na(name)) %>%
  filter(!gene_type == "piRNA") %>% 
  count(gene_type)

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

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

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

## checking for the piRNA values
is.na(concated_mm10_piBAnk$seq_id) %>% table

## checking for duplicates
concated_mm10_piBAnk %>% 
  filter(duplicated(seq_id)) %>% 
  arrange(name)

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

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

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

fasta_tr_ms_tbl %>% 
  filter(mm10 %in% dupl_seqs$dna)

```
### viii. RNAcentral + piRNABank. generation of GRanges

```{r GRanges generation}
concated_mm10_piBAnk  # df with combined sequences piRNAbank+RNAcentral
transcripts_pi_mm10_clean # has all alignments from piRNAbank 

c_piBNK_RCent <- concated_mm10_piBAnk %>%
  full_join(transcripts_pi_mm10_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("mouse_data", "GRCm38", "gencode.vM25.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:843712, CDS:528978

# 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= 760315,
                                   # CDS:protein_coding=527274

#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 = 932,645
  filter(gene_type == "piRNA") %>% # GR = 851,947
  find_overlaps_directed(gene_annot_exon_CDS_prot_cod, maxgap = -1L, 
                         minoverlap = 1L, suffix = c("_piRNAs", "_genes")) %>% 
  select(starts_with(c("gene_","typ"))) %>% 
  as_tibble() %>% 
  unite(col = "seq_piBNK", seqnames:strand, sep = "_") %>% 
  distinct(seq_piBNK, .keep_all = TRUE) # GR = 6,886 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) %>%  # 932,645 -> 925,759
  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 29 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 29 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 = 845,061
  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 = 10,700 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.Mmusculus.UCSC.mm10, 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("mouse_data", "sncRNA_piRNBnk_RNACent_GRCm38_v34.fa"))

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

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

gtf_piB_RCentr %>% 
  write_gff2(file.path("mouse_data", "sncRNA_piRNBnk_RNACent_GRCm38_v34.gtf"))

```
Until now, we have prepared the files for annotation and quantification (GTF, FASTA) 
of smallRNA sequencing samples of mouse 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("mouse_CM")
data_sets %>% 
  file.path("Datasets_analysis",.) %>% 
  map(~dir.create(., recursive = T))
```
## 6. Optional Make the indexes for STAR and salmon
### iv. Indexes for STAR and Salmon with Spike-ins

```{bash  mouse salmon index}
docker run --rm -ti -v $(pwd):/home/my_data  congelos/sncrna_workflow
mkdir my_data/mouse_data/indexes

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

sed -i.bak -e 's/>//g' my_data/mouse_data/decoys_GRCm38.txt

#concat genome transcriptome
fcat my_data/mouse_data/sncRNA_piRNBnk_RNACent_GRCm38_v34.fa my_data/mouse_data/GRCm38/GRCm38.primary_assembly.genome.fa > my_data/mouse_data/gentrome_sncRNA_piRNBnk_RNACent_GRCm38_v34.fasta

pigz --best -p 10 my_data/mouse_data/gentrome_sncRNA_piRNBnk_RNACent_GRCm38_v34.fasta 

mkdir  my_data/mouse_data/indexes/GRCm38_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/mouse_data/gentrome_sncRNA_piRNBnk_RNACent_GRCm38_v34.fasta.gz \
-d my_data/mouse_data/decoys_GRCm38.txt \
-i my_data/mouse_data/indexes/GRCm38_v34_salmon -k 15 -p 6

exit
```


############################  The following code need to be revised more

Following this workflow the files for annotation and quantification of small-RNA samples have been prepared. 
Afterwards, the steps in the  pre-processing of the samples, alignment, quantification and calculation of transcript abundances could be followed.

## 7. Provide extra information regarding genomic locations, genes, transcripts, for the gtf
### i. Load libraries

```{r bumphunter libraries}
suppressPackageStartupMessages({
library('TxDb.Mmusculus.UCSC.mm10.knownGene')
library('org.Mm.eg.db')
library('bumphunter')
library('BiocParallel')
library('stats')   
})
```
### ii. import regions of transcripts

```{r transcript regions bumphunter}
genes <- annotateTranscripts(TxDb.Mmusculus.UCSC.mm10.knownGene, annotation="org.Mm.eg.db") %>% 
  keepStandardChromosomes(pruning.mode="coarse") %>% arrange(seqnames)

piRNAbank_rCentral_gtf <- read_gff2("mouse_data/ncRNA_transcripts_100bp_RNA_Central_piRNAbank_mm10.gtf")

identical(genes %>% seqlevels(), piRNAbank_rCentral_gtf %>% seqlevels())

piRNAbank_rCentral_gtf %>% length()

map(piRNAbank_rCentral_gtf %>% seqlevels() %>% purrr::set_names(),
    ~piRNAbank_rCentral_gtf %>% 
        filter(seqnames == .x) %>% 
        length()) %>% bind_rows() %>%  
  pivot_longer(cols = chr1:chrM) %>% 
  arrange(desc(value))

# we will parallelize per chr.

 if(.Platform$OS.type == "windows") {
mt_param <- SnowParam()
 } else{
mt_param <- MulticoreParam()
 }
# we will work with 10 workers 
mt_param <- MulticoreParam(workers = 8)

# simple function which takes lists of Granges and the chromosome 
# name to select from each list
matchGenes_fun <- function(our_Grs, genes_GRs){
  suppressPackageStartupMessages({
    library('dplyr')
    library('bumphunter')
  })
  message("working on matchGenes")
  matchGenes(our_Grs, genes_GRs, 
    type = "any", promoterDist = 2500, 
    skipExons = FALSE, verbose = TRUE) %>% as_tibble()
}

genes_chr <- map(genes %>% seqlevels() %>% purrr::set_names(), 
  ~genes %>% filter(seqnames == .x))

gen_test <- genes_chr[c("chrM","chrY")]

piR_chr <- map(piRNAbank_rCentral_gtf %>% 
    seqlevels() %>% 
    purrr::set_names(), ~piRNAbank_rCentral_gtf %>% 
        filter(seqnames == .x))

piR_test <- piR_chr[c("chrM","chrY")]

res_chr <- bpmapply(matchGenes_fun,
  piR_chr, genes_chr, USE.NAMES=TRUE, SIMPLIFY = FALSE,
    BPREDO=list(), BPPARAM = mt_param) 

res_chr <- bind_rows(res_chr) %>% 
  bind_cols(as_tibble(piRNAbank_rCentral_gtf)) %>% 
  dplyr::select(name:subjectHits, gene_id, 
    gene_type, sRNA_id, source, seq_RNA) %>% 
  write_tsv("mouse_data/gene_regions_piRNAbank_rCentral.txt")
``` 
## 8. Find multimapping piRNAs

```{r multimapping piRNAs}
multi_test <- piRNAbank_rCentral_gtf %>%
    plyranges::select(gene_id, seq_RNA, gene_type) %>% join_overlap_inner_directed(plyranges::select(piRNAbank_rCentral_gtf, gene_id, gene_type, seq_RNA)) %>% 
  arrange(seqnames)

multi_test %>% 
  filter(gene_type.x == "piRNA", 
    !gene_id.x == gene_id.y ) %>%
  as_tibble() %>% 
  count(gene_type.x, gene_type.y, sort = T) %>% 
  write_tsv("genomic_locations_stats_multi.txt")

piRNAbank_rCentral_gtf %>% 
  filter(gene_type == "piRNA") %>% 
  plyranges::select(-c(score, phase, source, type)) %>% 
  as_tibble() %>% 
  unite(col = "seq_s",seqnames:strand, sep = "_") %>% 
  count(gene_id, sort = T) %>% write_tsv("genomic_locations_stats_multi_piRNA.txt")
```
## 9. Find how many piRNAs are in common and uncommon in piRNABank and RNAcentral in the new gtf 

```{r}
c_piBNK_RCent %>% distinct(seq_id, .keep_all = T)

c_piBNK_RCent %>% distinct(seq_id, .keep_all = T) %>% filter(is.na(sRNA_type),gene_type == "piRNA")

c_piBNK_RCent %>% distinct(seq_id, .keep_all = T) %>% filter(!is.na(sRNA_type),is.na(sRNA_id),gene_type == "piRNA")

c_piBNK_RCent %>% distinct(seq_id, .keep_all = T) %>% filter(!is.na(name),!is.na(sRNA_id), gene_type== "piRNA")
```
## 10. Find which smallRNAs are inside Trasposable Elements

We have downloaded a gtf file with the information about genomic regions of 
Transposable Elements for human genome: http://labshare.cshl.edu/shares/mhammelllab/www-data/TEtranscripts/TE_GTF/
more precisely: [GRCm38_Ensembl_rmsk_TE.gtf.gz](http://labshare.cshl.edu/shares/mhammelllab/www-data/TEtranscripts/TE_GTF/GRCm38_Ensembl_rmsk_TE.gtf.gz)

```{r Trasposable Elements annotation}
TEs <- read_gff2("GRCm38_Ensembl_rmsk_TE.gtf.gz") %>% 
  plyranges::select("TE_gene_id" = gene_id, "TE_transcript_id" = transcript_id, 
   "TE_family_id" = family_id, "TE_class_id" = class_id) %>% 
  keepStandardChromosomes(pruning.mode = "coarse") %>% 
  arrange(seqnames) 

piRNAbank_rCentral_gtf %>% 
  plyranges::select(gene_id, sRNA_id,gene_type, seq_RNA) %>% 
  find_overlaps_directed(TEs) %>% 
  write_gff2("TEs_piRNAbank_rCentral.gtf")

piRNAbank_rCentral_gtf %>% 
  join_overlap_left_directed( piRNAbank_rCentral_gtf %>% 
  find_overlaps_directed(TEs)) %>% length()

piRNAbank_rCentral_gtf %>% 
  find_overlaps_directed(TEs) %>% 
  plyranges::reduce_ranges_directed() %>% length()

```