%\VignetteIndexEntry{Overlap encodings}
%\VignetteDepends{pasillaBamSubset, GenomicAlignments, GenomicFeatures, BSgenome.Dmelanogaster.UCSC.dm3, TxDb.Dmelanogaster.UCSC.dm3.ensGene}
%\VignetteKeywords{sequence, sequencing, alignments}
%\VignettePackage{GenomicAlignments}

\documentclass{article}

<<style, eval=TRUE, echo=FALSE, results=tex>>=
BiocStyle::latex()
@

\title{Overlap encodings}
\author{Herv\'e Pag\`es}
\date{Last modified: May 2014; Compiled: \today}

\begin{document}

\maketitle

<<options,echo=FALSE>>=
options(width=100)
.precomputed_results_dir <- "precomputed_results"
.loadPrecomputed <- function(objname)
{
    filename <- paste0(objname, ".rda")
    path <- file.path(.precomputed_results_dir, filename)
    tempenv <- new.env(parent=emptyenv())
    load(path, envir=tempenv)
    get(objname, envir=tempenv)
}
.checkIdenticalToPrecomputed <- function(obj, objname, ignore.metadata=FALSE)
{
    precomputed_obj <- .loadPrecomputed(objname)
    if (ignore.metadata)
        metadata(obj) <- metadata(precomputed_obj) <- list()
    if (!identical(obj, precomputed_obj))
        stop("'", objname, "' is not identical to precomputed version")
}
@

\tableofcontents



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In the context of an RNA-seq experiment, encoding the overlaps between
the aligned reads and the transcripts can be used for detecting those
overlaps that are ``compatible'' with the splicing of the transcript.

Various tools are provided in the \Rpackage{GenomicAlignments} package for
working with {\it overlap encodings}.
In this vignette, we illustrate the use of these tools on the single-end
and paired-end reads of an RNA-seq experiment.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Load reads from a BAM file}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Load single-end reads from a BAM file}

BAM file {\tt untreated1\_chr4.bam} (located in the
\Rpackage{pasillaBamSubset} data package) contains single-end reads from
the ``Pasilla'' experiment and aligned against the dm3 genome (see
\Rcode{?untreated1\_chr4} in the \Rpackage{pasillaBamSubset} package for
more information about those reads):

<<untreated1_chr4>>=
library(pasillaBamSubset)
untreated1_chr4()
@

We use the \Rfunction{readGAlignments} function defined in the
\Rpackage{GenomicAlignments} package to load the reads into a
\Rclass{GAlignments} object.
It's probably a good idea to get rid of the PCR or optical duplicates
(flag bit 0x400 in the SAM format, see the SAM Spec
\footnote{\url{http://samtools.sourceforge.net/}} for the details), as well
as reads not passing quality controls (flag bit 0x200 in the SAM format).
We do this by creating a \Rclass{ScanBamParam} object that we pass to
\Rcode{readGAlignments} (see \Rcode{?ScanBamParam} in the
\Rpackage{Rsamtools} package for the details). Note that we also use
\Rcode{use.names=TRUE} in order to load the {\it query names} (aka {\it query
template names}, see QNAME field in the SAM Spec) from the BAM file
(\Rcode{readGAlignments} will use them to set the names of the returned
object):

<<readGAlignments>>=
library(GenomicAlignments)
flag0 <- scanBamFlag(isDuplicate=FALSE, isNotPassingQualityControls=FALSE)
param0 <- ScanBamParam(flag=flag0)
U1.GAL <- readGAlignments(untreated1_chr4(), use.names=TRUE, param=param0)
head(U1.GAL)
@

Because the aligner used to align those reads can report more than 1
alignment per {\it original query} (i.e. per read stored in the input file,
typically a FASTQ file), we shouldn't expect the names of \Rcode{U1.GAL} to
be unique:

<<U1.GAL_names_is_dup>>=
U1.GAL_names_is_dup <- duplicated(names(U1.GAL))
table(U1.GAL_names_is_dup)
@

Storing the {\it query names} in a factor will be useful as we will see later
in this document:

<<U1.GAL_qnames>>=
U1.uqnames <- unique(names(U1.GAL))
U1.GAL_qnames <- factor(names(U1.GAL), levels=U1.uqnames)
@

Note that we explicitely provide the levels of the factor to enforce their
order. Otherwise \Rcode{factor()} would put them in lexicographic order which
is not advisable because it depends on the locale in use.

Another object that will be useful to keep near at hand is the mapping between
each {\it query name} and its first occurence in \Rcode{U1.GAL\_qnames}:

<<U1.GAL_dup2unq>>=
U1.GAL_dup2unq <- match(U1.GAL_qnames, U1.GAL_qnames)
@

Our reads can have up to 2 gaps (a gap corresponds to an N operation in the
CIGAR):

<<gaps-in-U1.GAL>>=
head(unique(cigar(U1.GAL)))
table(ngap(U1.GAL))
@

Also, the following table indicates that indels were not allowed/supported 
during the alignment process (no I or D CIGAR operations):

<<no-indels-in-U1.GAL>>=
colSums(cigarOpTable(cigar(U1.GAL)))
@


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Load paired-end reads from a BAM file}

BAM file {\tt untreated3\_chr4.bam} (located in the
\Rpackage{pasillaBamSubset} data package) contains paired-end reads from
the ``Pasilla'' experiment and aligned against the dm3 genome (see
\Rcode{?untreated3\_chr4} in the \Rpackage{pasillaBamSubset} package for
more information about those reads).
We use the \Rfunction{readGAlignmentPairs} function to load them
into a \Rclass{GAlignmentPairs} object:

<<readGAlignmentPairs>>=
U3.galp <- readGAlignmentPairs(untreated3_chr4(), use.names=TRUE, param=param0)
head(U3.galp)
@

The \Rcode{show} method for \Rclass{GAlignmentPairs} objects displays
two {\tt ranges} columns, one for the {\it first} alignment in the pair (the
left column), and one for the {\it last} alignment in the pair (the right
column). The {\tt strand} column corresponds to the strand of the {\it first}
alignment.

<<first-and-last-U3.galp>>=
head(first(U3.galp))
head(last(U3.galp))
@

According to the SAM format specifications, the aligner is expected to mark
each alignment pair as {\it proper} or not (flag bit 0x2 in the SAM format).
The SAM Spec only says that a pair is {\it proper} if the {\it first} and
{\it 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 \Rcode{isProperPair} to extract this flag from the
\Rclass{GAlignmentPairs} object:

<<isProperPair>>=
table(isProperPair(U3.galp))
@

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

<<keep-only-proper-pairs>>=
U3.GALP <- U3.galp[isProperPair(U3.galp)]
@

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

<<U3.GALP_names_is_dup>>=
U3.GALP_names_is_dup <- duplicated(names(U3.GALP))
table(U3.GALP_names_is_dup)
@

Storing the {\it query template names} in a factor will be useful:

<<U3.GALP_qnames>>=
U3.uqnames <- unique(names(U3.GALP))
U3.GALP_qnames <- factor(names(U3.GALP), levels=U3.uqnames)
@

as well as having the mapping between each {\it query template name} and its
first occurence in \Rcode{U3.GALP\_qnames}:

<<U3.GALP_dup2unq>>=
U3.GALP_dup2unq <- match(U3.GALP_qnames, U3.GALP_qnames)
@

Our reads can have up to 1 gap per end:

<<gaps-in-U3.GALP>>=
head(unique(cigar(first(U3.GALP))))
head(unique(cigar(last(U3.GALP))))
table(ngap(first(U3.GALP)), ngap(last(U3.GALP)))
@

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

<<no-indels-in-U3.GALP>>=
colSums(cigarOpTable(cigar(first(U3.GALP))))
colSums(cigarOpTable(cigar(last(U3.GALP))))
@



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Find all the overlaps between the reads and transcripts}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Load the transcripts from a \Rclass{TranscriptDb} object}

In order to compute overlaps between reads and transcripts, we need access
to the genomic positions of a set of known transcripts and their exons.
It is essential that the reference genome of this set of transcripts and exons
be {\bf exactly} the same as the reference genome used to align the reads.

We could use the \Rfunction{makeTranscriptDbFromUCSC} function defined in
the \Rpackage{GenomicFeatures} package to make a \Rclass{TranscriptDb}
object containing the dm3 transcripts and their exons retrieved from the
UCSC Genome Browser\footnote{\url{http://genome.ucsc.edu/cgi-bin/hgGateway}}.
The Bioconductor project however provides a few annotation packages
containing \Rclass{TranscriptDb} objects for the most commonly studied
organisms (those data packages are sometimes called the {\it TxDb} packages).
One of them is the \Rpackage{TxDb.Dmelanogaster.\-UCSC.\-dm3.ensGene} package.
It contains a \Rclass{TranscriptDb} object that was made by pointing
the \Rfunction{makeTranscriptDbFromUCSC} function to the dm3 genome
and {\it Ensembl Genes} track \footnote{See
\url{http://genome.ucsc.edu/cgi-bin/hgTrackUi?hgsid=276880911&g=ensGene}
for a description of this track.}. We can use it here:

<<txdb>>=
library(TxDb.Dmelanogaster.UCSC.dm3.ensGene)
TxDb.Dmelanogaster.UCSC.dm3.ensGene
txdb <- TxDb.Dmelanogaster.UCSC.dm3.ensGene
@

We extract the exons grouped by transcript in a \Rclass{GRangesList} object:

<<exbytx>>=
exbytx <- exonsBy(txdb, by="tx", use.names=TRUE)
length(exbytx)  # nb of transcripts
@
<<CHECK_exbytx, echo=FALSE, results=hide>>=
.checkIdenticalToPrecomputed(exbytx, "exbytx", ignore.metadata=TRUE)
@

We check that all the exons in any given transcript belong to the same
chromosome and strand. Knowing that our set of transcripts is free of
this sort of trans-splicing events typically allows some significant
simplifications during the downstream analysis \footnote{Dealing with
trans-splicing events is not covered in this document.}.
A quick and easy way to check this is to take advantage of the fact
that \Rcode{seqnames} and \Rcode{strand} return \Rclass{RleList} objects.
So we can extract the number of Rle runs for each transcript and make
sure it's always 1:

<<check-for-trans-splicing-in-exbytx>>=
table(elementLengths(runLength(seqnames(exbytx))))
table(elementLengths(runLength(strand(exbytx))))
@

Therefore the strand of any given transcript is unambiguously defined
and can be extracted with:

<<exbytx_strand>>=
exbytx_strand <- unlist(runValue(strand(exbytx)), use.names=FALSE)
@

We will also need the mapping between the transcripts and their gene.
We start by using \Rfunction{transcripts} to extract this information
from our \Rclass{TranscriptDb} object \Rcode{txdb}, and
then we construct a named factor that represents the mapping:

<<exbytx2gene>>=
tx <- transcripts(txdb, columns=c("tx_name", "gene_id"))
head(tx)
df <- mcols(tx)
exbytx2gene <- as.character(df$gene_id)
exbytx2gene <- factor(exbytx2gene, levels=unique(exbytx2gene))
names(exbytx2gene) <- df$tx_name
exbytx2gene <- exbytx2gene[names(exbytx)]
head(exbytx2gene)
nlevels(exbytx2gene)  # nb of genes
@


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Single-end overlaps}

\subsubsection{Find the single-end overlaps}

We are ready to compute the overlaps with the \Rfunction{findOverlaps}
function. Note that the strand of the queries produced by the RNA-seq
experiment is typically unknown so we use \Rcode{ignore.strand=TRUE}:

<<U1.OV00>>=
U1.OV00 <- findOverlaps(U1.GAL, exbytx, ignore.strand=TRUE)
@

\Rcode{U1.OV00} is a \Rclass{Hits} object that contains 1 element per
overlap. Its length gives the number of overlaps:

<<length-of-U1.OV00>>=
length(U1.OV00)
@

\subsubsection{Tabulate the single-end overlaps}

We will repeatedly use the 2 following little helper functions to ``tabulate''
the overlaps in a given \Rclass{Hits} object (e.g. \Rcode{U1.OV00}), i.e. to
count the number of overlaps for each element in the query or for each
element in the subject:

<<nhitPerQuery-and-nhitPerSubject>>=
nhitPerQuery <- function(x) tabulate(queryHits(x), nbins=queryLength(x))
nhitPerSubject <- function(x) tabulate(subjectHits(x), nbins=subjectLength(x))
@

Number of transcripts for each alignment in \Rcode{U1.GAL}:

<<U1.GAL_ntx>>=
U1.GAL_ntx <- nhitPerQuery(U1.OV00)
mcols(U1.GAL)$ntx <- U1.GAL_ntx
head(U1.GAL)
table(U1.GAL_ntx)
mean(U1.GAL_ntx >= 1)
@

76\% of the alignments in \Rcode{U1.GAL} have an overlap with
at least 1 transcript in \Rcode{exbytx}.

Note that \Rfunction{countOverlaps} can be used directly on \Rcode{U1.GAL}
and \Rcode{exbytx} for computing \Rcode{U1.GAL\_ntx}:

<<U1.GAL_ntx_again, eval=FALSE>>=
U1.GAL_ntx_again <- countOverlaps(U1.GAL, exbytx, ignore.strand=TRUE)
stopifnot(identical(unname(U1.GAL_ntx_again), U1.GAL_ntx))
@

Because \Rcode{U1.GAL} can (and actually does) contain more than 1 alignment
per {\it original query} (aka read), we also count the number of transcripts
for each read:

<<U1.uqnames_ntx>>=
U1.OV10 <- remapHits(U1.OV00, query.map=U1.GAL_qnames)
U1.uqnames_ntx <- nhitPerQuery(U1.OV10)
names(U1.uqnames_ntx) <- U1.uqnames
table(U1.uqnames_ntx)
mean(U1.uqnames_ntx >= 1)
@

78.4\% of the reads have an overlap with
at least 1 transcript in \Rcode{exbytx}.

Number of reads for each transcript:

<<U1.exbytx_nOV10>>=
U1.exbytx_nOV10 <- nhitPerSubject(U1.OV10)
names(U1.exbytx_nOV10) <- names(exbytx)
mean(U1.exbytx_nOV10 >= 50)
@

Only 0.869\% of the transcripts in \Rcode{exbytx} have an overlap with
at least 50 reads.

Top 10 transcripts:

<<top-10-transcripts-based-on-U1.exbytx_nOV10>>=
head(sort(U1.exbytx_nOV10, decreasing=TRUE), n=10) 
@


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Paired-end overlaps}

\subsubsection{Find the paired-end overlaps}

Like with our single-end overlaps, we call \Rfunction{findOverlaps} with
\Rcode{ignore.strand=TRUE}:

<<U3.OV00>>=
U3.OV00 <- findOverlaps(U3.GALP, exbytx, ignore.strand=TRUE)
@

Like \Rcode{U1.OV00}, \Rcode{U3.OV00} is a \Rclass{Hits} object. Its length gives
the number of paired-end overlaps:

<<length-of-U3.OV00>>=
length(U3.OV00)
@

\subsubsection{Tabulate the paired-end overlaps}

Number of transcripts for each alignment pair in \Rcode{U3.GALP}:

<<U3.GALP_ntx>>=
U3.GALP_ntx <- nhitPerQuery(U3.OV00)
mcols(U3.GALP)$ntx <- U3.GALP_ntx
head(U3.GALP)
table(U3.GALP_ntx)
mean(U3.GALP_ntx >= 1)
@

71\% of the alignment pairs in \Rcode{U3.GALP} have an overlap
with at least 1 transcript in \Rcode{exbytx}.

Note that \Rfunction{countOverlaps} can be used directly on \Rcode{U3.GALP}
and \Rcode{exbytx} for computing \Rcode{U3.GALP\_ntx}:

<<U3.GALP_ntx_again, eval=FALSE>>=
U3.GALP_ntx_again <- countOverlaps(U3.GALP, exbytx, ignore.strand=TRUE)
stopifnot(identical(unname(U3.GALP_ntx_again), U3.GALP_ntx))
@

Because \Rcode{U3.GALP} can (and actually does) contain more than 1 alignment
pair per {\it original query template}, we also count the number of
transcripts for each template:

<<U3.uqnames_ntx>>=
U3.OV10 <- remapHits(U3.OV00, query.map=U3.GALP_qnames)
U3.uqnames_ntx <- nhitPerQuery(U3.OV10)
names(U3.uqnames_ntx) <- U3.uqnames
table(U3.uqnames_ntx)
mean(U3.uqnames_ntx >= 1)
@

72.3\% of the templates have an overlap with at least 1 transcript in
\Rcode{exbytx}.

Number of templates for each transcript:

<<U3.exbytx_nOV10>>=
U3.exbytx_nOV10 <- nhitPerSubject(U3.OV10)
names(U3.exbytx_nOV10) <- names(exbytx)
mean(U3.exbytx_nOV10 >= 50)
@

Only 0.756\% of the transcripts in \Rcode{exbytx} have an overlap with
at least 50 templates.

Top 10 transcripts:

<<top-10-transcripts-based-on-U3.exbytx_nOV10>>=
head(sort(U3.exbytx_nOV10, decreasing=TRUE), n=10)
@



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Encode the overlaps between the reads and transcripts}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Single-end encodings}

The {\it overlap encodings} are strand sensitive so we will compute them
twice, once for the ``original alignments'' (i.e. the alignments of the
{\it original queries}), and once again for the ``flipped alignments''
(i.e. the alignments of the ``flipped {\it original queries}''). We extract
the ranges of the ``original'' and ``flipped'' alignments in 2
\Rclass{GRangesList} objects with:

<<U1.grl_and_U1.grlf>>=
U1.grl <- grglist(U1.GAL, order.as.in.query=TRUE)
U1.grlf <- flipQuery(U1.grl)  # flipped
@

and encode their overlaps with the transcripts:

<<U1.ovencAB>>=
U1.ovencA <- encodeOverlaps(U1.grl, exbytx, hits=U1.OV00)
U1.ovencB <- encodeOverlaps(U1.grlf, exbytx, hits=U1.OV00)
@

\Rcode{U1.ovencA} and \Rcode{U1.ovencB} are 2 \Rclass{OverlapsEncodings}
objects of the same length as \Rclass{Hits} object \Rcode{U1.OV00}.
For each hit in \Rcode{U1.OV00}, we have 2 corresponding encodings, one
in \Rcode{U1.ovencA} and one in \Rcode{U1.ovencB}, but only one of them
encodes a hit between alignment ranges and exon ranges that are on the
same strand.
We use the \Rfunction{selectEncodingWithCompatibleStrand} function
to merge them into a single \Rclass{OverlapsEncodings} of the same
length. For each hit in \Rcode{U1.OV00}, this selects the encoding
corresponding to alignment ranges and exon ranges with compatible
strand:

<<U1.ovenc>>=
U1.grl_strand <- unlist(runValue(strand(U1.grl)), use.names=FALSE)
U1.ovenc <- selectEncodingWithCompatibleStrand(U1.ovencA, U1.ovencB,
                                               U1.grl_strand, exbytx_strand,
                                               hits=U1.OV00)
U1.ovenc
@

As a convenience, the 2 above calls to \Rfunction{encodeOverlaps} + merging
step can be replaced by a single call to \Rfunction{encodeOverlaps} on
\Rcode{U1.grl} (or \Rcode{U1.grlf}) with \Rcode{flip.query.if.wrong.strand=TRUE}:

<<U1.ovenc_again>>=
U1.ovenc_again <- encodeOverlaps(U1.grl, exbytx, hits=U1.OV00, flip.query.if.wrong.strand=TRUE)
stopifnot(identical(U1.ovenc_again, U1.ovenc))
@

Unique encodings in \Rcode{U1.ovenc}:

<<U1.ovenc_table>>=
U1.unique_encodings <- levels(U1.ovenc)
length(U1.unique_encodings)
head(U1.unique_encodings)
U1.ovenc_table <- table(encoding(U1.ovenc))
tail(sort(U1.ovenc_table))
@

Encodings are sort of cryptic but utilities are provided to extract
specific meaning from them. Use of these utilities is covered later in this
document.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Paired-end encodings}

Let's encode the overlaps in \Rcode{U3.OV00}:

<<U3.ovenc>>=
U3.grl <- grglist(U3.GALP, order.as.in.query=TRUE)
U3.ovenc <- encodeOverlaps(U3.grl, exbytx, hits=U3.OV00, flip.query.if.wrong.strand=TRUE)
U3.ovenc
@

Unique encodings in \Rcode{U3.ovenc}:

<<U3.ovenc_table>>=
U3.unique_encodings <- levels(U3.ovenc)
length(U3.unique_encodings)
head(U3.unique_encodings)
U3.ovenc_table <- table(encoding(U3.ovenc))
tail(sort(U3.ovenc_table))
@



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{``Compatible'' overlaps}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


We are interested in a particular type of overlap where the read
overlaps the transcript in a ``compatible'' way, that is, in a way
compatible with the splicing of the transcript.
The \Rfunction{isCompatibleWithSplicing} function can be used on an
\Rclass{OverlapEncodings} object to detect this type of overlap.
Note that \Rfunction{isCompatibleWithSplicing} can also be used on a
character vector or factor.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{``Compatible'' single-end overlaps}

\subsubsection{``Compatible'' single-end encodings}

\Rcode{U1.ovenc} contains 7 unique encodings ``compatible'' with the splicing
of the transcript:

<<U1-unique-compatible-encodings>>=
sort(U1.ovenc_table[isCompatibleWithSplicing(U1.unique_encodings)])
@

Encodings \Rcode{"1:i:"} (455176 occurences in \Rcode{U1.ovenc}),
\Rcode{"2:jm:af:"} (72929 occurences in \Rcode{U1.ovenc}),
and \Rcode{"3:jmm:agm:aaf:"} (488 occurences in \Rcode{U1.ovenc}),
correspond to the following overlaps:

\begin{itemize}
  \item \Rcode{"1:i:"}
\begin{verbatim}
  - read (no gap):            oooooooo
  - transcript:     ...   >>>>>>>>>>>>>>   ...
\end{verbatim}

  \item \Rcode{"2:jm:af:"}
\begin{verbatim}
  - read (1 gap):             ooooo---ooo
  - transcript:     ...   >>>>>>>>>   >>>>>>>>>   ...
\end{verbatim}

  \item \Rcode{"3:jmm:agm:aaf:"}
\begin{verbatim}
  - read (2 gaps):              oo---ooooo---o
  - transcript:     ...   >>>>>>>>   >>>>>   >>>>>>>   ...
\end{verbatim}
\end{itemize}

For clarity, only the exons involved in the overlap are represented. The
transcript can of course have more upstream and downstream exons, which
is denoted by the ... on the left side (5' end) and right side (3' end)
of each drawing. Note that the exons represented in the 2nd and 3rd drawings
are consecutive and adjacent in the processed transcript.

Encodings \Rcode{"1:f:"} and \Rcode{"1:j:"} are variations of the situation
described by encoding \Rcode{"1:i:"}. For \Rcode{"1:f:"}, the first
aligned base of the read (or ``flipped'' read) is aligned with the first
base of the exon. For \Rcode{"1:j:"}, the last aligned base of the read (or
``flipped'' read) is aligned with the last base of the exon:

\begin{itemize}
  \item \Rcode{"1:f:"}
\begin{verbatim}
  - read (no gap):        oooooooo
  - transcript:     ...   >>>>>>>>>>>>>>   ...
\end{verbatim}

  \item \Rcode{"1:j:"}
\begin{verbatim}
  - read (no gap):              oooooooo
  - transcript:     ...   >>>>>>>>>>>>>>   ...
\end{verbatim}
\end{itemize}

<<U1.OV00_is_comp>>=
U1.OV00_is_comp <- isCompatibleWithSplicing(U1.ovenc)
table(U1.OV00_is_comp)  # 531797 "compatible" overlaps
@

Finally, let's extract the ``compatible'' overlaps from \Rcode{U1.OV00}:

<<U1.compOV00>>=
U1.compOV00 <- U1.OV00[U1.OV00_is_comp]
@

Note that high-level convenience wrapper \Rfunction{findCompatibleOverlaps}
can be used for computing the ``compatible'' overlaps directly between
a \Rclass{GAlignments} object (containing reads) and a
\Rclass{GRangesList} object (containing transcripts):

<<U1.compOV00_again, eval=FALSE>>=
U1.compOV00_again <- findCompatibleOverlaps(U1.GAL, exbytx)
stopifnot(identical(U1.compOV00_again, U1.compOV00))
@

\subsubsection{Tabulate the ``compatible'' single-end overlaps}

Number of ``compatible'' transcripts for each alignment in \Rcode{U1.GAL}:

<<U1.GAL_ncomptx>>=
U1.GAL_ncomptx <- nhitPerQuery(U1.compOV00)
mcols(U1.GAL)$ncomptx <- U1.GAL_ncomptx
head(U1.GAL)
table(U1.GAL_ncomptx)
mean(U1.GAL_ncomptx >= 1)
@

75\% of the alignments in \Rcode{U1.GAL} are ``compatible''
with at least 1 transcript in \Rcode{exbytx}.

Note that high-level convenience wrapper \Rfunction{countCompatibleOverlaps}
can be used directly on \Rcode{U1.GAL} and \Rcode{exbytx} for computing
\Rcode{U1.GAL\_ncomptx}:

<<U1.GAL_ncomptx_again, eval=FALSE>>=
U1.GAL_ncomptx_again <- countCompatibleOverlaps(U1.GAL, exbytx)
stopifnot(identical(U1.GAL_ncomptx_again, U1.GAL_ncomptx))
@

Number of ``compatible'' transcripts for each read:

<<U1.uqnames_ncomptx>>=
U1.compOV10 <- remapHits(U1.compOV00, query.map=U1.GAL_qnames)
U1.uqnames_ncomptx <- nhitPerQuery(U1.compOV10)
names(U1.uqnames_ncomptx) <- U1.uqnames
table(U1.uqnames_ncomptx)
mean(U1.uqnames_ncomptx >= 1)
@

77.5\% of the reads are ``compatible'' with at least 1 transcript in
\Rcode{exbytx}.

Number of ``compatible'' reads for each transcript:

<<U1.exbytx_ncompOV10>>=
U1.exbytx_ncompOV10 <- nhitPerSubject(U1.compOV10)
names(U1.exbytx_ncompOV10) <- names(exbytx)
mean(U1.exbytx_ncompOV10 >= 50)
@

Only 0.87\% of the transcripts in \Rcode{exbytx} are ``compatible''
with at least 50 reads.

Top 10 transcripts:

<<top-10-transcripts-based-on-U1.exbytx_ncompOV10>>=
head(sort(U1.exbytx_ncompOV10, decreasing=TRUE), n=10)
@

Note that this ``top 10'' is slightly different from the ``top 10'' we
obtained earlier when we counted {\bf all} the overlaps.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{``Compatible'' paired-end overlaps}

\subsubsection{``Compatible'' paired-end encodings}

\Rcode{U3.ovenc} contains 13 unique paired-end encodings ``compatible'' with
the splicing of the transcript:

<<U3-unique-compatible-encodings>>=
sort(U3.ovenc_table[isCompatibleWithSplicing(U3.unique_encodings)])
@

Paired-end encodings \Rcode{"1{-}{-}1:i{-}{-}i:"} (100084 occurences in
\Rcode{U3.ovenc}), \Rcode{"2{-}{-}1:jm{-}{-}m:af{-}{-}i:"} (2700 occurences
in \Rcode{U3.ovenc}), \Rcode{"1{-}{-}2:i{-}{-}jm:a{-}{-}af:"} (2480 occurences
in \Rcode{U3.ovenc}), \Rcode{"1{-}{-}1:i{-}{-}m:a{-}{-}i:"} (287 occurences
in \Rcode{U3.ovenc}), and \Rcode{"2{-}{-}2:jm{-}{-}mm:af{-}{-}jm:aa{-}{-}af:"}
(153 occurences in \Rcode{U3.ovenc}), correspond to the following paired-end
overlaps:

\begin{itemize}
  \item \Rcode{"1{-}{-}1:i{-}{-}i:"}
\begin{verbatim}
  - paired-end read (no gap on the first end, no gap on the
    last end):               oooo   oooo
  - transcript:     ...   >>>>>>>>>>>>>>>>   ...
\end{verbatim}

  \item \Rcode{"2{-}{-}1:jm{-}{-}m:af{-}{-}i:"}
\begin{verbatim}
  - paired-end read (1 gap on the first end, no gap on the
    last end):                 ooo---o    oooo
  - transcript:     ...   >>>>>>>>   >>>>>>>>>>>   ...
\end{verbatim}

  \item \Rcode{"1{-}{-}2:i{-}{-}jm:a{-}{-}af:"}
\begin{verbatim}
  - paired-end read (no gap on the first end, 1 gap on the
    last end):             oooo      oo---oo
  - transcript:     ...  >>>>>>>>>>>>>>   >>>>>>>>>  ...
\end{verbatim}

  \item \Rcode{"1{-}{-}1:i{-}{-}m:a{-}{-}i:"}
\begin{verbatim}
  - paired-end read (no gap on the first end, no gap on the
    last end):               oooo      oooo
  - transcript:     ...   >>>>>>>>>   >>>>>>>   ...
\end{verbatim}

  \item \Rcode{"2{-}{-}2:jm{-}{-}mm:af{-}{-}jm:aa{-}{-}af:"}
\begin{verbatim}
  - paired-end read (1 gap on the first end, 1 gap on the
    last end):               ooo---o    oo---oo
  - transcript:     ...   >>>>>>   >>>>>>>   >>>>>   ...
\end{verbatim}
\end{itemize}

Note: switch use of ``first'' and ``last'' above if the read was ``flipped''.

<<U3.OV00_is_comp>>=
U3.OV00_is_comp <- isCompatibleWithSplicing(U3.ovenc)
table(U3.OV00_is_comp)  # 106835 "compatible" paired-end overlaps
@

Finally, let's extract the ``compatible'' paired-end overlaps from
\Rcode{U3.OV00}:

<<U3.compOV00>>=
U3.compOV00 <- U3.OV00[U3.OV00_is_comp]
@

Note that, like with our single-end reads, high-level convenience
wrapper \Rfunction{findCompatibleOverlaps} can be used for computing
the ``compatible'' paired-end overlaps directly between a
\Rclass{GAlignmentPairs} object (containing paired-end reads)
and a \Rclass{GRangesList} object (containing transcripts):

<<U3.compOV00_again, eval=FALSE>>=
U3.compOV00_again <- findCompatibleOverlaps(U3.GALP, exbytx)
stopifnot(identical(U3.compOV00_again, U3.compOV00))
@

\subsubsection{Tabulate the ``compatible'' paired-end overlaps}

Number of ``compatible'' transcripts for each alignment pair in
\Rcode{U3.GALP}:

<<U3.GALP_ncomptx>>=
U3.GALP_ncomptx <- nhitPerQuery(U3.compOV00)
mcols(U3.GALP)$ncomptx <- U3.GALP_ncomptx
head(U3.GALP)
table(U3.GALP_ncomptx)
mean(U3.GALP_ncomptx >= 1)
@

69.7\% of the alignment pairs in \Rcode{U3.GALP} are ``compatible''
with at least 1 transcript in \Rcode{exbytx}.

Note that high-level convenience wrapper \Rfunction{countCompatibleOverlaps}
can be used directly on \Rcode{U3.GALP} and \Rcode{exbytx} for computing
\Rcode{U3.GALP\_ncomptx}:

<<U3.GALP_ncomptx_again, eval=FALSE>>=
U3.GALP_ncomptx_again <- countCompatibleOverlaps(U3.GALP, exbytx)
stopifnot(identical(U3.GALP_ncomptx_again, U3.GALP_ncomptx))
@

Number of ``compatible'' transcripts for each template:

<<U3.uqnames_ncomptx>>=
U3.compOV10 <- remapHits(U3.compOV00, query.map=U3.GALP_qnames)
U3.uqnames_ncomptx <- nhitPerQuery(U3.compOV10)
names(U3.uqnames_ncomptx) <- U3.uqnames
table(U3.uqnames_ncomptx)
mean(U3.uqnames_ncomptx >= 1)
@

70.7\% of the templates are ``compatible'' with at least 1 transcript in
\Rcode{exbytx}.

Number of ``compatible'' templates for each transcript:

<<U3.exbytx_ncompOV10>>=
U3.exbytx_ncompOV10 <- nhitPerSubject(U3.compOV10)
names(U3.exbytx_ncompOV10) <- names(exbytx)
mean(U3.exbytx_ncompOV10 >= 50)
@

Only 0.7\% of the transcripts in \Rcode{exbytx} are ``compatible'' with at
least 50 templates.

Top 10 transcripts:

<<top-10-transcripts-based-on-U3.exbytx_ncompOV10>>=
head(sort(U3.exbytx_ncompOV10, decreasing=TRUE), n=10)
@

Note that this ``top 10'' is slightly different from the ``top 10'' we
obtained earlier when we counted {\bf all} the paired-end overlaps.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Compute the {\it reference query sequences} and project them on
         the transcriptome}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Compute the {\it reference query sequences}}

The {\it reference query sequences} are the query sequences {\bf after}
alignment, by opposition to the {\it original query sequences} (aka ``true''
or ``real'' query sequences) which are the query sequences {\bf before}
alignment.

The {\it reference query sequences} can easily be computed by extracting
the nucleotides mapped to each read from the reference genome. This of
course requires that we have access to the reference genome used by the
aligner. In Bioconductor, the full genome sequence for the dm3 assembly
is stored in the \Rpackage{BSgenome.Dmelanogaster.UCSC.dm3} data package
\footnote{See \url{http://bioconductor.org/packages/release/data/annotation/}
for the full list of annotation packages available in the current release
of Bioconductor.}:

<<Dmelanogaster>>=
library(BSgenome.Dmelanogaster.UCSC.dm3)
Dmelanogaster
@

To extract the portions of the reference genome corresponding to the ranges
in \Rcode{U1.grl}, we can use the \Rfunction{extractTranscriptSeqs}
function defined in the \Rpackage{GenomicFeatures} package:

<<U1-reference-query-sequences>>=
library(GenomicFeatures)
U1.GAL_rqseq <- extractTranscriptSeqs(Dmelanogaster, U1.grl)
head(U1.GAL_rqseq)
@

When reads are paired-end, we need to extract separately the ranges
corresponding to their {\it first} ends (aka {\it first} segments in BAM
jargon) and those corresponding to their {\it last} ends (aka {\it last}
segments in BAM jargon):

<<U3.grl_first-and-U3.grl_last>>=
U3.grl_first <- grglist(first(U3.GALP), order.as.in.query=TRUE)
U3.grl_last <- grglist(last(U3.GALP, invert.strand=TRUE), order.as.in.query=TRUE)
@

Then we extract the portions of the reference genome corresponding to the
ranges in \Rclass{GRangesList} objects \Rcode{U3.grl\_first} and
\Rcode{U3.grl\_last}:

<<U3-reference-query-sequences>>=
U3.GALP_rqseq1 <- extractTranscriptSeqs(Dmelanogaster, U3.grl_first)
U3.GALP_rqseq2 <- extractTranscriptSeqs(Dmelanogaster, U3.grl_last)
@


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Project the single-end alignments on the transcriptome}

The \Rfunction{extractQueryStartInTranscript} function computes for each
overlap the position of the {\it query start} in the transcript:

<<U1.OV00_qstart>>=
U1.OV00_qstart <- extractQueryStartInTranscript(U1.grl, exbytx,
                                                hits=U1.OV00, ovenc=U1.ovenc)
head(subset(U1.OV00_qstart, U1.OV00_is_comp))
@

\Rcode{U1.OV00\_qstart} is a data frame with 1 row per overlap and 3 columns:
\begin{enumerate}
    \item \Rcode{startInTranscript}: the 1-based start position of the
    read with respect to the transcript. Position 1 always corresponds to
    the first base on the 5' end of the transcript sequence.
    \item \Rcode{firstSpannedExonRank}: the rank of the first exon spanned
    by the read, that is, the rank of the exon found at position
    \Rcode{startInTranscript} in the transcript.
    \item \Rcode{startInFirstSpannedExon}: the 1-based start position of the
    read with respect to the first exon spanned by the read.
\end{enumerate}

Having this information allows us for example to compare the read and
transcript nucleotide sequences for each ``compatible'' overlap.
If we use the {\it reference query sequence} instead of the {\it original
query sequence} for this comparison, then it should match {\bf exactly}
the sequence found at the {\it query start} in the transcript.

Let's start by using \Rfunction{extractTranscriptSeqs} again
to extract the transcript sequences (aka transcriptome) from the dm3
reference genome:

<<txseq>>=
txseq <- extractTranscriptSeqs(Dmelanogaster, exbytx)
@

For each ``compatible'' overlap, the read sequence in \Rcode{U1.GAL\_rqseq}
must be an {\it exact} substring of the transcript sequence in
\Rcode{exbytx\_seq}:

<<U1.OV00_rqseq-vs-U1.OV00_txseq>>=
U1.OV00_rqseq <- U1.GAL_rqseq[queryHits(U1.OV00)]
U1.OV00_rqseq[flippedQuery(U1.ovenc)] <- reverseComplement(U1.OV00_rqseq[flippedQuery(U1.ovenc)])
U1.OV00_txseq <- txseq[subjectHits(U1.OV00)]
stopifnot(all(
    U1.OV00_rqseq[U1.OV00_is_comp] ==
        narrow(U1.OV00_txseq[U1.OV00_is_comp],
               start=U1.OV00_qstart$startInTranscript[U1.OV00_is_comp],
               width=width(U1.OV00_rqseq)[U1.OV00_is_comp])
))
@

Because of this relationship between the {\it reference query sequence}
and the transcript sequence of a ``compatible'' overlap, and because of
the relationship between the {\it original query sequences} and the
{\it reference query sequences}, then the edit distance reported in the
NM tag is actually the edit distance between the {\it original query} and
the transcript of a ``compatible'' overlap.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Project the paired-end alignments on the transcriptome}

For a paired-end read, the {\it query start} is the start of its ``left end''.

<<U3.OV00_Lqstart>>=
U3.OV00_Lqstart <- extractQueryStartInTranscript(U3.grl, exbytx,
                                                 hits=U3.OV00, ovenc=U3.ovenc)
head(subset(U3.OV00_Lqstart, U3.OV00_is_comp))
@

Note that \Rfunction{extractQueryStartInTranscript} can be called with
\Rcode{for.query.right.end=TRUE} if we want this information for the ``right
ends'' of the reads:

<<U3.OV00_Rqstart>>=
U3.OV00_Rqstart <- extractQueryStartInTranscript(U3.grl, exbytx,
                                                 hits=U3.OV00, ovenc=U3.ovenc,
                                                 for.query.right.end=TRUE)
head(subset(U3.OV00_Rqstart, U3.OV00_is_comp))
@

Like with single-end reads, having this information allows us for example
to compare the read and transcript nucleotide sequences for each
``compatible'' overlap. If we use the {\it reference query sequence} instead
of the {\it original query sequence} for this comparison, then it should
match {\bf exactly} the sequences of the ``left'' and ``right'' ends of the
read in the transcript.

Let's assign the ``left and right reference query sequences'' to each
overlap:

<<U3.OV00_Lrqseq_and_Rrqseq>>=
U3.OV00_Lrqseq <- U3.GALP_rqseq1[queryHits(U3.OV00)]
U3.OV00_Rrqseq <- U3.GALP_rqseq2[queryHits(U3.OV00)]
@

For the single-end reads, the sequence associated with a ``flipped query''
just needed to be ``reverse complemented''. For paired-end reads, we also
need to swap the 2 sequences in the pair:

<<U3.OV00_Lrqseq_and_Rrqseq>>=
flip_idx <- which(flippedQuery(U3.ovenc))
tmp <- U3.OV00_Lrqseq[flip_idx]
U3.OV00_Lrqseq[flip_idx] <- reverseComplement(U3.OV00_Rrqseq[flip_idx])
U3.OV00_Rrqseq[flip_idx] <- reverseComplement(tmp)
@

Let's assign the transcript sequence to each overlap:

<<U3.OV00_txseq>>=
U3.OV00_txseq <- txseq[subjectHits(U3.OV00)]
@

For each ``compatible'' overlap, we expect the ``left and right reference
query sequences'' of the read to be {\it exact} substrings of the transcript
sequence. Let's check the ``left reference query sequences'':

<<U3.OV00_Lrqseq-vs-U3.OV00_txseq>>=
stopifnot(all(
    U3.OV00_Lrqseq[U3.OV00_is_comp] ==
        narrow(U3.OV00_txseq[U3.OV00_is_comp],
               start=U3.OV00_Lqstart$startInTranscript[U3.OV00_is_comp],
               width=width(U3.OV00_Lrqseq)[U3.OV00_is_comp])
))
@

and the ``right reference query sequences'':

<<U3.OV00_Rrqseq-vs-U3.OV00_txseq>>=
stopifnot(all(
    U3.OV00_Rrqseq[U3.OV00_is_comp] ==
        narrow(U3.OV00_txseq[U3.OV00_is_comp],
               start=U3.OV00_Rqstart$startInTranscript[U3.OV00_is_comp],
               width=width(U3.OV00_Rrqseq)[U3.OV00_is_comp])
))
@



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Align the reads to the transcriptome}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Aligning the reads to the reference genome is not the most efficient nor
accurate way to count the number of ``compatible'' overlaps per {\it original
query}. Supporting junction reads (i.e. reads that align with at least 1 gap)
introduces a significant computational cost during the alignment process.
Then, as we've seen in the previous sections, each alignment produced by the
aligner needs to be broken into a set of ranges (based on its CIGAR) and
those ranges compared to the ranges of the exons grouped by transcript.

A more straightforward and accurate approach is to align the reads directly
to the transcriptome, and without allowing the typical gap that the aligner
needs to introduce when aligning a junction read to the reference genome.
With this approach, a ``hit'' between a read and a transcript is necessarily
compatible with the splicing of the transcript. In case of a ``hit'', we'll
say that the read and the transcript are ``string-based compatible'' (to
differentiate from our previous notion of ``compatible'' overlaps that we
will call ``encoding-based compatible'' from now on, unless the context is
clear).


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Align the single-end reads to the transcriptome}

\subsubsection{Find the ``hits''}

The single-end reads are in \Rcode{U1.oqseq}, the transcriptome is in
\Rcode{exbytx\_seq}.

Since indels were not allowed/supported during the alignment of the reads to
the reference genome, we don't need to allow/support them either for aligning
the reads to the transcriptome. Also since our goal is to find (and count)
``compatible'' overlaps between reads and transcripts, we don't need to keep
track of the details of the alignments between the reads and the transcripts.
Finally, since BAM file {\tt untreated1\_chr4.bam} is not the full output of
the aligner but the subset obtained by keeping only the alignments located
on chr4, we don't need to align \Rcode{U1.oqseq} to the full transcriptome,
but only to the subset of \Rcode{exbytx\_seq} made of the transcripts located
on chr4.

With those simplifications in mind, we write the following function that we
will use to find the ``hits'' between the reads and the transcriptome:

<<findSequenceHits>>=
### A wrapper to vwhichPDict() that supports IUPAC ambiguity codes in 'qseq'
### and 'txseq', and treats them as such.
findSequenceHits <- function(qseq, txseq, which.txseq=NULL, max.mismatch=0)
{
    .asHits <- function(x, pattern_length)
    {
        query_hits <- unlist(x)
        if (is.null(query_hits))
            query_hits <- integer(0)
        subject_hits <- rep.int(seq_len(length(x)), elementLengths(x))
        new("Hits", queryHits=query_hits, subjectHits=subject_hits,
                    queryLength=pattern_length, subjectLength=length(x))
    }

    .isHitInTranscriptBounds <- function(hits, qseq, txseq)
    {
        sapply(seq_len(length(hits)),
            function(i) {
                pattern <- qseq[[queryHits(hits)[i]]]
                subject <- txseq[[subjectHits(hits)[i]]]
                v <- matchPattern(pattern, subject,
                                  max.mismatch=max.mismatch, fixed=FALSE)
                any(1L <= start(v) & end(v) <= length(subject))
            })
    }

    if (!is.null(which.txseq)) {
        txseq0 <- txseq
        txseq <- txseq[which.txseq]
    }

    names(qseq) <- NULL
    other <- alphabetFrequency(qseq, baseOnly=TRUE)[ , "other"]
    is_clean <- other == 0L  # "clean" means "no IUPAC ambiguity code"

    ## Find hits for "clean" original queries.
    qseq0 <- qseq[is_clean]
    pdict0 <- PDict(qseq0, max.mismatch=max.mismatch)
    m0 <- vwhichPDict(pdict0, txseq,
                      max.mismatch=max.mismatch, fixed="pattern")
    hits0 <- .asHits(m0, length(qseq0))
    hits0@queryLength <- length(qseq)
    hits0@queryHits <- which(is_clean)[hits0@queryHits]

    ## Find hits for non "clean" original queries.
    qseq1 <- qseq[!is_clean]
    m1 <- vwhichPDict(qseq1, txseq,
                      max.mismatch=max.mismatch, fixed=FALSE)
    hits1 <- .asHits(m1, length(qseq1))
    hits1@queryLength <- length(qseq)
    hits1@queryHits <- which(!is_clean)[hits1@queryHits]

    ## Combine the hits.
    query_hits <- c(queryHits(hits0), queryHits(hits1))
    subject_hits <- c(subjectHits(hits0), subjectHits(hits1))

    if (!is.null(which.txseq)) {
        ## Remap the hits.
        txseq <- txseq0
        subject_hits <- which.txseq[subject_hits]
        hits0@subjectLength <- length(txseq)
    }

    ## Order the hits.
    oo <- IRanges:::orderIntegerPairs(query_hits, subject_hits)
    hits0@queryHits <- query_hits[oo]
    hits0@subjectHits <- subject_hits[oo]

    if (max.mismatch != 0L) {
        ## Keep only "in bounds" hits.
        is_in_bounds <- .isHitInTranscriptBounds(hits0, qseq, txseq)
        hits0 <- hits0[is_in_bounds]
    }
    hits0
}
@

Let's compute the index of the transcripts in \Rcode{exbytx\_seq} located on
chr4 (\Rfunction{findSequenceHits} will restrict the search to those
transcripts):

<<which.txseq, eval=FALSE>>=
chr4tx <- transcripts(txdb, vals=list(tx_chrom="chr4"))
chr4txnames <- mcols(chr4tx)$tx_name
which.txseq <- match(chr4txnames, names(txseq))
@

We know that the aligner tolerated up to 6 mismatches per read.
The 3 following commands find the ``hits'' for each {\it original
query}, then find the ``hits'' for each ``flipped {\it original query}'',
and finally merge all the ``hits'' (note that the 3 commands take about
1 hour to complete on a modern laptop):

<<U1.sbcompHITS, eval=FALSE>>=
U1.sbcompHITSa <- findSequenceHits(U1.oqseq, txseq,
                                   which.txseq=which.txseq, max.mismatch=6)
U1.sbcompHITSb <- findSequenceHits(reverseComplement(U1.oqseq), txseq,
                                   which.txseq=which.txseq, max.mismatch=6)
U1.sbcompHITS <- union(U1.sbcompHITSa, U1.sbcompHITSb)
@
<<LOAD_U1.sbcompHITS, echo=FALSE, results=hide>>=
U1.sbcompHITSa <- .loadPrecomputed("U1.sbcompHITSa")
U1.sbcompHITSb <- .loadPrecomputed("U1.sbcompHITSb")
U1.sbcompHITS <- union(U1.sbcompHITSa, U1.sbcompHITSb)
@

\subsubsection{Tabulate the ``hits''}

Number of ``string-based compatible'' transcripts for each read:

<<U1.uqnames_nsbcomptx>>=
U1.uqnames_nsbcomptx <- nhitPerQuery(U1.sbcompHITS)
names(U1.uqnames_nsbcomptx) <- U1.uqnames
table(U1.uqnames_nsbcomptx)
mean(U1.uqnames_nsbcomptx >= 1)
@

77.7\% of the reads are ``string-based compatible'' with at least 1
transcript in \Rcode{exbytx}.

Number of ``string-based compatible'' reads for each transcript:

<<U1.exbytx_nsbcompHITS>>=
U1.exbytx_nsbcompHITS <- nhitPerSubject(U1.sbcompHITS)
names(U1.exbytx_nsbcompHITS) <- names(exbytx)
mean(U1.exbytx_nsbcompHITS >= 50)
@

Only 0.865\% of the transcripts in \Rcode{exbytx} are ``string-based
compatible'' with at least 50 reads.

Top 10 transcripts:

<<top-10-transcripts-based-on-U1.exbytx_nsbcompHITS>>=
head(sort(U1.exbytx_nsbcompHITS, decreasing=TRUE), n=10)
@

\subsubsection{A closer look at the ``hits''}

[WORK IN PROGRESS, might be removed or replaced soon...]

Any ``encoding-based compatible'' overlap is of course ``string-based
compatible'':

<<encoding-based-compatible-implies-string-based-compatible>>=
stopifnot(length(setdiff(U1.compOV10, U1.sbcompHITS)) == 0)
@

but the reverse is not true:

<<string-based-compatible-does-NOT-imply-encoding-based-compatible>>=
length(setdiff(U1.sbcompHITS, U1.compOV10))
@

%To understand why the {\it overlap encodings} approach doesn't find all
%the ``string-based compatible'' hits, let's look at the second hit in
%\Rcode{setdiff(U1.sbcompHITS, U1.compOV10)}. This is a perfect hit between
%read SRR031728.4692406 and transcript 18924:
%
%<<>>=
%matchPattern(U1.oqseq[[6306]], txseq[[18924]])
%U1.GAL_idx <- which(U1.GAL_qnames == "SRR031728.4692406")
%U1.GAL[U1.GAL_idx]
%U1.GAL_idx %in% queryHits(U1.OV00)
%U1.GAL[12636]
%which(queryHits(U1.OV00) == 12636)
%U1.OV00[305]
%as.character(encoding(U1.ovenc)[305])
%@


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Align the paired-end reads to the transcriptome}

[COMING SOON...]



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{``Almost compatible'' overlaps}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


In many aspects, ``compatible'' overlaps can be seen as perfect.
We are now insterested in a less perfect type of overlap where the
read overlaps the transcript in a way that {\it would} be ``compatible''
if 1 or more exons were removed from the transcript. In that case we say
that the overlap is ``almost compatible'' with the transcript.
The \Rfunction{isCompatibleWithSkippedExons} function can be used on an
\Rclass{OverlapEncodings} object to detect this type of overlap.
Note that \Rfunction{isCompatibleWithSkippedExons} can also be used on a
character vector of factor.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{``Almost compatible'' single-end overlaps}

\subsubsection{``Almost compatible'' single-end encodings}

\Rcode{U1.ovenc} contains 7 unique encodings ``almost compatible'' with the
splicing of the transcript:

<<U1-unique-almost-compatible-encodings>>=
sort(U1.ovenc_table[isCompatibleWithSkippedExons(U1.unique_encodings)])
@

Encodings \Rcode{"2:jm:am:af:"} (1015 occurences in \Rcode{U1.ovenc}),
\Rcode{"2:jm:am:am:af:"} (144 occurences in \Rcode{U1.ovenc}),
and \Rcode{"3:jmm:agm:aam:aaf:"} (21 occurences in \Rcode{U1.ovenc}),
correspond to the following overlaps:

\begin{itemize}
  \item \Rcode{"2:jm:am:af:"}
\begin{verbatim}
  - read (1 gap):           ooooo----------ooo
  - transcript:     ...   >>>>>>>   >>>>   >>>>>>>>   ...
\end{verbatim}

  \item \Rcode{"2:jm:am:am:af:"}
\begin{verbatim}
  - read (1 gap):           ooooo------------------ooo
  - transcript:     ...   >>>>>>>   >>>>   >>>>>   >>>>>>>>   ...
\end{verbatim}

  \item \Rcode{"3:jmm:agm:aam:aaf:"}
\begin{verbatim}
  - read (2 gaps):             oo---oooo-----------oo
  - transcript:     ...   >>>>>>>   >>>>   >>>>>   >>>>>>>>   ...
\end{verbatim}
\end{itemize}

<<U1.OV00_is_acomp>>=
U1.OV00_is_acomp <- isCompatibleWithSkippedExons(U1.ovenc)
table(U1.OV00_is_acomp)  # 1202 "almost compatible" overlaps
@

Finally, let's extract the ``almost compatible'' overlaps from \Rcode{U1.OV00}:

<<U1.acompOV00>>=
U1.acompOV00 <- U1.OV00[U1.OV00_is_acomp]
@

\subsubsection{Tabulate the ``almost compatible'' single-end overlaps}

Number of ``almost compatible'' transcripts for each alignment in
\Rcode{U1.GAL}:

<<U1.GAL_nacomptx>>=
U1.GAL_nacomptx <- nhitPerQuery(U1.acompOV00)
mcols(U1.GAL)$nacomptx <- U1.GAL_nacomptx
head(U1.GAL)
table(U1.GAL_nacomptx)
mean(U1.GAL_nacomptx >= 1)
@

Only 0.27\% of the alignments in \Rcode{U1.GAL} are ``almost compatible''
with at least 1 transcript in \Rcode{exbytx}.

Number of ``almost compatible'' alignments for each transcript:

<<U1.exbytx_nacompOV00>>=
U1.exbytx_nacompOV00 <- nhitPerSubject(U1.acompOV00)
names(U1.exbytx_nacompOV00) <- names(exbytx)
table(U1.exbytx_nacompOV00)
mean(U1.exbytx_nacompOV00 >= 50)
@

Only 0.017\% of the transcripts in \Rcode{exbytx} are ``almost compatible''
with at least 50 alignments in \Rcode{U1.GAL}.

Finally note that the ``query start in transcript'' values returned by
\Rfunction{extractQueryStartInTranscript} are also defined for ``almost
compatible'' overlaps:

<<U1.OV00_qstart>>=
head(subset(U1.OV00_qstart, U1.OV00_is_acomp))
@


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{``Almost compatible'' paired-end overlaps}

\subsubsection{``Almost compatible'' paired-end encodings}

\Rcode{U3.ovenc} contains 5 unique paired-end encodings ``almost compatible''
with the splicing of the transcript:

<<U3-unique-almost-compatible-encodings>>=
sort(U3.ovenc_table[isCompatibleWithSkippedExons(U3.unique_encodings)])
@

Paired-end encodings \Rcode{"2{-}{-}1:jm{-}{-}m:am{-}{-}m:af{-}{-}i:"}
(73 occurences in \Rcode{U3.ovenc}),
\Rcode{"1{-}{-}2:i{-}{-}jm:a{-}{-}am:a{-}{-}af:"} (53 occurences in
\Rcode{U3.ovenc}), and
\Rcode{"2{-}{-}2:jm{-}{-}mm:am{-}{-}mm:af{-}{-}jm:aa{-}{-}af:"}
(9 occurences in \Rcode{U3.ovenc}), correspond to the following paired-end
overlaps:

\begin{itemize}
  \item \Rcode{"2{-}{-}1:jm{-}{-}m:am{-}{-}m:af{-}{-}i:"}
\begin{verbatim}
  - paired-end read (1 gap on the first end, no gap on the
    last end):              ooo----------o  oooo
  - transcript:     ...   >>>>>   >>>>   >>>>>>>>>   ...
\end{verbatim}

  \item \Rcode{"1{-}{-}2:i{-}{-}jm:a{-}{-}am:a{-}{-}af:"}
\begin{verbatim}
  - paired-end read (no gap on the first end, 1 gap on the
    last end):              oooo   oo---------oo
  - transcript:     ...   >>>>>>>>>>>   >>>   >>>>>>   ...
\end{verbatim}

  \item \Rcode{"2{-}{-}2:jm{-}{-}mm:am{-}{-}mm:af{-}{-}jm:aa{-}{-}af:"}
\begin{verbatim}
  - paired-end read (1 gap on the first end, 1 gap on the
    last end):                o----------ooo   oo---oo
  - transcript:     ...   >>>>>   >>>>   >>>>>>>>   >>>>>>   ...
\end{verbatim}
\end{itemize}

Note: switch use of ``first'' and ``last'' above if the read was ``flipped''.

<<U3.OV00_is_acomp>>=
U3.OV00_is_acomp <- isCompatibleWithSkippedExons(U3.ovenc)
table(U3.OV00_is_acomp)  # 141 "almost compatible" paired-end overlaps
@

Finally, let's extract the ``almost compatible'' paired-end overlaps from
\Rcode{U3.OV00}:

<<U3.acompOV00>>=
U3.acompOV00 <- U3.OV00[U3.OV00_is_acomp]
@

\subsubsection{Tabulate the ``almost compatible'' paired-end overlaps}

Number of ``almost compatible'' transcripts for each alignment pair in
\Rcode{U3.GALP}:

<<U3.GALP_nacomptx>>=
U3.GALP_nacomptx <- nhitPerQuery(U3.acompOV00)
mcols(U3.GALP)$nacomptx <- U3.GALP_nacomptx
head(U3.GALP)
table(U3.GALP_nacomptx)
mean(U3.GALP_nacomptx >= 1)
@

Only 0.2\% of the alignment pairs in \Rcode{U3.GALP} are ``almost compatible''
with at least 1 transcript in \Rcode{exbytx}.

Number of ``almost compatible'' alignment pairs for each transcript:

<<U3.exbytx_nacompOV00>>=
U3.exbytx_nacompOV00 <- nhitPerSubject(U3.acompOV00)
names(U3.exbytx_nacompOV00) <- names(exbytx)
table(U3.exbytx_nacompOV00)
mean(U3.exbytx_nacompOV00 >= 50)
@

Only 0.0034\% of the transcripts in \Rcode{exbytx} are ``almost compatible''
with at least 50 alignment pairs in \Rcode{U3.GALP}.

Finally note that the ``query start in transcript'' values returned by
\Rfunction{extractQueryStartInTranscript} are also defined for ``almost
compatible'' paired-end overlaps:

<<U3.OV00_Lqstart-and-U3.OV00_Rqstart>>=
head(subset(U3.OV00_Lqstart, U3.OV00_is_acomp))
head(subset(U3.OV00_Rqstart, U3.OV00_is_acomp))
@



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Detect novel splice junctions}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{By looking at single-end overlaps}

An alignment in \Rcode{U1.GAL} with ``almost compatible'' overlaps but
no ``compatible'' overlaps suggests the presence of one or more transcripts
that are not in our annotations.

First we extract the index of those alignments ({\it nsj} here stands
for ``{\bf n}ovel {\bf s}plice {\bf j}unction''):

<<U1.GAL_is_nsj>>=
U1.GAL_is_nsj <- U1.GAL_nacomptx != 0L & U1.GAL_ncomptx == 0L
head(which(U1.GAL_is_nsj))
@

We make this an index into \Rcode{U1.OV00}:

<<U1.OV00_is_nsj>>=
U1.OV00_is_nsj <- queryHits(U1.OV00) %in% which(U1.GAL_is_nsj)
@

We intersect with \Rcode{U1.OV00\_is\_acomp} and then subset \Rcode{U1.OV00}
to keep only the overlaps that suggest novel splicing:

<<narrow-U1.OV00_is_nsj>>=
U1.OV00_is_nsj <- U1.OV00_is_nsj & U1.OV00_is_acomp
U1.nsjOV00 <- U1.OV00[U1.OV00_is_nsj]
@

For each overlap in \Rcode{U1.nsjOV00}, we extract the ranks of the
skipped exons (we use a list for this as there might be more than 1 skipped
exon per overlap):

<<U1.nsjOV00_skippedex>>=
U1.nsjOV00_skippedex <- extractSkippedExonRanks(U1.ovenc)[U1.OV00_is_nsj]
names(U1.nsjOV00_skippedex) <- queryHits(U1.nsjOV00)
table(elementLengths(U1.nsjOV00_skippedex))
@

Finally, we split \Rcode{U1.nsjOV00\_skippedex} by transcript names:

<<U1.exbytx_skippedex>>=
f <- factor(names(exbytx)[subjectHits(U1.nsjOV00)], levels=names(exbytx))
U1.exbytx_skippedex <- split(U1.nsjOV00_skippedex, f)
@

\Rcode{U1.exbytx\_skippedex} is a named list of named lists of integer
vectors. The first level of names (outer names) are transcript names and
the second level of names (inner names) are alignment indices into
\Rcode{U1.GAL}:

<<names-of-U1.exbytx_skippedex>>=
head(names(U1.exbytx_skippedex))  # transcript names
@

Transcript FBtr0089124 receives 7 hits. All of them skip exons 9 and 10:

<<FBtr0089124-skipped-exons>>=
U1.exbytx_skippedex$FBtr0089124
@

Transcript FBtr0089147 receives 4 hits. Two of them skip exon 2, one of them
skips exons 2 to 6, and one of them skips exon 10:

<<FBtr0089147-skipped-exons>>=
U1.exbytx_skippedex$FBtr0089147
@

A few words about the interpretation of \Rcode{U1.exbytx\_skippedex}:
Because of how we've conducted this analysis, the aligments reported in
\Rcode{U1.exbytx\_skippedex} are guaranteed to not have any ``compatible''
overlaps with other known transcripts. All we can say, for example in the
case of transcript FBtr0089124, is that the 7 reported hits that skip exons 9
and 10 show evidence of one or more unknown transcripts with a splice junction
that corresponds to the gap between exons 8 and 11. But without further
analysis, we can't make any assumption about the exons structure of those
unknown transcripts. In particular, we cannot assume the existence of an
unknown transcript made of the same exons as transcript FBtr0089124 minus
exons 9 and 10!


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{By looking at paired-end overlaps}

[COMING SOON...]



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{\Rcode{sessionInfo()}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

<<sessionInfo>>=
sessionInfo()
@

\end{document}