Last data update: 2014.03.03

 R: Plot a set of transcripts along genomic coordinates.
plotTranscripts-methodsR Documentation

Plot a set of transcripts along genomic coordinates.

Description

plotTranscripts uses the Gviz package to plot the exon structure of a set of transcripts along genomic coordinates.

Usage

plotTranscripts(x, reads=NULL, from=NA, to=NA, max.plot.reads=200)

Arguments

x

A GRangesList object containing the genomic ranges of a set of exons grouped by transcript. Alternatively, x can be a TxDb object, or a SplicingGraphs object of length 1.

reads

A GAlignments or GAlignmentPairs object containing single-end or paired-end reads.

from, to

Single numeric values, giving the range of genomic coordinates to plot the tracks in. By default (i.e. from=NA and to=NA), the plot covers the range spanned by the transcripts. If from=NULL and to=NULL, then the plot covers the range spanned by the transcripts and the reads.

max.plot.reads

The maximum number of reads that will be plotted. When the number of reads that fall in the region being plotted is very large, plotting them all would take a long time and result in a plot that is not very useful. If that number is greater than max.plot.reads, then only max.plot.reads randomly chosen reads are plotted.

Author(s)

H. Pages

See Also

This man page is part of the SplicingGraphs package. Please see ?`SplicingGraphs-package` for an overview of the package and for an index of its man pages.

Other topics related to this man page and documented in other packages:

  • plotTranscripts is based on the plotTracks function defined in the Gviz package.

  • The GRangesList class defined in the GenomicRanges package.

  • The GAlignments and GAlignmentPairs classes defined in the GenomicAlignments package.

  • The TxDb class defined in the GenomicFeatures package.

Examples

## ---------------------------------------------------------------------
## A. PLOT TRANSCRIPTS
## ---------------------------------------------------------------------
example(SplicingGraphs)  # create SplicingGraphs object 'sg'
sg

## 'sg' has 1 element per gene and 'names(sg)' gives the gene ids.
names(sg)

## The transcripts of a given gene can be extracted with [[. The result
## is an *unnamed* GRangesList object containing the exons grouped by
## transcript:
sg[["geneD"]]
plotTranscripts(sg[["geneD"]])  # requires the Gviz package

## The transcripts of all the genes can be extracted with unlist(). The
## result is a *named* GRangesList object containing the exons grouped
## by transcript. The names on the object are the gene ids:
ex_by_tx <- unlist(sg)
ex_by_tx
plotTranscripts(ex_by_tx)

## ---------------------------------------------------------------------
## B. PLOT TRANSCRIPTS AND READS
## ---------------------------------------------------------------------
gal <- readGAlignments(toy_reads_bam(), use.names=TRUE)
plotTranscripts(sg[["geneA"]], reads=gal)
plotTranscripts(ex_by_tx, reads=gal)
plotTranscripts(ex_by_tx, reads=gal, from=1, to=320)
plotTranscripts(ex_by_tx, reads=gal[21:26], from=1, to=320)

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(SplicingGraphs)
Loading required package: GenomicFeatures
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, cbind, colnames, do.call, duplicated, eval, evalq,
    get, grep, grepl, intersect, is.unsorted, lapply, lengths, mapply,
    match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, rank,
    rbind, rownames, sapply, setdiff, sort, table, tapply, union,
    unique, unsplit

Loading required package: S4Vectors
Loading required package: stats4

Attaching package: 'S4Vectors'

The following objects are masked from 'package:base':

    colMeans, colSums, expand.grid, rowMeans, rowSums

Loading required package: IRanges
Loading required package: GenomeInfoDb
Loading required package: GenomicRanges
Loading required package: AnnotationDbi
Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.

Loading required package: GenomicAlignments
Loading required package: SummarizedExperiment
Loading required package: Biostrings
Loading required package: XVector
Loading required package: Rsamtools
Loading required package: Rgraphviz
Loading required package: graph

Attaching package: 'graph'

The following object is masked from 'package:Biostrings':

    complement

Loading required package: grid

Attaching package: 'Rgraphviz'

The following objects are masked from 'package:IRanges':

    from, to

The following objects are masked from 'package:S4Vectors':

    from, to

Warning messages:
1: replacing previous import 'IRanges::from' by 'Rgraphviz::from' when loading 'SplicingGraphs' 
2: replacing previous import 'IRanges::to' by 'Rgraphviz::to' when loading 'SplicingGraphs' 
> png(filename="/home/ddbj/snapshot/RGM3/R_BC/result/SplicingGraphs/plotTranscripts-methods.Rd_%03d_medium.png", width=480, height=480)
> ### Name: plotTranscripts-methods
> ### Title: Plot a set of transcripts along genomic coordinates.
> ### Aliases: plotTranscripts-methods plotTranscripts
> ###   plotTranscripts,GRangesList-method plotTranscripts,TxDb-method
> ###   plotTranscripts,SplicingGraphs-method
> 
> ### ** Examples
> 
> ## ---------------------------------------------------------------------
> ## A. PLOT TRANSCRIPTS
> ## ---------------------------------------------------------------------
> example(SplicingGraphs)  # create SplicingGraphs object 'sg'

SplcnG> ## ---------------------------------------------------------------------
SplcnG> ## 1. Load a toy gene model as a TxDb object
SplcnG> ## ---------------------------------------------------------------------
SplcnG> 
SplcnG> library(GenomicFeatures)

SplcnG> suppressWarnings(
SplcnG+   toy_genes_txdb <- makeTxDbFromGFF(toy_genes_gff())
SplcnG+ )
Import genomic features from the file as a GRanges object ... OK
Prepare the 'metadata' data frame ... OK
Make the TxDb object ... OK

SplcnG> ## ---------------------------------------------------------------------
SplcnG> ## 2. Compute all the splicing graphs (1 graph per gene) and return them
SplcnG> ##    in a SplicingGraphs object
SplcnG> ## ---------------------------------------------------------------------
SplcnG> 
SplcnG> ## Extract the exons grouped by transcript:
SplcnG> ex_by_tx <- exonsBy(toy_genes_txdb, by="tx", use.names=TRUE)

SplcnG> ## Extract the transcripts grouped by gene:
SplcnG> tx_by_gn <- transcriptsBy(toy_genes_txdb, by="gene")

SplcnG> sg <- SplicingGraphs(ex_by_tx, tx_by_gn)

SplcnG> sg
SplicingGraphs object with 5 gene(s) and 13 transcript(s)

SplcnG> ## Alternatively 'sg' can be constructed directly from the TxDb
SplcnG> ## object:
SplcnG> sg2 <- SplicingGraphs(toy_genes_txdb)  # same as 'sg'

SplcnG> sg2
SplicingGraphs object with 5 gene(s) and 13 transcript(s)

SplcnG> ## Note that because SplicingGraphs objects have a slot that is an
SplcnG> ## environment (for caching the bubbles), they cannot be compared with
SplcnG> ## 'identical()' (will always return FALSE). 'all.equal()' should be
SplcnG> ## used instead:
SplcnG> stopifnot(isTRUE(all.equal(sg2, sg)))

SplcnG> ## 'sg' has 1 element per gene and 'names(sg)' gives the gene ids:
SplcnG> length(sg)
[1] 5

SplcnG> names(sg)
[1] "geneA" "geneB" "geneC" "geneD" "geneE"

SplcnG> ## ---------------------------------------------------------------------
SplcnG> ## 3. Basic manipulation of a SplicingGraphs object
SplcnG> ## ---------------------------------------------------------------------
SplcnG> 
SplcnG> ## Basic accessors:
SplcnG> seqnames(sg)
geneA geneB geneC geneD geneE 
 chrX  chrX  chrX  chrX  chrX 
Levels: chrX

SplcnG> strand(sg)
geneA geneB geneC geneD geneE 
    +     -     +     +     + 
Levels: + - *

SplcnG> seqinfo(sg)
Seqinfo object with 1 sequence from an unspecified genome; no seqlengths:
  seqnames seqlengths isCircular genome
  chrX             NA         NA   <NA>

SplcnG> ## Number of transcripts per gene:
SplcnG> elementNROWS(sg)
geneA geneB geneC geneD geneE 
    2     2     3     4     2 

SplcnG> ## The transcripts of a given gene can be extracted with [[. The result
SplcnG> ## is an *unnamed* GRangesList object containing the exons grouped by
SplcnG> ## transcript:
SplcnG> sg[["geneD"]]
GRangesList object of length 4:
[[1]] 
GRanges object with 2 ranges and 5 metadata columns:
      seqnames     ranges strand |   exon_id   exon_name exon_rank start_SSid
         <Rle>  <IRanges>  <Rle> | <integer> <character> <integer>  <integer>
  [1]     chrX [601, 630]      + |        10         Dx2         1          1
  [2]     chrX [666, 675]      + |        12         Dx4         2          5
       end_SSid
      <integer>
  [1]         3
  [2]         6

[[2]] 
GRanges object with 2 ranges and 5 metadata columns:
      seqnames     ranges strand | exon_id exon_name exon_rank start_SSid
  [1]     chrX [601, 620]      + |       9       Dx1         1          1
  [2]     chrX [651, 700]      + |      11       Dx3         2          4
      end_SSid
  [1]        2
  [2]        8

[[3]] 
GRanges object with 3 ranges and 5 metadata columns:
      seqnames     ranges strand | exon_id exon_name exon_rank start_SSid
  [1]     chrX [601, 620]      + |       9       Dx1         1          1
  [2]     chrX [666, 675]      + |      12       Dx4         2          5
  [3]     chrX [691, 700]      + |      13       Dx5         3          7
      end_SSid
  [1]        2
  [2]        6
  [3]        8

...
<1 more element>
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths

SplcnG> ## See '?plotTranscripts' for how to plot those transcripts.
SplcnG> 
SplcnG> ## The transcripts of all the genes can be extracted with unlist(). The
SplcnG> ## result is a *named* GRangesList object containing the exons grouped
SplcnG> ## by transcript. The names on the object are the gene ids:
SplcnG> ex_by_tx <- unlist(sg)

SplcnG> ex_by_tx
GRangesList object of length 13:
$geneA 
GRanges object with 1 range and 5 metadata columns:
      seqnames    ranges strand |   exon_id   exon_name exon_rank start_SSid
         <Rle> <IRanges>  <Rle> | <integer> <character> <integer>  <integer>
  [1]     chrX  [11, 50]      + |         2         Ax2         1          1
       end_SSid
      <integer>
  [1]         3

$geneA 
GRanges object with 2 ranges and 5 metadata columns:
      seqnames    ranges strand | exon_id exon_name exon_rank start_SSid
  [1]     chrX [11,  40]      + |       1       Ax1         1          1
  [2]     chrX [71, 100]      + |       3       Ax3         2          4
      end_SSid
  [1]        2
  [2]        5

$geneB 
GRanges object with 2 ranges and 5 metadata columns:
      seqnames     ranges strand | exon_id exon_name exon_rank start_SSid
  [1]     chrX [251, 300]      - |      23       Bx1         1          3
  [2]     chrX [201, 230]      - |      20       Bx2         2          6
      end_SSid
  [1]        1
  [2]        4

...
<10 more elements>
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
> sg
SplicingGraphs object with 5 gene(s) and 13 transcript(s)
> 
> ## 'sg' has 1 element per gene and 'names(sg)' gives the gene ids.
> names(sg)
[1] "geneA" "geneB" "geneC" "geneD" "geneE"
> 
> ## The transcripts of a given gene can be extracted with [[. The result
> ## is an *unnamed* GRangesList object containing the exons grouped by
> ## transcript:
> sg[["geneD"]]
GRangesList object of length 4:
[[1]] 
GRanges object with 2 ranges and 5 metadata columns:
      seqnames     ranges strand |   exon_id   exon_name exon_rank start_SSid
         <Rle>  <IRanges>  <Rle> | <integer> <character> <integer>  <integer>
  [1]     chrX [601, 630]      + |        10         Dx2         1          1
  [2]     chrX [666, 675]      + |        12         Dx4         2          5
       end_SSid
      <integer>
  [1]         3
  [2]         6

[[2]] 
GRanges object with 2 ranges and 5 metadata columns:
      seqnames     ranges strand | exon_id exon_name exon_rank start_SSid
  [1]     chrX [601, 620]      + |       9       Dx1         1          1
  [2]     chrX [651, 700]      + |      11       Dx3         2          4
      end_SSid
  [1]        2
  [2]        8

[[3]] 
GRanges object with 3 ranges and 5 metadata columns:
      seqnames     ranges strand | exon_id exon_name exon_rank start_SSid
  [1]     chrX [601, 620]      + |       9       Dx1         1          1
  [2]     chrX [666, 675]      + |      12       Dx4         2          5
  [3]     chrX [691, 700]      + |      13       Dx5         3          7
      end_SSid
  [1]        2
  [2]        6
  [3]        8

...
<1 more element>
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
> plotTranscripts(sg[["geneD"]])  # requires the Gviz package
  - plotting genomic range: from=591.1, to=709.9
> 
> ## The transcripts of all the genes can be extracted with unlist(). The
> ## result is a *named* GRangesList object containing the exons grouped
> ## by transcript. The names on the object are the gene ids:
> ex_by_tx <- unlist(sg)
> ex_by_tx
GRangesList object of length 13:
$geneA 
GRanges object with 1 range and 5 metadata columns:
      seqnames    ranges strand |   exon_id   exon_name exon_rank start_SSid
         <Rle> <IRanges>  <Rle> | <integer> <character> <integer>  <integer>
  [1]     chrX  [11, 50]      + |         2         Ax2         1          1
       end_SSid
      <integer>
  [1]         3

$geneA 
GRanges object with 2 ranges and 5 metadata columns:
      seqnames    ranges strand | exon_id exon_name exon_rank start_SSid
  [1]     chrX [11,  40]      + |       1       Ax1         1          1
  [2]     chrX [71, 100]      + |       3       Ax3         2          4
      end_SSid
  [1]        2
  [2]        5

$geneB 
GRanges object with 2 ranges and 5 metadata columns:
      seqnames     ranges strand | exon_id exon_name exon_rank start_SSid
  [1]     chrX [251, 300]      - |      23       Bx1         1          3
  [2]     chrX [201, 230]      - |      20       Bx2         2          6
      end_SSid
  [1]        1
  [2]        4

...
<10 more elements>
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
> plotTranscripts(ex_by_tx)
  - plotting genomic range: from=-77.9, to=988.9
> 
> ## ---------------------------------------------------------------------
> ## B. PLOT TRANSCRIPTS AND READS
> ## ---------------------------------------------------------------------
> gal <- readGAlignments(toy_reads_bam(), use.names=TRUE)
> plotTranscripts(sg[["geneA"]], reads=gal)
  - plotting genomic range: from=2.1, to=108.9
  - nb of reads to plot (i.e. overlapping with that range): 11
Warning message:
'compare' is deprecated.
Use 'pcompare' instead.
See help("Deprecated") 
> plotTranscripts(ex_by_tx, reads=gal)
  - plotting genomic range: from=-77.9, to=988.9
  - nb of reads to plot (i.e. overlapping with that range): 42
Warning message:
'compare' is deprecated.
Use 'pcompare' instead.
See help("Deprecated") 
> plotTranscripts(ex_by_tx, reads=gal, from=1, to=320)
  - nb of reads to plot (i.e. overlapping with that range): 42
Warning message:
'compare' is deprecated.
Use 'pcompare' instead.
See help("Deprecated") 
> plotTranscripts(ex_by_tx, reads=gal[21:26], from=1, to=320)
  - nb of reads to plot (i.e. overlapping with that range): 6
Warning message:
'compare' is deprecated.
Use 'pcompare' instead.
See help("Deprecated") 
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>


Created & Maintained by Osamu Ogasawara (osamu.ogasawara@gmail.com) and IMSBIO