Argument Description
object A profileplyr object
con Connection to write/read deeptools data to/from.
decreasing If object@params$mcolToOrderBy has been set and not NULL, then the ranges will be ordered by the column indicated in this slot of the metadata. By default, the order will be increasing for the factor or numeric value. For decreasing order, choose decreasing = TRUE.
overwrite Logical specifying whether to overwrite output if it exists.

Example

Here we load a pre-compiled profileplyr object. This profileplyr object is then grouped by gene lists, and exported out to deeptools for visualization.

Group by user-supplied gene list

A profileplyr object can also be subset or annotated based on overlap of its gene annotations (from annotateRanges() or annotateRanges_great()) with a list of genes from the user. This argument must be a list(), and the elements of list can either be a character vector of gene symbols or a list of data frames with the gene symbols as the rownames. Additional columns in these data frames will be inculded as columns in the range metadata and can also be used for additional range annotation.

Read in gene lists

First we will read in two genes lists, one for genes going up in the mesoderm compared to the ESCs, and one for the genes that go down (or up in ESC). This analysis was done with DESeq2 using the following counts tables downloaded from ENCODE: * HUES64 ESC: Replicate 1 and Replicate 2 * Mesoderm: Replicate 1 and Replicate 2

library(dplyr)

UP_in_meso_genes <- read.table("UP_in_meso_genes.txt",
                               header = FALSE,
                               stringsAsFactors = FALSE) %>%
  .[,1]

head(UP_in_meso_genes)
## [1] "HAPLN1" "COL3A1" "ANXA1"  "ITGA1"  "LOX"    "COL1A1"
DOWN_in_meso_genes <- read.table("DOWN_in_meso_genes.txt",
                               header = FALSE,
                               stringsAsFactors = FALSE) %>%
  .[,1]

head(UP_in_meso_genes)
## [1] "HAPLN1" "COL3A1" "ANXA1"  "ITGA1"  "LOX"    "COL1A1"

groupBy the gene lists

We combine these gene sets into a single list, and again use the groupBy() function to group the ranges. This will create a column in the range metadata shoing these overlaps, and the’rowgroupsInUse’ will by default be switched to this column. Remember that you can always switch this grouping to any column you’d like with the same groupBy() function, as shown above.

# IMPORTANT: The profileplyr object must have been annotated with either annotateRanges() or annotateRanges_great() before subsetting by a gene list

proplyr_object_subset <- annotateRanges_great(proplyr_object_subset,
                                              species = "hg19")

geneList_for_grouping <- list(meso_UP= UP_in_meso_genes,
                              meso_DOWN = DOWN_in_meso_genes)

proplyr_object_subset <- groupBy(proplyr_object_subset,
                                 group = geneList_for_grouping)
## A column has been added to the range metadata with the column name 'GL_overlap_names' that specifies the GRanges each range overlaps with, but the inherited groups are not included.
mcols(proplyr_object_subset)
## DataFrame with 843 rows and 9 columns
##                 name     score                sgGroup     giID    names
##          <character> <numeric>               <factor> <factor> <factor>
## giID808           NA         0 K27ac_top10_HUES64.bed  giID808       NA
## giID785           NA         0 K27ac_top10_HUES64.bed  giID785       NA
## giID736           NA         0 K27ac_top10_HUES64.bed  giID736       NA
## giID714           NA         0 K27ac_top10_HUES64.bed  giID714       NA
## giID714           NA         0 K27ac_top10_HUES64.bed  giID714       NA
## ...              ...       ...                    ...      ...      ...
## giID983           NA         0 K27ac_top10_HUES64.bed  giID983       NA
## giID983           NA         0 K27ac_top10_HUES64.bed  giID983       NA
## giID999           NA         0 K27ac_top10_HUES64.bed  giID999       NA
## giID999           NA         0 K27ac_top10_HUES64.bed  giID999       NA
## giID1000          NA         0 K27ac_top10_HUES64.bed giID1000       NA
##               SYMBOL distanceToTSS overlap_matrix GL_overlap_names
##          <character>     <numeric>       <matrix>        <ordered>
## giID808        PUSL1          6780            0:1        meso_DOWN
## giID785         RER1       -112760            0:1        meso_DOWN
## giID736         RERE        146478            0:1        meso_DOWN
## giID714     SLC25A33       -167481            0:1        meso_DOWN
## giID714        SPSB1         79121            1:0          meso_UP
## ...              ...           ...            ...              ...
## giID983        HCFC1          -690            0:1        meso_DOWN
## giID983      TMEM187           170            1:0          meso_UP
## giID999       NLGN4Y       -721757            0:1        meso_DOWN
## giID999       TMSB4Y         98422            0:1        meso_DOWN
## giID1000       KDM5D           105            0:1        meso_DOWN

Export to deepTools

The profileplyr object can always be exported to deepTools using the ‘export_deepToolsMat’ function, however, remember that while the grouping will be preserved in deepTools, the extra annotation columns will not.

export_deepToolsMat(proplyr_object_subset,
                    con = "groupBy_geneList_deeptools",
                    overwrite = TRUE)
## [1] "groupBy_geneList_deeptools"
##system("plotHeatmap -m groupBy_geneList_deeptools.gz -o groupBy_geneList_deeptools.pdf --yMax 50 5 100 50 5 150 --zMax 100 20 180 100 20 180")

igv

Visualize same figure within R and profileplyr

generateEnrichedHeatmap(proplyr_object_subset,
                        include_group_annotation = TRUE,
                        color_by_sample_group = "chip",
                        ylim = "chip",
                        show_heatmap_legend = rep(FALSE, 6)
                        )