---
title: "Importing & Processing Data"
package: BRGenomics
output:
  BiocStyle::html_document:
    toc: true
    toc_float: true
  BiocStyle::pdf_document:
    toc: true
vignette: |
  %\VignetteIndexEntry{Importing and Processing Data}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
---

BRGenomics provides several functions for conveniently importing and processing 
BAM, bigWig, bedGraph files.

# Importing BAM Files

The `import_bam()` function provides a number of options for filtering and 
processing bam files. BRGenomics includes an example BAM file with a small 
number of reads from the included PRO-seq data. The file's local location can 
be found (on your computer) as follows:

```{r, message = FALSE}
library(BRGenomics)
```

```{r}
bfile <- system.file("extdata", "PROseq_dm6_chr4.bam", 
                     package = "BRGenomics")
```

Because PRO-seq data is sequenced in the 3'-to-5' direction of the original RNA 
molecule, we'll use the `revcomp` option to reverse-complement all the input 
reads. We'll also set a minimum MAPQ score of 20:

```{r}
ps_reads <- import_bam(bfile, mapq = 20, revcomp = TRUE, paired_end = FALSE)
ps_reads
```

By default, `import_bam()` combines identical reads into the same range, and 
the `score` metadata column indicates the number of perfectly-overlapping 
alignments. This means that the total number of alignments (reads) is equal to 
the sum of the score:

```{r}
sum(score(ps_reads))
```

Alternatively, you can import each read as its own range by setting `field = 
NULL`:

```{r}
reads_expanded <- import_bam(bfile, mapq = 20, revcomp = TRUE, 
                             field = NULL, paired_end = FALSE)
ps_reads[1:8]
reads_expanded[1:8]
```

Notice that reads 5-7 are now identical, rather than combined into a single 
range with a score = 3.

```{r}
length(reads_expanded) == sum(score(ps_reads))
```

Many BRGenomics function have a `field` argument, and setting `field = NULL` 
will treat each range has a single read. 

## Example: Importing PRO-seq BAM files at Basepair Resolution

We can use the `import_bam()` function to extract the positions of interest 
from BAM files. Below, we construct an import function for PRO-seq data that 
returns a basepair-resolution GRanges object.

In PRO-seq, a "run-on" reaction is performed in which actively engaged RNA 
polymerases incorporate a biotinylated nucleotide at the 3' end of a nascent 
RNA. Our base of interest is therefore the base immediately preceding the RNA 3' 
end, as this was the original position of a polymerase active site.

The processing options in `import_bam()` are applied in the same order that 
they're listed in the documentation page. Following this order, we will apply 
the options:

1. Filter reads by a minimum MAPQ score 
2. Take the reverse complement
3. Shift reads upstream by 1 base
4. Extract the 3' base

```{r}
ps <- import_bam(bfile, 
                 mapq = 20, 
                 revcomp = TRUE,
                 shift = -1,
                 trim.to = "3p",
                 paired_end = FALSE)
ps
```

_Note that for paired-end data, `import_bam()` will automatically filter 
unmatched read pairs._

Notice that the number of ranges in `ps` is not the same as for `ps_reads`, in 
which we imported the entire read lengths:

```{r, collapse = TRUE}
length(ps_reads)
length(ps)
```

This is because identical positions are collapsed together after applying the 
processing options. However, we can check that all of the same reads are 
represented:

```{r, collapse = TRUE}
sum(score(ps)) == sum(score(ps_reads))
```

And we can check that collapsing identical positions has created a 
basepair-resolution GRanges object:

```{r, collapse = TRUE}
isBRG(ps)
```

## Pre-formatted Input Functions

For convenience, we've included several functions with default options for 
several kinds of data, including `import_bam_PROseq()`, `import_bam_PROcap()`, 
and `import_bam_ATACseq()`, the latter of which corrects for the 9 bp offset 
between Tn5 insertion sites.^[Jason D. Buenrostro, Paul G. Giresi, Lisa C. 
Zaba, Howard Y. Chang, William J. Greenleaf (2013). Transposition of native 
chromatin for fast and sensitive epigenomic profiling of open chromatin, 
dna-binding proteins and nucleosome position. \emph{Nature Methods} 10: 
1213–1218. \url{https://doi.org/10.1038/nmeth.2688}]

## Example: Converting BAMs to bigWigs

In conjunction with export functions from the `r Biocpkg("rtracklayer")` 
package, we can use the functions described above to write a post-alignment 
pipeline for generating bigWig files for PRO-seq data:

```{r, eval=FALSE}
# import bam, automatically applying processing steps for PRO-seq
ps <- import_bam_PROseq(bfile, mapq = 30, paired_end = FALSE)

# separate strands, and make minus-strand scores negative
ps_plus <- subset(ps, strand == "+")
ps_minus <- subset(ps, strand == "-")
score(ps_minus) <- -score(ps_minus)

# use rtracklayer to export bigWig files
export.bw(ps_plus, "~/Data/PROseq_plus.bw")
export.bw(ps_minus, "~/Data/PROseq_minus.bw")
```

## Performance Considerations

For single-ended bam files, import is much faster if the bam files are sorted 
and indexed (i.e. by `samtools index`).

For paired-end files, we assume that collating (`samtools collate`) or sorting
by name is faster.

Additionally, while single-ended files can often be imported "all at once"
(particularly if sorted and indexed), processing paired-end data is more memory
intensive, and requires breaking up the file into chunks for processing. For 
this, use the `yieldSize` argument.

For example, to process 500 thousands reads at a time, set the 
`yieldSize = 5e5`.

# Importing bedGraphs and bigWigs

bedGraph and bigWig files are efficient and portable, but unstranded 
representations of basepair-resolution genomics data. 

As compared to `rtracklayer::import.bedGraph()`, the BRGenomics function 
`import_bedGraph()` imports both plus-strand and minus-strand files as a single 
object, and has options for filtering out odd chromosomes, mitochondrial 
chromosomes, and sex chromosomes.

```{r}
# local paths to included bedGraph files
bg.p <- system.file("extdata", "PROseq_dm6_chr4_plus.bedGraph",
                    package = "BRGenomics")
bg.m <- system.file("extdata", "PROseq_dm6_chr4_minus.bedGraph",
                    package = "BRGenomics")

import_bedGraph(bg.p, bg.m, genome = "dm6")
```

The `import_bigWig()` function provides the same added functionality as 
compared to `rtracklayer::import.bw()`, but also removes run-length compression 
and returns a basepair-resolution GRanges object by default.

```{r}
# local paths to included bigWig files
bw.p <- system.file("extdata", "PROseq_dm6_chr4_plus.bw",
                    package = "BRGenomics")
bw.m <- system.file("extdata", "PROseq_dm6_chr4_minus.bw",
                    package = "BRGenomics")

import_bigWig(bw.p, bw.m, genome = "dm6")
```

Conversion to a basepair-resolution GRanges object can be turned off by setting
`makeBRG = FALSE`.

# Merging GRanges Data

Biological replicates are best used to independently reproduce and measure 
effects, and therefore we often want to handle them separately. However, there 
are times when combining replicates can allow for more sensitive measurements, 
assuming that the replicates are highly concordant.

The `mergeGRangesData()` function can be used to combine basepair-resolution 
GRanges objects.

We'll break up the included PRO-seq data into a list of toy datasets:

```{r}
ps_list <- lapply(1:6, function(i) ps[seq(i, length(ps), 6)])
names(ps_list) <- c("A_rep1", "A_rep2", 
                    "B_rep1", "B_rep2",
                    "C_rep1", "C_rep2")
```

```{r}
ps_list[1:2]
names(ps_list)
```

We can pass a list of GRanges objects directly as an argument:

```{r}
mergeGRangesData(ps_list, ncores = 1)
```

Or we can pass any number of individual GRanges objects as arguments:

```{r}
merge_ps <- mergeGRangesData(ps_list[[1]], ps_list[[2]], ps, ncores = 1)
merge_ps
```

_Note that the output is also a basepair-resolution GRanges object:_

```{r, collapse = TRUE}
isBRG(merge_ps)
```

\ 

## Merging Replicates

The `mergeReplicates()` function makes combining replicates particularly 
simple:

```{r}
mergeReplicates(ps_list, ncores = 1)
```

