01_snv_analysis.Rmd

---
title: "CLiP postprocessing - deletion/SNV based"
author: "Katharina Hembach, updated by Izaskun Mallona"
date: "16.08.2021"
output: 
  html_document:
    toc: true
    toc_float: true
---
  
```{r, echo = FALSE, warning = FALSE, message = FALSE}
knitr::opts_chunk$set(autodep = TRUE, cache = FALSE,
                      cache.lazy = FALSE,
                      dev = "png",
                      dev.args = list(png = list(type = "cairo")),
                      warning = FALSE, message = FALSE)
```


### Load packages

```{r, message = FALSE}
## library(here)
library(Rsamtools)
library(rtracklayer)
library(dplyr)
library(BSgenome.Hsapiens.UCSC.hg38)
library(eulerr)
library(ggplot2)
library(GenomicFeatures)
library(ggrepel)
library(ggpubr)
library(stringr)
library(RColorBrewer)
library(GenomicAlignments)
library(memes)
```

```{r}
WD <- file.path('/home', 'imallona', 'polymenidou_manu_clip/', 'armor_output')
```

# Load data
```{r load-files}
samples <- c("WT",  "6M", "RBDm")
bams <- lapply(samples, function(x) 
  list.files(file.path(WD, "BAM_deduplicated_clip5"), 
             pattern = paste0("20200123.A-", x, "2*_R2_deduplicated.bam$")))
names(bams) <- samples

## reference genome
genome <- BSgenome.Hsapiens.UCSC.hg38
seqlevelsStyle(genome) <- 'UCSC' ## "Ensembl"
## gene annotations
gtf_file <- file.path(WD, "..", "reference", "gencode.v39.annotation.gtf")
gtf <- import(gtf_file)
table(gtf$type)
genes <- gtf[gtf$type =="gene"]
```

# Compute pileup

We generate genomewide pileup tables for each bam file. The pileup is used to identify the location and frequency of mutations (including deletions).

```{r pileup-mut-table}

if (!file.exists('pu_res.rds')) {
    pu_res <- lapply(samples,function(x) NULL)  %>% setNames(samples)
    min_allel_depth <- 1 ## at least 1 read with deletion per sample

    for(s in samples){
                                        # sbp <- ScanBamParam(which = mp[[s]])
        ## require at least 3 reads with any mutation per position
        ## show all mutations with at least 1 supporting read
        pp <- PileupParam(max_depth = 30000, include_deletions = TRUE, 
                          include_insertions=FALSE,
                          min_nucleotide_depth=1, min_minor_allele_depth = min_allel_depth)
        ## compute pileup for both of the replicates
        
        bf <- open(BamFile(file.path(WD, "BAM_deduplicated_clip5", bams[[s]][1]), 
                           yieldSize=5e4, 
                           index = file.path(WD, "BAM_deduplicated_clip5", 
                                             paste0(bams[[s]][1], ".bai"))))
        
        pu <- list("rep1" = data.frame(), "rep2" = data.frame())
        repeat {
            res <- pileup(bf)
            if(nrow(res) == 0L)
                break
            res <- res %>% dplyr::filter(nucleotide == "-")
            pu[["rep1"]] <- rbind(pu[["rep1"]], res)
        }
        close(bf)
        
        bf <- open(BamFile(file.path(WD, "BAM_deduplicated_clip5", bams[[s]][2]), 
                           yieldSize=5e4,
                           index = file.path(WD, "BAM_deduplicated_clip5", 
                                             paste0(bams[[s]][2], ".bai"))))
        
        repeat {
            res <- pileup(bf)
            if(nrow(res) == 0L)
                break
            res <- res %>% dplyr::filter(nucleotide == "-")
            pu[["rep2"]] <- rbind(pu[["rep2"]], res)
        }
        close(bf)
        lapply(pu, dim)

        pu_res[[s]] <- pu
    }

    saveRDS(pu_res, 'pu_res.rds')
} else {
    pu_res <- readRDS('pu_res.rds')
}
```

We compute the fraction of reads with deletion, and remove positions with >50% deletions, these are either true deletions compared to the reference genome or positions with very low read coverage.

```{r read-genome-coverage}
## import BAM files to get genome read coverage
if (!file.exists( "clip5_BAM_coverage.rds")) {
  fs <- dir(path = file.path(WD, "BAM_deduplicated_clip5"), pattern = "_R2_deduplicated.bam$", 
            full.names = TRUE)
  names(fs) <- str_extract(basename(fs), 
                         "(?<=20200123.A-).*(?=_R2_deduplicated.bam)")
  cov <- list()
  for (sample in names(fs)){
    ga <- readGAlignments(fs[sample])
    cov[[sample]] <- coverage(ga)
    rm(ga)
  }
  gc()
  saveRDS(cov,  "clip5_BAM_coverage.rds")
} else{
  cov <- readRDS("clip5_BAM_coverage.rds")
}
```

```{r add-del-perc}
lapply(pu_res, function(x) lapply(x, nrow))

pu_res_cov <- list() 
for(i in names(pu_res)){
  m <- match(i, names(cov))
  a <- pu_res[[i]]
  a[["rep1"]]$nreads <- cov[[m]][GRanges(seqnames = a[["rep1"]]$seqnames, 
                               range = IRanges(start = a[["rep1"]]$pos, 
                                               end = a[["rep1"]]$pos), 
                               strand = a[["rep1"]]$strand)] %>% 
    unlist %>% as.vector 
  a[["rep1"]]$del_frac <- a[["rep1"]]$count / a[["rep1"]]$nreads 
  a[["rep2"]]$nreads <- cov[[m+1]][GRanges(seqnames = a[["rep2"]]$seqnames, 
                               range = IRanges(start = a[["rep2"]]$pos, 
                                               end = a[["rep2"]]$pos), 
                               strand = a[["rep2"]]$strand)] %>% 
    unlist %>% as.vector
  a[["rep2"]]$del_frac <- a[["rep2"]]$count / a[["rep2"]]$nreads 
  pu_res_cov[[i]] <- a
}
lapply(pu_res_cov, function(x) lapply(x, function(y) summary(y$del_frac)))
## fraction of positions with > 50% deletions
lapply(pu_res_cov, function(x) lapply(x, function(y) sum(y$del_frac > 0.5)/nrow(y)))
## what is the number of reads at these positions?
lapply(pu_res_cov, function(x) lapply(x, function(y) y[y$del_frac > 0.5, "nreads"] %>% summary))

## remove all positions with > 50% deletions
pu_res_cov <- lapply(pu_res_cov, function(x) lapply(x, function(y) y[y$del_frac <= 0.5,]))
lapply(pu_res_cov, function(x) lapply(x, nrow))
```



```{r deletion-pos}
# ## check the number of reads supporting each deletion in the two replicates and filter positions
# ## that only have deletions in one replicate but not the other!lapply(pu_res,nrow)
pu_res_cov[["RBDm"]][["rep1"]] %>% head
pu_res_cov[["RBDm"]][["rep2"]] %>% head

pu_merge <- lapply(pu_res_cov, function(x){
  x[["rep1"]] %>% 
  full_join(x[["rep2"]], by = c("seqnames", "pos", "strand", "nucleotide"), 
            suffix = c("1", "2")) 
  })
lapply(pu_merge, nrow)

for(s in names(pu_merge)){
  print(s)
  print(paste0("number of positions: ", nrow(pu_merge[[s]])))
  ## positions with deletions in both replicates
  print(paste0("deletion in both reps: ", 
          pu_merge[[s]] %>% 
            dplyr::filter(!is.na(count1) & !is.na(count2)) %>% nrow))
  ## positions with deletions supported by at least 2 reads in one of the replicates
  print(paste(">=2 reads in one rep: ", 
          pu_merge[[s]] %>% dplyr::filter((count1 > 1 & !is.na(count2)) |
                                            (!is.na(count1) & count2 > 1)) %>% nrow))
  ## positions with deletions supported by at least 2 reads in both replicates
  print(paste0(">=2 reads in both reps: ",