2.1 Starting with output from deepTools ‘computeMatrix’

If you do do not currently have deepTools installed, instructions for installation are in the deepTools manual, with the easiest method likely being through Bioconda.

This direct output from the ‘computeMatrix’ function can be imported into R as a profileplyr object using the import_deepToolsMat() function. This function takes as its only argument the path to the matrix from deepTools. The information contained within this matrix is stored in a profileplyr object.

library(profileplyr)

example <- system.file("extdata", 
                       "example_deepTools_MAT", 
                       package = "profileplyr") 
proplyrObject <- import_deepToolsMat(example) 

2.2 Starting with output from soGGi

The soGGi function plotRegion() allows for the quantification of signal over genomic intervals within R, and the information from this quantification is stored in a ChIPprofile object. profileplyr can take the ChIPprofile object as an input to generate a profileplyr object using the as_profileplyr() function. While the plotRegion function from soGGi does not allow for inputting multiple signal files or multiple BED files, profileplyr contains a function to do this, BamBigwig_to_ChIPprofile(). This function takes a character vector of paths to bigWig files or BAM files in the ‘signalFiles’ argument. Note that the files for each function call must be all BAMs or all bigWigs and not a combination of the two. The ‘testRanges’ argument is must be a character vector of paths to BED files. The corresponding names for each BED file that will appear on visualizations involving the resulting profileplyr object can be set with the ‘testRanges_names’ argument.

signalFiles <- c(system.file("extdata","Sorted_Hindbrain_day_12_1_filtered.bam",package = "profileplyr"),
                 system.file("extdata","Sorted_Liver_day_12_1_filtered.bam",package = "profileplyr"))

# BAMs must be indexed
require(Rsamtools)
for (i in seq_along(signalFiles)){
  indexBam(signalFiles[i])
}

testRanges <- system.file("extdata", 
                          "newranges_small.bed", 
                          package = "profileplyr")

chipProfile <- BamBigwig_to_chipProfile(signalFiles, 
                         testRanges, 
                         format = "bam",
                         paired=FALSE,
                         style="percentOfRegion",
                         nOfWindows=20,
                         distanceAround=40
                         )
chipProfile

## class: ChIPprofile 
## dim: 3 36 
## metadata(2): names AlignedReadsInBam
## assays(2): '' ''
## rownames(3): giID1 giID2 giID3
## rowData names(4): name score sgGroup giID
## colnames(36): Start-1 Start-2 ... End+7 End+8
## colData names(0):

Once we have the ChIPprofile object, the as_profileplyr() function will generate the profileplyr object that can be used in the many downstream functions described below. Similar to when you import from deepTools, the grouping of the ranges after using the as_profileplyr() function is stored in the ‘rowGroupsInUse’ section of params(proplyrObject). Further, if the object is derived from a soGGi ChIPprofile object, the inherited group column will be called ‘sgGroup’.

proplyrObject <- as_profileplyr(chipProfile)
params(proplyrObject)$rowGroupsInUse

## [1] "sgGroup"

3.1 Key components of the object

The matrices that represent the bigWig or BAM signal within each bin (bin size specified within deepTools, soGGi by default uses base pair resolution, so bin = 1) are contained in a list within the profileplyr object that can be accessed with the assays() function. The columns of each matrix are the bins, while the rows are the ranges from the BED files that were used as the input to deepTools or soGGi, and the dimensions of each matrix can be seen with dim().

library(SummarizedExperiment)

assays(proplyrObject)

## List of length 3
## names(3): hindbrain1_binned liver1_binned kidney1_binned

dim(proplyrObject)

## [1] 400  25

A key feature of the RangedSummarizedExperiment class is the ability to link the rows of the assay matrices to genomic ranges. Those ranges are stored in a GRanges object allowing for efficient integration with other Bioconductor packages. The GRanges object within the profileplyr object can be accessed with the rowRanges() function. Standard GRanges accessor functions such as start(), end(), or ranges() can be used on the entire profileplyr object. Furthermore, as these ranges are annotated by the various functions of profileplyr, information and classifications of these ranges will be stored in this GRanges object in the metadata columns, which can be accessed as a Dataframe with the mcols() function (this same information is also contained in the rowData section of the object, accessed with rowData()).

rowRanges(proplyrObject)[1:3]

## GRanges object with 3 ranges and 3 metadata columns:
##       seqnames              ranges strand |    names     score  dpGroup
##          <Rle>           <IRanges>  <Rle> | <factor> <numeric> <factor>
##   [1]    chr14   34880170-34881532      * |        .         0    genes
##   [2]    chr15   41155088-41156065      * |     ._r1         0    genes
##   [3]     chr1 135258238-135259293      * |     ._r2         0    genes
##   -------
##   seqinfo: 20 sequences from an unspecified genome; no seqlengths

The information associated with each sample is stored as a Dataframe within with the parameters of the profileplyr object, and is accessed with sampleData().

sampleData(proplyrObject)

## DataFrame with 3 rows and 19 columns
##                    upstream downstream      body  bin.size unscaled.5.prime
##                   <numeric>  <numeric> <numeric> <numeric>        <numeric>
## hindbrain1_binned         0          0      1000        40                0
## liver1_binned             0          0      1000        40                0
## kidney1_binned            0          0      1000        40                0
##                   unscaled.3.prime     sample_labels   verbose bin.avg.type
##                          <numeric>          <factor> <logical>     <factor>
## hindbrain1_binned                0 hindbrain1_binned     FALSE         mean
## liver1_binned                    0     liver1_binned     FALSE         mean
## kidney1_binned                   0    kidney1_binned     FALSE         mean
##                   missing.data.as.zero     scale skip.zeros nan.after.end
##                              <logical> <numeric>  <logical>     <logical>
## hindbrain1_binned                FALSE         1      FALSE         FALSE
## liver1_binned                    FALSE         1      FALSE         FALSE
## kidney1_binned                   FALSE         1      FALSE         FALSE
##                   proc.number sort.regions sort.using ref.point
##                     <numeric>     <factor>   <factor>    <list>
## hindbrain1_binned           1         keep       mean      NULL
## liver1_binned               1         keep       mean      NULL
## kidney1_binned              1         keep       mean      NULL
##                   min.threshold max.threshold
##                          <list>        <list>
## hindbrain1_binned          NULL          NULL
## liver1_binned              NULL          NULL
## kidney1_binned             NULL          NULL

The ‘rowGroupsInUse’ element for the params slot indicates which column of the range metadata (mcols() or rowRanges()) that will be used for grouping in the final output if that object were used for visualization. For a profileplyr object created from a deepTools computeMatrix output or from a soGGi ChIPprofile object, the inherited groups correspond to the BED files over which the signal was quantified, and these groups are contained in the ‘dpGroup’ column.

params(proplyrObject)$rowGroupsInUse

## [1] "dpGroup"

The ‘mcolToOrderBy’ slot of the profileplyr object parameters indicates which column of the range metadata will be used for ordering the ranges as they are exported to either deepTools or to EnrichedHeatmap (with the generateEnrichedHeatmap() function within profileplyr). This can be set using the orderBy() function, which requires a profileplyr object and a character string matching a column name of the range metadata as arguments. If groupBy is never used on an object, the default will be to order by the mean signal of each range (within each group). In addition, until groupBy() has been used on a profileplyr object, the ‘mcolToOrderBy’ slot of the parameters will be NULL and will not be seen with params(proplyrObject).

NOTE: If you are exporting to deepTools and want to order by something other than mean range signal (or other deepTools ordering options), then make sure to set the –sortRegions argument of deepTools ‘plotHeatmap’ to “no” so your custom ordering will be used. generateEnrichedHeatmap() will always use the column indicated by params(proplyrObject)$mcolToOrderBy, or the mean signal if that value is NULL.

params(proplyrObject)$mcolToOrderBy

## NULL

proplyrObject <- orderBy(proplyrObject, "score")
params(proplyrObject)$mcolToOrderBy

## [1] "score"

3.2 Subsetting the profileplyr object

The profileplyr object can be subset either by sample, or by rows and columns of the matrix for each sample. This is done using the ‘[ ]’ brackets, with the first position being assay matrix rows, the second position being assay matrix columns, and the third position being the entire matrix for each sample (in addition to the rest of the parameters and range data). For example if you wanted to get the first ten rows and columns of each sample matrix:

proplyrObject[1:10,1:10]

## class: profileplyr 
## dim: 10 10 
## metadata(0):
## assays(3): hindbrain1_binned liver1_binned kidney1_binned
## rownames: NULL
## rowData names(3): names score dpGroup
## colnames: NULL
## colData names(0):

The more useful subsetting functionality is likely the ability to subset by samples using the third position of the bracket. For example, if you only wanted the first two samples to export to a deepTools matrix or for further downstream analysis, you could simply just subset with the third position:

proplyrObject[ , , 1:2]

## class: profileplyr 
## dim: 400 25 
## metadata(0):
## assays(2): hindbrain1_binned liver1_binned
## rownames: NULL
## rowData names(3): names score dpGroup
## colnames: NULL
## colData names(0):

3.3 Changing sample names

It might be helpful to change the sample names to something that is shorter to make labeling of figures clearer. This can be done separately in most of the visualization packages that a profileplyr object can be used in, but you can change the names within in the sampleData it will be changed for all downstream analyses. After importing either a deepTools matrix or a soGGi object, the sample names are stored as the rownames of the sampleData Dataframe.

rownames(sampleData(proplyrObject))

## [1] "hindbrain1_binned" "liver1_binned"     "kidney1_binned"

The names can simply be changed by reassigning the rownames of this Dataframe.

rownames(sampleData(proplyrObject)) <- c("Hindbrain",
                                         "Liver",
                                         "Kidney")

Importantly, this also changes the names of the matricies stored in the assays() section of the profileplyr object, and the names used for any output method will also be changed (e.g. deepTools, EnrichedHeatmap, ggplot)

3.4 Connecting profileplyr functions using the pipe (%>%) operator

Code involving profileplyr objects can be made clearer using the pipe operator from the magrittr package. This is aided by the fact that the profileplyr object is both the output and input of most functions within the package. This is demonstrated throughout the vignette, including the more complex visualizations later in this vignette. The user can go all the way from importing a deepTools matrix or soGGi object, through annotation and grouping, and ending with export to either the generateEnrichedHeatmap() (shown below) or export_deepToolsMat() function

4.1 Export to a deepTools matrix

Any profileplyr object can be exported as a matrix formatted for the deepTools ‘plotHeatmap’ function using the export_deepToolsMat() function. This simply requires the name of the profileplyr object and the path to the desired destination of the gzipped matrix. This function will build the metadata for the deepTools matrix based on the parameters of the profileplyr object. Importantly, the range metadata column of profileplyr object (accessed with rowRanges(object)) that is specified in the ‘rowGroupsInUse’ section of the parameters for the whole object (found using params(object)$rowGroupsInUse) is automatically used to set the grouping of the exported deepTools matrix. The ‘rowGroupsInUse’ can be changed using the groupBy() function (see more details below).

output_path <- file.path(tempdir(),"ATAC_example.MAT.gz")
export_deepToolsMat(proplyrObject, con = output_path)

4.2 Directly generate a customized EnrichedHeatmap with annotated ranges

To generate a heatmap within R directly from the profileplyr object, the generateEnrichedHeatmap() function should be used. This function takes a profileplyr object and produces a heatmap using the EnrichedHeatmap package. It allows for easy export from the profileplyr object and simple inclusion of the metadata as heatmap annotations. So far the profileplyr object used in this vignette is relatively basic (no groups or additional metadata), but this function generates a multipanel heatmap with a variety of arguments that have been tailored to visualizing the profileplyr object. As we explore more functionality within profileplyr, some of these features will be demonstrated.

# heatmap not printed, see below for examples of rendered heatmaps from this function
heatmap <- generateEnrichedHeatmap(proplyrObject)
class(heatmap)

## [1] "HeatmapList"
## attr(,"package")
## [1] "ComplexHeatmap"

4.3 Convert to an EnrichedHeatmap matrix

profileplyr also facilitates direct use of the EnrichedHeatmap() function within the EnrichedHeatmap package, allowing the user further flexibility in the heatmaps that can be generated from a profileplyr object.EnrichedHeatmap uses a special type of matrix of the ‘normalizedMatrix’ class to generate heatmaps. The convertToEnrichedHeatmapMat() function within profileplyr takes a profileplyr object as the only required input, and then converts the matrices contained within the assays slot to a list of ‘normalizedMatrix’ class objects that can be used in the EnrichedHeatmap() function. See the EnrichedHeatmap vignette for detailed examples for how heatmaps can be concatenated together to visualize all of the data in one figure. Importantly, the metadata stored in the rowRanges slot of the profileplyr object can also be utilized to annotate these heatmaps. NOTE: While this function provides the most flexibility for the user to build a custom set of heatmaps, a more useful function for quick visualization of data stored in a profileplyr object is the generateEnrichedHeatmap() function, which builds these heatmaps with one command using the profileplyr object (see the previous section).

EH_mat <- convertToEnrichedHeatmapMat(proplyrObject)
EH_mat[[1]]

## Normalize Hindbrain to target:
##   Upstream 0 bp (0 window)
##   Downstream 0 bp (0 window)
##   Include target regions (25 windows)
##   400 target regions

# the matricies of this list can be inputs for EnrichedHeatmap()

5.1 Matrix output for heatmaps

If the ‘output’ argument is set to ‘matrix’, then only a matrix will be returned with a single column for each sample containing the bins summarized as indicated with the ‘fun’ argument. The row names of this matrix is a unique identifier for each range containing the chromosome, start, end, and group.

object_sumMat <- summarize(proplyrObject, 
                           fun = rowMeans, 
                           output = "matrix") 
object_sumMat[1:3, ]

##                                Hindbrain      Liver    Kidney
## chr14_34880170_34881532_genes  0.5134429 0.04936532 0.4721223
## chr15_41155088_41156065_genes  0.3062093 0.04938936 0.2161117
## chr1_135258238_135259293_genes 0.4648443 0.37157612 0.7157744

This matrix can be used directly in other heatmap generating packages, including heatmap or pheatmap.

5.2 Long output for ggplot

If the ‘output’ argument is set to ‘long’, then the output will be a long data frame that can be used for plotting with ggplot. The grouping column of the range metadata as specified by ‘params(proplyrObject)$rowGroupsInUse’ will automatically be included in the data frame. If the other range metadata columns should be included in the data frame, then the ‘keep_all_mcols’ argument should be set to TRUE. Additionally, columns specifying the range, as well as the sample and the summarized signal that correspond to that range are included by default.

proplyrObject_long <- summarize(proplyrObject, 
                                fun = rowMeans, 
                                output = "long") 
proplyrObject_long[1:3, ]

##   dpGroup          combined_ranges    Sample    Signal
## 1   genes  chr14_34880170_34881532 Hindbrain 0.5134429
## 2   genes  chr15_41155088_41156065 Hindbrain 0.3062093
## 3   genes chr1_135258238_135259293 Hindbrain 0.4648443

This data frame can then be used directly with ggplot for plotting.

It is often helpful to log transform the signal to more clearly see trends in the signal that is quantified.

library(ggplot2)
ggplot(proplyrObject_long, aes(x = Sample, y = log(Signal))) + 
       geom_boxplot()

5.3 profileplyr object output with summarized matrix

Lastly, if the ‘output’ argument is set to ‘object’, then a profileplyr object containing the summarized matrix will be returned. This will allow for further grouping or manipulation of the summarized ranges with other profileplyr functions, as opposed to using the binned ranges that are often used in the examples below.

proplyrObject_summ <- summarize(proplyrObject, 
                                fun = rowMeans, 
                                output = "object")

assays(proplyrObject_summ)[[1]][1:3, ]

##      Hindbrain      Liver    Kidney
## [1,] 0.5134429 0.04936532 0.4721223
## [2,] 0.3062093 0.04938936 0.2161117
## [3,] 0.4648443 0.37157612 0.7157744

rowRanges(proplyrObject_summ)[1:3]

## GRanges object with 3 ranges and 3 metadata columns:
##       seqnames              ranges strand |    names     score  dpGroup
##          <Rle>           <IRanges>  <Rle> | <factor> <numeric> <factor>
##   [1]    chr14   34880170-34881532      * |        .         0    genes
##   [2]    chr15   41155088-41156065      * |     ._r1         0    genes
##   [3]     chr1 135258238-135259293      * |     ._r2         0    genes
##   -------
##   seqinfo: 20 sequences from an unspecified genome; no seqlengths

6.1 K-means and hierarchical clustering of the genomic ranges

Clustering of the signal across the genomic intervals of interest can shed light on patterns that aren’t immediately apparent when looking at the data as a whole. The clusterRanges() function provides a framework to cluster ranges that are contained in a profileplyr object. It should be noted that clustering can also be performed outside of R within the deepTools ‘plotHeatmap’ function, however, profileplyr clustering allows for integration of this clustering with further grouping mechanisms and visualizations within R and profileplyr.

clusterRanges() takes a profileplyr object as its first argument as well as a function to summarize the signal in each range (similar to the summarize() function above). The pheatmap package is then used to cluster the ranges, and the type of clustering depends on whether the user inputs a value for ‘kmeans_k’ (for kmeans clustering) or ‘cutree_rows’ (for hierarchical clustering using hclust). An integer value entered for either of these arguments will specify the number of clusters used for each method. If both ‘kmeans_k’ or ‘cutree_rows’ are left as NULL (default), then a heatmap will be printed with hierarchical clustering, but no distinct clusters defined, and no profileplyr object will be returned. This might be helpful as an initial and quick look at the ranges or as a means to determine the number of clusters to try.

clusterRanges(proplyrObject, 
              fun = rowMeans)

# this code prints heatmap (does not return), but heatmap not shown here to save space

If an integer is entered for either ‘kmeans_k’ or ‘cutree_rows’, then a profileplyr object will be returned with a cluster assigned to each range, and a column is added to the range metadata with the cluster for each range. Further, the ‘rowGroupsInUse’ is changed to this column, meaning that if this object is exported as a deepTools matrix with the export_deepToolsMat() function, the heatmap generated by ‘plotHeatmap’ will be grouped by the clusters.

set.seed(0)
kmeans <- clusterRanges(proplyrObject, 
                        fun = rowMeans, 
                        kmeans_k = 4)

rowRanges(kmeans)[1:3]

## GRanges object with 3 ranges and 4 metadata columns:
##       seqnames              ranges strand |    names     score  dpGroup
##          <Rle>           <IRanges>  <Rle> | <factor> <numeric> <factor>
##   [1]    chr14   34880170-34881532      * |        .         0    genes
##   [2]    chr15   41155088-41156065      * |     ._r1         0    genes
##   [3]     chr1 135258238-135259293      * |     ._r2         0    genes
##         cluster
##       <ordered>
##   [1]         1
##   [2]         1
##   [3]         3
##   -------
##   seqinfo: 20 sequences from an unspecified genome; no seqlengths

params(kmeans)$rowGroupsInUse

## [1] "cluster"

To visually inspect the clusters, the heatmap can also be printed (though it will not be returned) if the ‘silent’ argument is set to FALSE. A profileplyr object will still be returned even if silent is set to FALSE.

Kmeans clustering:

set.seed(0)
kmeans <- clusterRanges(proplyrObject, 
                        fun = rowMeans, 
                        kmeans_k = 4, 
                        silent = FALSE)

Hierarchical Clustering:

hclust <- clusterRanges(proplyrObject, 
                        fun = rowMeans, 
                        cutree_rows = 4, 
                        silent = FALSE)

6.2 Output matrix with cluster information to deepTools

The profileplyr object that is the output of clusterRanges() can be exported as a deepTools matrix, which can then be used as an input for plotHeatmap.

output_path <- file.path(tempdir(),"kmeans_cluster_mat.gz")
export_deepToolsMat(kmeans, con = output_path)

This matrix can then be passed directly to ‘plotHeatmap’, either by using the terminal or by using the system() function in R.

Code from command line :

plotHeatmap -m kmeans_cluster_mat.gz -o kmeans_cluster_mat.jpg –startLabel start –endLabel end –xAxisLabel “distance (bp)”

6.3 Generate group-annotated heatmap in R directly with generateEnrichedHeatmap()

After clustering the profileplyr object can be passed directly as an argument into the generateEnrichedHeatmap() function, and by default the heatmap will be grouped and annotated by these clusters, which were automatically set as the ‘rowGroupsInUse’ in the clusterRanges() function. Assuming that the ‘include_group_annotation’ argument of this function is set to the default value of TRUE, whichever metadata column is set to the ‘rowGroupsInUse’ slot will be used for the grouping and annotation seen below to the left of the heatmap. If there are no range groups, then the user can set this argument to be FALSE, and those color annotations will be absent.

Further, the maximum value for the y-axis in the line plots at the top of each heatmap are also automatically set based on the highest mean range signal from all of the samples. This can be set manually with the ‘ylim’ argument, or if ylim = NULL, the maximum will be inferred (i.e. be different) for each individual heatmap.

generateEnrichedHeatmap(kmeans)

6.4 Visualize mean range signal for each cluster with ggplot

An alternative way to visualize the clusters is by summarizing the signal over each range and then plotting those means by cluster. The flexibility with grouping and annotating within profileplyr makes this relatively easy. This is also a good demonstration of using piped code, going all the way from import to clustering, then to summarizing by mean range signal, and finishing with ggplot visualization:

library(magrittr)

set.seed(0)
import_deepToolsMat(example) %>%
  clusterRanges(fun = rowMeans, 
                kmeans_k = 4, 
                silent = TRUE) %>% 
  summarize(fun = rowMeans, 
            output = "long") %>% 
  ggplot(aes(x = Sample, y = log(Signal))) + 
         geom_boxplot() + 
         facet_grid(~cluster) + 
         theme(axis.text.x = element_text(angle = 90, hjust = 1)) + 
         scale_x_discrete(labels= c("Hindbrain", "Liver", "Kidney"))

6.5 Gene annotation of ranges

Understanding the genes that are in the proximity to a particular set of ranges allows for potential understanding of functional consequences of signal patterns within these ranges. By annotating the ranges within a profileplyr object with genes, the user can connect signal of particular samples to specific genes and the functions attributed to those genes. profileplyr provides two methods of annotating the ranges within the object using two established Bioconductor packages, ChIPseeker and rGREAT.

anno <- annotateRanges(proplyrObject, 
                       TxDb = "mm10")

## >> preparing features information...      2019-05-12 11:37:21 PM 
## >> identifying nearest features...        2019-05-12 11:37:22 PM 
## >> calculating distance from peak to TSS...   2019-05-12 11:37:23 PM 
## >> assigning genomic annotation...        2019-05-12 11:37:23 PM 
## >> adding gene annotation...          2019-05-12 11:37:44 PM 
## >> assigning chromosome lengths           2019-05-12 11:37:44 PM 
## >> done...                    2019-05-12 11:37:44 PM

rowRanges(anno)[1:3]

## GRanges object with 3 ranges and 16 metadata columns:
##       seqnames              ranges strand |    names     score  dpGroup
##          <Rle>           <IRanges>  <Rle> | <factor> <numeric> <factor>
##   [1]    chr14   34880170-34881532      * |        .         0    genes
##   [2]    chr15   41155088-41156065      * |     ._r1         0    genes
##   [3]     chr1 135258238-135259293      * |     ._r2         0    genes
##                                                annotation   geneChr
##                                               <character> <integer>
##   [1] Intron (ENSMUST00000043349.6/14803, intron 2 of 15)        14
##   [2]                                   Distal Intergenic        15
##   [3]                                    Promoter (<=1kb)         1
##       geneStart   geneEnd geneLength geneStrand      geneId
##       <integer> <integer>  <integer>  <integer> <character>
##   [1]  34894609  34894706         98          1      723847
##   [2]  41173487  41173868        382          1   100418302
##   [3] 135253575 135258568       4994          2       13710
##                transcriptId distanceToTSS            ENSEMBL        SYMBOL
##                 <character>     <numeric>        <character>   <character>
##   [1]  ENSMUST00000083547.1        -13077 ENSMUSG00000065481        Mir346
##   [2]  ENSMUST00000060066.5        -17422               <NA> 4930555K19Rik
##   [3] ENSMUST00000003135.13             0 ENSMUSG00000003051          Elf3
##                                                                  GENENAME
##                                                               <character>
##   [1]                                                        microRNA 346
##   [2] leucine rich repeat and sterile alpha motif containing 1 pseudogene
##   [3]                                                   E74-like factor 3
##        annotation_short
##               <ordered>
##   [1]            Intron
##   [2] Distal Intergenic
##   [3]          Promoter
##   -------
##   seqinfo: 20 sequences from an unspecified genome; no seqlengths

anno_great <- annotateRanges_great(proplyrObject, 
                                   species = "mm10")
rowRanges(anno_great)[1:3]

## GRanges object with 3 ranges and 5 metadata columns:
##       seqnames              ranges strand |    names     score  dpGroup
##          <Rle>           <IRanges>  <Rle> | <factor> <numeric> <factor>
##   [1]     chr1 135258238-135259293      * |     ._r2         0    genes
##   [2]     chr4   97648929-97650011      * |     ._r4         0    genes
##   [3]     chr4   97648929-97650011      * |     ._r4         0    genes
##            SYMBOL distanceToTSS
##       <character>     <numeric>
##   [1]        Elf3          -294
##   [2]     Gm10192       -466304
##   [3]        Nfia       -128138
##   -------
##   seqinfo: 20 sequences from an unspecified genome; no seqlengths

7.1 Switch output grouping columns within existing range metadata

If the ‘group’ argument is a character string then it must match a name of one of the columns in mcols(proplyrObject). In this case the groupBy function will change the column that will be used for grouping a profileplyr object during visualization. While many of the profileplyr functions shown above will set the grouping column (specified in the ‘rowGroupsInUse’ section of params(proplyrObject)) to the appropriate column for that function, the user may want to use another column for the grouping of the deepTools matrix to be output. The user may have added an additional column to the range metadata that can be used, or a column that is generated in one of the above functions, but is not the default grouping column, might be useful. Here we add a column of random numbers to demonstate both adding a column and then swithcing to use this column in grouping.

mcols(proplyrObject)$newColumn <- rnorm(n = nrow(proplyrObject))

params(proplyrObject)$rowGroupsInUse

## [1] "dpGroup"

switchGroup <- groupBy(proplyrObject, 
                            group = "newColumn")
params(switchGroup)$rowGroupsInUse

## [1] "newColumn"

7.2 Group by user-supplied GRanges

If the ‘group’ argument is a GRanges or a GRangesList then the overlap between each GRanges in this list and the ranges of the profileplyr object will be determined. A column will be added to the range metadata (column name is ‘group_and_overlap’) that shows a combination of the inherited group and whether the range overlaps with any of the GRanges. The ‘rowGroupsInUse’ will be changed to this column by default. If the inherit_groups argument is set to FALSE, then only the overlaps with the GRanges will determine the grouping.

One consideration with this function is what should be done with ranges that overlap multiple GRanges. The default behavior of this function is to duplicate the range and place it in each overlap group for each GRanges. If separateDuplicated = TRUE, then separate groups will be generated for the ranges that overlap multiple GRanges so that every range is in the resulting profileplyr object exactly once. This option is discouraged with high numbers of input GRanges, as the combinations of overlapping groups can get quite large. Lastly, the non-overlapping regions can be included in the heatmap if the ‘include_nonoverlapping’ argument is set to TRUE.

Here we import BED files containing the top 5000 H3K27ac peaks from hindbrain and liver using the rtracklayer package. These peaks sets were adapted from BED files that were downloaded from ENCODE using the provided links. The GRanges are then combined into a GRangesList to be used in the ‘group’ argument for the groupBy() function.

#import pre-made GRangesList 
data("K27ac_GRlist_hind_liver_top5000")
K27ac_GRlist_names <- c("hindbrain_K27ac", 
                        "liver_K27ac")

# subset the profileplyr object to just hindbrain and liver, 
# and group by GRanges
K27ac_groupByGR <- proplyrObject[,,grepl("Hindbrain|Liver", 
                                         names(assays(proplyrObject)))] %>%
                   groupBy(group = K27ac_GRlist_hind_liver_top5000,
                           GRanges_names = K27ac_GRlist_names, 
                           include_nonoverlapping = FALSE,
                           inherit_groups = FALSE)

rowRanges(K27ac_groupByGR)[1:3]

## GRanges object with 3 ranges and 7 metadata columns:
##       seqnames            ranges strand |    names     score  dpGroup
##          <Rle>         <IRanges>  <Rle> | <factor> <numeric> <factor>
##   [1]    chr13 75707002-75708631      * |     ._r6         0    genes
##   [2]    chr10 98623831-98625379      * |    ._r17         0    genes
##   [3]     chr7 27490506-27491305      * |    ._r18         0    genes
##       overlap_matrix GR_overlap_names        name   score.1
##             <matrix>        <ordered> <character> <numeric>
##   [1]            0:1            liver        <NA>       233
##   [2]            0:1            liver        <NA>       247
##   [3]            0:1            liver        <NA>       281
##   -------
##   seqinfo: 20 sequences from an unspecified genome; no seqlengths

As with any profileplyr object, this can be exported to deepTools ‘plotHeatmap’ using the profileplyr function export_deepToolsMat()

output_path <- file.path(tempdir(),"K27ac_GRoverlap.gz")
export_deepToolsMat(K27ac_groupByGR, con = output_path)

This matrix can then be passed directly to ‘plotHeatmap’, either by using the terminal or by using the system() function in R.

Code from command line:

plotHeatmap -m K27ac_GRoverlap.gz -o K27ac_GRoverlap.jpg –startLabel start –endLabel end –xAxisLabel “distance (bp)”

The profileplyr object can also be passed directly to the generateEnrichedHeatmap() function to be visualized in R with the following command: generateEnrichedHeatmap(K27ac_groupByGR)

7.3 Group by user-supplied gene list

If the ‘group’ argument is a list (i.e. class(group) returns “list”), then it is assumed that this list contains lists of gene sets. These gene sets within this list can either be character vectors of gene symbols, or data frames with the gene symbols as rownames. Importantly, if all of the elements of the list are data frames AND they all have the same columns (as determined by having matching column names), then those columns will be added to the range metadata to potentially be used to annotate the ranges. This is especially useful for annotating the ranges with additional gene-level metadata (e.g. gene expression). However, if any of the elements of the list are character vectors, or if they are data frames with columns that are different, then the additional columns will not be included in the range metadata and only the overlap information will be stored.

The gene symbols will be used to determine overlap with the genes that are associated with each range of the profileplyr object. This type of grouping will typically follow one of the annotation functions: annotateRanges() or annotateRanges_great(). The names of the list elements (names(list)) will be used as the names of each set of genes in the range metadata and any exported deepTools matrix. A new column will be added to the range metadata that shows a combination of the inherited group and whether the range overlaps with any of the gene sets.

data("gene_list_character")

names(gene_list_character)

## [1] "upInHindbrain_VS_Liver"   "downInHindbrain_VS_Liver"

head(gene_list_character[[1]])

## [1] "Tuba1a" "Agap1"  "Nnat"   "Col2a1" "Rtn1"   "Map1b"

data("gene_list_dataframe")
names(gene_list_dataframe)

## [1] "upInHindbrain_VS_Liver"   "downInHindbrain_VS_Liver"

head(gene_list_dataframe[[1]])

##        hindbrain_liv_stat
## Tuba1a           75.80852
## Agap1            53.32196
## Nnat             51.68476
## Col2a1           49.94588
## Rtn1             48.51060
## Map1b            44.71733

signalFiles <- c(system.file("extdata","Sorted_Hindbrain_day_12_1_filtered.bam",
                             package = "profileplyr"),
                 system.file("extdata","Sorted_Liver_day_12_1_filtered.bam",
                             package = "profileplyr"))

# BAMs must be indexed
require(Rsamtools)
for (i in seq_along(signalFiles)){
  indexBam(signalFiles[i])
}

testRanges <- system.file("extdata", 
                          "newranges.bed", 
                          package = "profileplyr")

BamBigwig_to_chipProfile(signalFiles, 
                         testRanges, 
                         format = "bam",
                         paired=FALSE,
                         style="percentOfRegion",
                         nOfWindows=20,
                         distanceAround=40) %>% 
  as_profileplyr %>% 
  annotateRanges(TxDb = "mm10")  %>%
  groupBy(group = gene_list_dataframe) %>% 
  generateEnrichedHeatmap(extra_annotation_columns = c("hindbrain_liv_stat"),
                          sample_names = c("Hindbrain",
                                           "Liver"))