---
title: "Working with aligned nucleotides (WORK-IN-PROGRESS!)"
author:
- name: "Hervé Pagès"
- name: "Halimat C. Atanda"
  affiliation: "Vignette translation from Sweave to Rmarkdown / HTML"
date: "Last modified: January 2014; Compiled:`r format(Sys.Date(), '%B %d %Y')`"
package: GenomicAlignments
output: BiocStyle::html_document
vignette: >
  %\VignetteIndexEntry{Working with aligned nucleotides (WORK-IN-PROGRESS!)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

# Introduction

This vignette belongs to the `r Biocpkg("GenomicAlignments")` package. It
illustrates how to use the package for working with the nucleotide content of
aligned reads.

After the reads generated by a high-throughput sequencing experiment have been
aligned to a reference genome, the questions that are being asked about these
alignments typically fall in two broad categories: **positional only** and
**nucleotide-related**.

**Positional only** questions are about the position of the alignments with
respect to the reference genome. Note that the position of an alignment is
actually better described in terms of genomic ranges (1 range for an alignment
with no gaps, 2 or more ranges for an alignment with gaps). Knowing the ranges
of the alignments is sufficient to perform common tasks like *read counting* or
for *computing the coverage*. *Read counting* is the process of counting the
number of aligned reads per gene or exon and is typically performed in the
context of a differential analysis. This task can be accomplished with the
`summarizeOverlaps` function provided in the package and is explained in details
in the "Counting reads with summarizeOverlaps" vignette (also located in this
package). *Computing the coverage* is often the preliminary step to peak
detection (ChIP-seq analysis) or to a copy number analysis. It can be
accomplished with the `coverage` function. See `?'coverage-methods'` for more
information.

**Nucleotide-related** questions are about the nucleotide content of the
alignments. In particular, how this content compares to the corresponding
nucleotides in the reference genome. These questions typically arise in the
context of small genetic variation detection between one or more samples and a
reference genome.

The `r Biocpkg("GenomicAlignments")` package provides a suite of low- to
mid-level tools for dealing with **nucleotide-related** questions about the
alignments. In this vignette, we illustrate their use on the single-end and
paired-end reads of an RNA-seq experiment.

Note that these tools do NOT constitute a complete variant toolbox. If this is
what you're looking for, other *Bioconductor* packages might be more
appropriate. See the GeneticVariability and SNP views at this URL
<https://bioconductor.org/packages/release/BiocViews.html#___AssayDomain> for a
complete list of packages that deal with small genetic variations. Most of them
provide tools of higher level than the tools described in this vignette. See for
example the `r Biocpkg("VariantTools")` and `r Biocpkg("VariantAnnotation")`
packages for a complete variant toolbox (including variant calling
capabilities).

# Load the aligned reads and their sequences from a BAM file

In this section, we illustrate how aligned reads and their sequences can be
loaded from a BAM file. The reads we're going to use for this are paired-end
reads from a published study by Zarnack et al., 2012. A subset of these reads
are stored in the BAM files located in the data `r Biocpkg("RNAseqData.HNRNPC.bam.chr14")` data package. The package contains 8 BAM
files, 1 per sequencing run:

```{r message=FALSE}
library(RNAseqData.HNRNPC.bam.chr14)
bamfiles <- RNAseqData.HNRNPC.bam.chr14_BAMFILES
names(bamfiles) # the names of the runs
```

Each BAM file was obtained by (1) aligning the reads (paired-end) to the full
hg19 genome with TopHat2, and then (2) subsetting to keep only alignments on
chr14. See `?RNAseqData.HNRNPS.bam.chr14` for more information about this data
set.

As a preliminary step, we check whether the BAM files contain single- or
paired-end alignments. This can be done with the `quickBamFlagSummary` utility
from the `r Biocpkg("Rsamtools")` package:

```{r message=FALSE}
library(Rsamtools)
quickBamFlagSummary(bamfiles[1], main.groups.only=TRUE)
```

This confirms that all the alignments in the 1st BAM file (run ERR127306) are
paired-end. This means that we should preferably load them with the
`readGAlignmentsPairs` function from the `r Biocpkg("GenomicAlignments")`
package. However for the purpose of keeping things simple, we will ignore the
pairing for now and load only the alignments corresponding to the first segment
of the pairs. We will use the `readGAlignments` function from the `r Biocpkg("GenomicAlignments")` package for this, together with a *ScanBamParam*
object for the filtering. See `?ScanBamParam` in the `r Biocpkg("Rsamtools")`
package for the details.

Furthermore, while preparing the *ScanBamParam* object to perform the filtering,
we'll also get rid of the PCR or optical duplicates (flag bit 0x400 in the SAM
format, see the SAM Spec [^1] for the details), as well as reads not passing
quality controls (flag bit 0x200 in the SAM format).

Finally we also request the read sequences (a.k.a. the *segment sequences* in the SAM Spec, stored in the SEQ field) via the *ScanBamParam* object:

```{r}
flag1 <- scanBamFlag(isFirstMateRead=TRUE, isSecondMateRead=FALSE, isDuplicate=FALSE, isNotPassingQualityControls=FALSE)
param1 <- ScanBamParam(flag=flag1, what="seq")
```

We're now ready to load the alignments and their sequences. Note that we use
`use.names=TRUE` in order to also load the *query names* (a.k.a. the *query
template names* in the SAM Spec, stored in the QNAME field) from the BAM file.
`readGAlignments` will use them to set the names of the returned object:

```{r message=FALSE}
library(GenomicAlignments)
gal1 <- readGAlignments(bamfiles[1], use.names=TRUE, param=param1)
```

This returns a *GAlignments* object. The read sequences are stored in the `seq` metadata column of the object:

```{r}
mcols(gal1)$seq
```

# Compute the *original query* sequences

Because the BAM format imposes that the read sequence is "reverse complemented"
when a read aligns to the minus strand, we need to "reverse complement" it again
to restore the *original query sequences* (i.e. the sequences before alignment,
that is, as they can be seen in the FASTQ file assuming that the aligner didn't
perform any hard-clipping on them):

```{r}
oqseq1 <- mcols(gal1)$seq
is_on_minus <- as.logical(strand(gal1) == "-")
oqseq1[is_on_minus] <- reverseComplement(oqseq1[is_on_minus])
```

Because the aligner used to align the reads can report more than 1 alignment per
read (i.e. per sequence stored in the FASTQ file), we shouldn't expect the names
of `gal1` to be unique:

```{r}
is_dup <- duplicated(names(gal1))
table(is_dup)
```

However, sequences with the same *query name* should correspond to the same
*original query* and therefore should be the same. Let's do a quick sanity
check:

```{r}
dup2unq <- match(names(gal1), names(gal1))
stopifnot(all(oqseq1 == oqseq1[dup2unq]))
```

Finally, let's `oqseq1` reduce to one *original query sequence* per unique
*query name* (like in the FASTQ file containing the 1st end of the unaligned
reads):

```{r}
oqseq1 <- oqseq1[!is_dup]
```

# Mismatches, indels, and gaps

Because the aligner possibly tolerated a small number of mismatches, indels,
and/or gaps during the alignment process, the sequences in `mcols(gal1)$seq`
generally don't match exactly the reference genome.

The information of where indels and/or gaps occur in the alignments is
represented in the CIGAR strings. Let's have a look at these string. First the
most frequent cigars:

```{r}
head(sort(table(cigar(gal1)), decreasing=TRUE))
```

Then a summary of the total number of insertions (I), deletions (D), andcgaps (N):

```{r}
colSums(cigarOpTable(cigar(gal1)))
```

This tells us that the aligner that was used supports indels (I/D) and junction reads (N).

Finally we count and summarize the number of gaps per alignment:

```{r}
table(njunc(gal1))
```

Some reads contain up to 3 gaps (i.e. span 3 splice junctions).

# Put the read sequences and reference sequences "side by side"

TODO (with `sequenceLayer`)

# OLD STUFF (needs to be recycled/updated)

## Load paired-end reads from a BAM file

BAM file untreated3_chr4.bam (located in the `r Biocpkg("pasillaBamSubset")`
data package) contains paired-end reads from the "Pasilla" experiment and
aligned against the dm3 genome (see `?untreated3_chr4` in the `r Biocpkg("pasillaBamSubset")` package for more information about those reads).
We use the `readGAlignmentPairs` function to load them into a *GAlignmentPairs* object:

```{r message=FALSE}
library(pasillaBamSubset)
flag0 <- scanBamFlag(isDuplicate=FALSE, isNotPassingQualityControls=FALSE)
param0 <- ScanBamParam(flag=flag0)
U3.galp <- readGAlignmentPairs(untreated3_chr4(), use.names=TRUE, param=param0)
head(U3.galp)
```

The `show` method for *GAlignmentPairs* objects displays two ranges columns, one
for the *first* alignment in the pair (the left column), and one for the *last*
alignment in the pair (the right column). The strand column corresponds to the
strand of the *first* alignment.

```{r}
head(first(U3.galp))
head(last(U3.galp))
```

According to the SAM format specifications, the aligner is expected to mark each
alignment pair as *proper* or not (flag bit 0x2 in the SAM format). The SAM Spec
only says that a pair is *proper* if the *first* and *last* alignments in the
pair are "properly aligned according to the aligner". So the exact criteria used
for setting this flag is left to the aligner.

We use `isProperPair` to extract this flag from the object:

```{r}
table(isProperPair(U3.galp))
```

Even though we could do *overlap encodings* with the full object, we keep only the *proper* pairs for our downstream analysis:

```{r}
U3.GALP <- U3.galp[isProperPair(U3.galp)]
```

Because the aligner used to align those reads can report more than 1 alignment
per *original query template* (i.e. per pair of sequences stored in the input
files, typically 1 FASTQ file for the *first* ends and 1 FASTQ file for the
*last* ends), we shouldn't expect the names of `U3.GALP` to be unique:

```{r}
U3.GALP_names_is_dup <- duplicated(names(U3.GALP))
table(U3.GALP_names_is_dup)
```

Storing the *query template names* in a factor will be useful:

```{r}
U3.uqnames <- unique(names(U3.GALP))
U3.GALP_qnames <- factor(names(U3.GALP), levels=U3.uqnames)
```

as well as having the mapping between each *query template name* and its first occurence in `U3.GALP_qnames`:

```{r}
U3.GALP_dup2unq <- match(U3.GALP_qnames, U3.GALP_qnames)
```

Our reads can have up to 1 gap per end:

```{r}
head(unique(cigar(first(U3.GALP))))
head(unique(cigar(last(U3.GALP))))
table(njunc(first(U3.GALP)), njunc(last(U3.GALP)))
```

Like for our single-end reads, the following tables indicate that indels were not allowed/supported during the alignment process:

```{r}
colSums(cigarOpTable(cigar(first(U3.GALP))))
colSums(cigarOpTable(cigar(last(U3.GALP))))
```

# `sessionInfo()`

```{r}
sessionInfo()
```

[^1]: <https://samtools.sourceforge.net/>
