man/txpath-methods.Rd
7ff5c4a7
 \name{txpath-methods}
 
 \alias{txpath-methods}
 
 \alias{txpath}
 \alias{txpath,GRangesList-method}
 \alias{txpath,SplicingGraphs-method}
 
475e1010
 \alias{txweight}
 \alias{txweight,SplicingGraphs-method}
 \alias{txweight<-}
 \alias{txweight<-,SplicingGraphs-method}
7ff5c4a7
 
 
 \title{
   Extract the transcript paths of a splicing graph
 }
 
 \description{
   \code{txpath} extracts the transcript paths of the splicing graph
   of a given gene from a \link{SplicingGraphs} object.
 }
 
 \usage{
 txpath(x, as.matrix=FALSE)
 
 ## Related utility:
475e1010
 txweight(x)
 txweight(x) <- value
7ff5c4a7
 }
 
 \arguments{
   \item{x}{
     A \link{SplicingGraphs} object of length 1
475e1010
     or a \link[GenomicRanges]{GRangesList} object for \code{txpath}.
 
     A \link{SplicingGraphs} object for \code{txweight}.
7ff5c4a7
   }
   \item{as.matrix}{
     TODO
   }
475e1010
   \item{value}{
     A numeric vector containing the weights to assign to each
     transcript in \code{x}.
   }
7ff5c4a7
 }
 
 \details{
   TODO
 }
 
 \value{
   A named list-like object with one list element per transcript in the gene.
   Each list element is an integer vector that describes the \emph{path}
   of the transcript i.e. the \emph{Splicing Site ids} that it goes thru.
 }
 
 \author{
   H. Pages
 }
 
 \seealso{
7bc34093
   This man page is part of the \pkg{SplicingGraphs} package.
20221b45
   Please see \code{?`\link{SplicingGraphs-package}`} for an overview of the
   package and for an index of its man pages.
7ff5c4a7
 
20221b45
   Other topics related to this man page and documented in other packages:
   \itemize{
05644e55
     \item The \link[GenomicRanges]{GRangesList} class defined in the
           \pkg{GenomicRanges} package.
 
     \item The \link[GenomicAlignments]{GAlignments} and
           \link[GenomicAlignments]{GAlignmentPairs} classes
           defined in the \pkg{GenomicAlignments} package.
7ff5c4a7
 
ab342602
     \item \link[GenomicRanges]{findOverlaps-methods} in the
7bc34093
           \pkg{GenomicRanges} package.
7ff5c4a7
 
ab342602
     \item \link[GenomicAlignments]{encodeOverlaps-methods} in the
           \pkg{GenomicAlignments} package.
 
05644e55
     \item The \code{\link[Rsamtools]{ScanBamParam}} function defined in the
7bc34093
           \pkg{Rsamtools} package.
7ff5c4a7
   }
 }
 
 \examples{
475e1010
 ## ---------------------------------------------------------------------
 ## 1. Make SplicingGraphs object 'sg' from toy gene model (see
 ##    '?SplicingGraphs')
 ## ---------------------------------------------------------------------
 example(SplicingGraphs)
7ff5c4a7
 sg
 
 ## 'sg' has 1 element per gene and 'names(sg)' gives the gene ids.
 names(sg)
 
475e1010
 ## ---------------------------------------------------------------------
 ## 2. txpath()
 ## ---------------------------------------------------------------------
7ff5c4a7
 ## Note that the list elements in the returned IntegerList object
 ## always consist of an even number of Splicing Site ids in ascending
 ## order.
 txpath(sg["geneB"])
 txpath(sg["geneD"])
 strand(sg)
 
 txpath(sg["geneD"], as.matrix=TRUE)  # splicing matrix
475e1010
 
 ## ---------------------------------------------------------------------
 ## 3. txweight()
 ## ---------------------------------------------------------------------
 txweight(sg)
 plot(sg["geneD"])
 
 txweight(sg) <- 5
 txweight(sg)
 plot(sg["geneD"])  # FIXME: Edges not rendered with correct width!
 plot(sgraph(sg["geneD"], as.igraph=TRUE))  # need to use this for now
 
 txweight(sg)[8:11] <- 5 * (4:1)
 txweight(sg)
 plot(sgraph(sg["geneD"], tx_id.as.edge.label=TRUE, as.igraph=TRUE))
 
 ## ---------------------------------------------------------------------
 ## 4. An advanced example
 ## ---------------------------------------------------------------------
e1bddf5a
 ## [TODO: Counting "unambiguous compatible hits" per transcript should be
 ##  supported by countReads(). Simplify the code below when countReads()
 ##  supports this.]
475e1010
 ## Here we show how to find "unambiguous compatible hits" between a set
 ## of RNA-seq reads and a set of transcripts, that is, hits that are
 ## compatible with the splicing of exactly 1 transcript. Then we set the
 ## transcript weights based on the number of unambiguous compatible hits
 ## they received and we finally plot some splicing graphs that show
 ## the weighted transcripts.
 ## Note that the code we use to find the unambiguous compatible hits
 ## uses findOverlaps() and encodeOverlaps() defined in the IRanges and
 ## GenomicRanges packages. It only requires that the transcripts are
3ac6c8f8
 ## represented as a GRangesList object and the reads as a GAlignments
 ## (single-end) or GAlignmentPairs (paired-end) object, and therefore is
 ## not specific to SplicingGraphs.
475e1010
 
 ## First we load toy reads (single-end) from a BAM file. We filter out
 ## secondary alignments, reads not passing quality controls, and PCR or
 ## optical duplicates (see ?scanBamFlag in the Rsamtools package for
 ## more information):
d0a37504
 flag0 <- scanBamFlag(isSecondaryAlignment=FALSE,
475e1010
                      isNotPassingQualityControls=FALSE,
                      isDuplicate=FALSE)
 param0 <- ScanBamParam(flag=flag0)
3ac6c8f8
 gal <- readGAlignments(toy_reads_bam(), use.names=TRUE, param=param0)
475e1010
 gal
 
 ## Put the reads in a GRangesList object:
 grl <- grglist(gal, order.as.in.query=TRUE)
 
 ## Put the transcripts in a GRangesList object (made of exons grouped
 ## by transcript):
 ex_by_tx <- unlist(sg)
 
 ## Most of the times the RNA-seq protocol is unstranded so the strand
 ## reported in the BAM file (and propagated to 'grl') for each alignment
 ## is meaningless. Thus we should call findOverlaps() with
 ## 'ignore.strand=TRUE':
 ov0 <- findOverlaps(grl, ex_by_tx, ignore.strand=TRUE)
 
 ## Encode the overlaps (this compare the fragmentation of the reads with
 ## the splicing of the transcripts):
ee2387b8
 ovenc0 <- encodeOverlaps(grl, ex_by_tx, hits=ov0,
                          flip.query.if.wrong.strand=TRUE)
 ov0_is_compat <- isCompatibleWithSplicing(ovenc0)
475e1010
 
 ## Keep compatible overlaps only:
ee2387b8
 ov1 <- ov0[ov0_is_compat]
475e1010
 
 ## Only keep overlaps that are compatible with exactly 1 transcript:
 ov2 <- ov1[queryHits(ov1) \%in\% which(countQueryHits(ov1) == 1L)]
 nhit_per_tx <- countSubjectHits(ov2)
 names(nhit_per_tx) <- names(txweight(sg))
 nhit_per_tx
 
 txweight(sg) <- 2 * nhit_per_tx
 plot(sgraph(sg["geneA"], tx_id.as.edge.label=TRUE, as.igraph=TRUE))
 plot(sgraph(sg["geneB"], tx_id.as.edge.label=TRUE, as.igraph=TRUE))
7ff5c4a7
 }