Visualization and annotation of read signal over genomic ranges with profileplyr
Introduction
Quick and straightforward visualization of read signal over genomic intervals is key for generating hypotheses from sequencing data sets (e.g. ChIP-seq, ATAC-seq, bisulfite/methyl-seq). Many tools both inside and outside of R and Bioconductor are available to explore these types of data, and they typically start with a bigWig or BAM file and end with some representation of the signal (e.g. heatmap). profileplyr leverages many Bioconductor tools to allow for both flexibility and additional functionality in workflows that end with visualization of the read signal. Signal over select genomic intervals can be quantified with Bioconductor packages in R (e.g. soGGi, EnrichedHeatmap, genomation or via the command line (e.g. deepTools, ngs.plot). profileplyr takes the signal quantification over ranges coming from either soGGi or deepTools, as well as the metadata associated with the ranges and the input bigWig or BAM files and stores this information in a RangedSummarizedExperiment object. profileplyr then uses this object to summarize and plot the signal within these ranges, and also to group and subset these ranges based on clustering, overlap with GRanges objects, overlap with gene lists, and genomic annotations. After manipulation with established Bioconductor tools, the profileplyr object can then be either directly visualized as a customized heatmap in R using EnrichedHeatmap or exported in a format that can be visualized with deepTools. While many users prefer to use deepTools certain steps of sequencing analysis (the handling of BAM files and utilizing computeMatrix/plotHeatmap), by staying outside of R one cannot take advantage of the many Bioconductor packages used to analyze genomic range data. The profileplyr object serves as an intermediate that can be easily manipulated and annotated using Bioconductor tools, and the integration of this object with deepTools, soGGi, and EnrichedHeatmap functions enhances the user’s ability to visualize genomic signal.
Import signal quantification from deepTools or soGGi
There are two main ways to go from a bigWig of BAM file to a profileplyr object, which is a version of the RangedSummarizedExperiment class within the SummarizedExperiment package. A profileplyr object can be generated from the command line with the output of the deepTools ‘computeMatrix’ function, or from within R using the output of the soGGi regionPlot() function.
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.
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’.
## [1] "sgGroup"The profileplyr object
The profileplyr object is a form of the RangedSummarizedExperiment class, which allows us to both store all of the relevant information that is imported from soGGi or deepTools, and have flexibility in how we manipulate this object in downstream analysis. We will go through many of these features by looking at the object we just generated from the import functions.
## class: profileplyr
## dim: 400 25
## metadata(0):
## assays(3): hindbrain1_binned liver1_binned kidney1_binned
## rownames: NULL
## rowData names(3): names score dpGroup
## colnames: NULL
## colData names(0):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().
## List of length 3
## names(3): hindbrain1_binned liver1_binned kidney1_binned## [1] 400 25A 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()).
## GRanges object with 3 ranges and 3 metadata columns:
## seqnames ranges strand | names score dpGroup
## <Rle> <IRanges> <Rle> | <character> <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 seqlengthsThe information associated with each sample is stored as a Dataframe within with the parameters of the profileplyr object, and is accessed with sampleData().
## DataFrame with 3 rows and 20 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> <character> <logical> <character>
## 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 min.threshold
## <numeric> <character> <character> <list> <list>
## hindbrain1_binned 1 keep mean
## liver1_binned 1 keep mean
## kidney1_binned 1 keep mean
## max.threshold generation.method
## <list> <character>
## hindbrain1_binned deepTools
## liver1_binned deepTools
## kidney1_binned deepToolsThe ‘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.
## [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.
## NULL## [1] "score"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:
## 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:
## 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):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.
## [1] "hindbrain1_binned" "liver1_binned" "kidney1_binned"The names can simply be changed by reassigning the rownames of this Dataframe.
Importantly, this also changes the names of the matrices 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)
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
Export/Conversion of profileplyr object for heatmap visualization of ranges
The profileplyr object allows for the flexibility to be exported to various types of visualization, both inside and outside of R. This flexibility coupled with the annotation and subsetting functions described later in this vignette provide great potential for quickly navigating through data sets. There are two main ways to visualize this data once the data is stored in a profileplyr object, first as a heatmap of signal over the individual bins within these genomic ranges (e.g. deepTools and EnrichedHeatmap), and second as summarized signal (e.g. mean of the entire range) in ggplot or clustered heatmaps (e.g. pheatmap).
For heatmap visualization of the signal over the genomic intervals this object can be converted to:
- A matrix that can be used as an input for the ‘plotHeatmap’ function of deepTools
- A heatmap containing all range grouping and annotation information (using a modified implementation of the EnrichedHeatmap constructor)
- A list of matrices that can be used in the EnrichedHeatmap() function within R
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).
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"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).
## Normalize Hindbrain to target:
## Upstream 0 bp (0 window)
## Downstream 0 bp (0 window)
## Include target regions (25 windows)
## 400 target regionsSummarize signal for ggplot or heatmap visualization
Oftentimes it will be important to use some kind of summary statistic for the entire range when trying to understand and visualize differences between samples. This is often going to be done using the summarize() function. The summarize function requires the name of a profileplyr object, the function used to summarize the ranges (e.g. rowMeans or rowMax), and the type of output. Here the basic options will be demonstrated with the example data which contains no groups, however, note how it is used later in the vignette when additional grouping options are discussed and summarizing the data in this way becomes especially useful.
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.
## 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.7157744This matrix can be used directly in other heatmap generating packages, including heatmap or pheatmap.
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.4648443This 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.
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.
## Hindbrain Liver Kidney
## [1,] 0.5134429 0.04936532 0.4721223
## [2,] 0.3062093 0.04938936 0.2161117
## [3,] 0.4648443 0.37157612 0.7157744## GRanges object with 3 ranges and 3 metadata columns:
## seqnames ranges strand | names score dpGroup
## <Rle> <IRanges> <Rle> | <character> <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 seqlengthsAnnotating of genomic ranges with clusters, genomic regions, and genes
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 spaceIf 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> | <character> <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## [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:
Hierarchical Clustering:
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)”
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.
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"))
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.
Annotation of ranges with genes and genomic regions using ChIPseeker
The annotateRanges() function passes the ranges of a profileplyr object to the annotatePeak() function of ChIPseeker, and then integrates the information from this output into the profileplyr object. The annotation of the ranges, including the annotation type (promoter, exon, etc) and closest gene, are then compiled in the range metadata. In addition to the profileplyr object, a TxDb object must be specified with the ‘TxDb’ argument. The TxDb argument can be one of three things:
- a character string that is either “hg19”, “hg38”, “mm9”, or “mm10”, which will automatically retrieve the relevant TxDb object
- a TxDb object
- a character string that is a path to a GTF or GFF file, in which case a TxDb object will be made using the makeTxDbFromGFF() function from the Genomic Features package.
Note that the default action of this function will not change the ‘rowGroupsInUse’ element of the params slot. However, if changeGroupToAnnotation = TRUE, a newly generated metadata column will be added that combines the inherited grouping of the ranges and the genomic annotation that was determined by ChIPseeker. If this object is exported as a deepTools matrix, the heatmap will group the ranges first by the inherited grouping, followed by the annotation groups (e.g. promoter, intron, etc). This can be modified by setting the ‘heatmap_grouping’ argument to ‘annotation’, which will force the ranges to be grouped first by annotation, followed by the inherited grouping.
## >> preparing features information... 2024-10-31 03:45:29 AM
## >> identifying nearest features... 2024-10-31 03:45:29 AM
## >> calculating distance from peak to TSS... 2024-10-31 03:45:30 AM
## >> assigning genomic annotation... 2024-10-31 03:45:30 AM
## >> adding gene annotation... 2024-10-31 03:45:43 AM## >> assigning chromosome lengths 2024-10-31 03:45:43 AM
## >> done... 2024-10-31 03:45:43 AM## GRanges object with 3 ranges and 16 metadata columns:
## seqnames ranges strand | names score dpGroup
## <Rle> <IRanges> <Rle> | <character> <numeric> <factor>
## [1] chr14 34880170-34881532 * | . 0 genes
## [2] chr15 41155088-41156065 * | ._r1 0 genes
## [3] chr1 135258238-135259293 * | ._r2 0 genes
## annotation geneChr geneStart geneEnd geneLength
## <character> <integer> <integer> <integer> <integer>
## [1] Intron (ENSMUST00000.. 14 34894609 34894706 98
## [2] Distal Intergenic 15 41173487 41173868 382
## [3] Promoter (<=1kb) 1 135253575 135258568 4994
## geneStrand geneId transcriptId distanceToTSS
## <integer> <character> <character> <numeric>
## [1] 1 723847 ENSMUST00000083547.1 -13077
## [2] 1 100418302 ENSMUST00000060066.5 -17422
## [3] 2 13710 ENSMUST00000003135.13 0
## ENSEMBL SYMBOL GENENAME annotation_short
## <character> <character> <character> <ordered>
## [1] ENSMUSG00000065481 Mir346 microRNA 346 Intron
## [2] <NA> 4930555K19Rik leucine rich repeat .. Distal Intergenic
## [3] ENSMUSG00000003051 Elf3 E74-like factor 3 Promoter
## -------
## seqinfo: 20 sequences from an unspecified genome; no seqlengthsLastly, if only a subset of annotation types are desired (e.g. just those ranges within promoters) then the ‘annotation_subset’ argument can be used and the profileplyr object will be subset accordingly.
Adding metadata information beyond the grouping column to the EnrichedHeatmap
The EnrichedHeatmap package brings a lot of heatmap annotation options, and the profileplyr generateEnrichedHeatmap() function takes advantage of this by giving the user the ability to add both categorical and numeric annotations on the right side of the matrix heatmaps. For example, the ChIPseeker annotations added to the range metadata in the previous example can be visualized. Adding metadata columns is done using the ‘extra_annotation_columns’ argument of generateEnrichedHeatmap(). This argument takes a character vector with strings that match column names in the range metadata, with no limit on the number of columns that can be visualized. For example, to generate a heatmap directly including annotation, and using the profileplyr object from the previous chunk, it would be a simple one line command: generateEnrichedHeatmap(anno, extra_annotation_columns = “annotation_short”)
In addition to being exported to a deepTools matrix or EnrichedHeatmap, these annotations from annotateRanges() can be further utilized with the groupBy() function, discussed below.
Annotation of ranges with genes using GREAT
Another common method of annotating genomic intervals with nearby genes is by using the GREAT method, which defines regulatory regions for each gene and then determines whether a certain genomic interval overlaps with that genic regulatory region. The profileplyr function annotateRanges_great() takes a profileplyr object and the genome that is to be used for the annotation, and will add a column to the range metadata to indicate whether that range overlaps with the regulatory region of a gene. If a region is in regulatory regions of multiple genes, then it will have multiple rows within the profileplyr object. This function will not change the grouping of the profileplyr object.
## GRanges object with 3 ranges and 5 metadata columns:
## seqnames ranges strand | names score dpGroup
## <Rle> <IRanges> <Rle> | <character> <numeric> <factor>
## [1] chr14 34880170-34881532 * | . 0 genes
## [2] chr15 41155088-41156065 * | ._r1 0 genes
## [3] chr1 135258238-135259293 * | ._r2 0 genes
## SYMBOL distanceToTSS
## <list> <list>
## [1] Grid1,4930596D02Rik 61026,931127
## [2] Oxr1,Zfpm2 -291905,500535
## [3] Elf3 -198
## -------
## seqinfo: 20 sequences from an unspecified genome; no seqlengthsGrouping ranges by range metadata, gene list, or additional GRanges
The groupBy() function within profileplyr allows for a range of grouping mechanisms. The groupBy() function always requires a profileplyr object and the ‘group’ argument, which will determine how the ranges are to be grouped. There are three options for grouping:
- Switch output grouping columns within existing range metadata
- Group by user-supplied GRanges
- Group by user-supplied gene list
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 demonstrate both adding a column and then switching to use this column in grouping.
mcols(proplyrObject)$newColumn <- rnorm(n = nrow(proplyrObject))
params(proplyrObject)$rowGroupsInUse## [1] "dpGroup"## [1] "newColumn"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)## GRanges object with 3 ranges and 7 metadata columns:
## seqnames ranges strand | names score dpGroup
## <Rle> <IRanges> <Rle> | <character> <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 seqlengthsAs 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)
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.
Gene set list contains a character data
The first option is to read in a list of gene sets that are character vectors. This list contains the top 1000 differentially expressed genes that are higher in the liver compared to hindbrain, and the top 1000 that are higher in the hindbrain compared to the liver.
## [1] "upInHindbrain_VS_Liver" "downInHindbrain_VS_Liver"## [1] "Tuba1a" "Agap1" "Nnat" "Col2a1" "Rtn1" "Map1b"Gene set list with data frames for heatmap annotation (using pipe - %>%)
Using the same list of genes as above, we can also use them as part of a list of data frames (the gene symbols must be the rownames) with a column containing gene expression changes for hindbrain vs liver (they are values from the ‘stat’ column of a DESeq2 results table). This allows us to add a heatmap to our visualization showing these expression changes.
## [1] "upInHindbrain_VS_Liver" "downInHindbrain_VS_Liver"## hindbrain_liv_stat
## Tuba1a 75.80852
## Agap1 53.32196
## Nnat 51.68476
## Col2a1 49.94588
## Rtn1 48.51060
## Map1b 44.71733signalFiles <- 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")Using the list of data frames as our gene list, we produce a heatmap with an additional column showing gene expression changes in hindbrain vs liver for each range. This is achieved by annotating these ranges with genes using the ‘annotateRanges’ function which utilizes the ChIPseeker package. This can also be done with the ‘annotateRanges_great’ function, which uses the rGREAT package for annotation. The resulting profileplyr object has range metadata columns that have gene annotation (‘SYMBOL’ column) and a column with the gene set that each range overlapped combined with the inherited grouping (‘group_and_overlap’ column).
Also note that the generateEnrichedHeatmap() function gives the ability to manually change the sample names (previously we had shortened them by changing the rownames of sampleData(proplyrObject))
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) %>% #subseting by our gene list
generateEnrichedHeatmap(extra_annotation_columns = c("hindbrain_liv_stat"),
sample_names = c("Hindbrain",
"Liver"))
Combining multiple profileplyr functions for heatmap annotation
As a demonstration, we can combine many of the profileplyr functions and visualize more than one metadata column. From this figure we can see a number of interesting things from the data. Many of the regions with the highest ATAC-seq signal are in promoters for all clusters (the ‘annotation_short’ column). Cluster 2 appears to have many ranges with specific signal in the liver sample which is supported by the main heatmap, the high density of top liver K27ac peaks (red signal in the ‘GRoverlap_names’ column), and the high density of genes that have reduced gene expression in the liver vs the hindbrain (red signal in the ‘stat’ column). On the other hand Cluster 1 has specific hindbrain peaks, high overlap with hindbrain K27ac peaks, and higher expression of genes that have higher expression in the hindbrain. By extracting all of this information from the profileplyr object and displaying in a heatmap with one command, a lot can be gleaned from these data.
The generateEnrichedHeatmap() function also allows for a great deal of customization of the heatmap that is built around the profileplyr object. Notice that in the heatmap below we modify colors. The color of the heatmaps from the assay matrices is changed to a two color spectrum, and the colors of the annotation heatmaps were customized as well. The user can also easily label ranges that have been annotated with specific genes with the ‘genes_to_label’ argument. This argument requires a character vector of genes that are matched against the column of the range metadata that contains the gene annotations from either annotateRanges() or annotateRanges_great().
# NOTE: could also use BamBigwig_to_chipProfile() and as_profileplyr()
# within R to start from bigwig/BAM and BED file.
set.seed(0)
annotated_object <- import_deepToolsMat(example) %>%
annotateRanges(TxDb = "mm10") %>% # annotate with ChIPseeker
groupBy(group = gene_list_dataframe,
include_nonoverlapping = TRUE,
inherit_groups = FALSE) %>% # group by gene overlap
groupBy(group = K27ac_GRlist_hind_liver_top5000,
GRanges_names = K27ac_GRlist_names,
include_nonoverlapping = TRUE,
inherit_groups = FALSE) %>% # group by GRanges
clusterRanges(kmeans_k = 3) %>%
generateEnrichedHeatmap(extra_annotation_columns=c("annotation_short",
"GR_overlap_names",
"hindbrain_liv_stat"),
sample_names = c("Hindbrain",
"Liver",
"Kidney"),
extra_anno_width = c(5,5,8),
matrices_color = c("white", "red"),
extra_anno_color = list(NULL,
c("black",
"green",
"red",
"white"),
c("#e7e1ef",
"#c994c7",
"#dd1c77")),
genes_to_label = c("Myh9", "Sp4", "Rcan2"),
gene_label_font_size = 10)
The color arguments for both the matrices and the extra annotation
columns can either be a character vector if all heatmaps should have the
same color scheme, or a list if each individual heatmap requires
different colors. If a list is used, the list must be the same length as
the number of heatmaps, and each element of the list maps to the
corresponding matrix in the profileplyr object (for ‘matrices_color’) or
to the corresponding element in the ‘extra_annotation_columns’ argument
(for ‘extra_anno_color’). If the color of a particular heatmap is to be
left unchanged, then that element of the list should be NULL.
Acknowledgements
Thank you to Ji-Dung Luo for testing/vignette review/critical feedback, Elitsa Stoyanova for critical feedback/vignette review and David Allis and Ziwei Liang for their support.
Session info
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] org.Mm.eg.db_3.20.0 AnnotationDbi_1.69.0
## [3] magrittr_2.0.3 ggplot2_3.5.1
## [5] Rsamtools_2.22.0 Biostrings_2.75.0
## [7] XVector_0.46.0 profileplyr_1.23.0
## [9] SummarizedExperiment_1.36.0 Biobase_2.67.0
## [11] GenomicRanges_1.59.0 GenomeInfoDb_1.43.0
## [13] IRanges_2.41.0 S4Vectors_0.44.0
## [15] MatrixGenerics_1.19.0 matrixStats_1.4.1
## [17] BiocGenerics_0.53.0 BiocStyle_2.35.0
##
## loaded via a namespace (and not attached):
## [1] splines_4.4.1
## [2] later_1.3.2
## [3] BiocIO_1.17.0
## [4] bitops_1.0-9
## [5] ggplotify_0.1.2
## [6] tibble_3.2.1
## [7] R.oo_1.26.0
## [8] preprocessCore_1.68.0
## [9] XML_3.99-0.17
## [10] lifecycle_1.0.4
## [11] pwalign_1.2.0
## [12] doParallel_1.0.17
## [13] lattice_0.22-6
## [14] sass_0.4.9
## [15] rmarkdown_2.28
## [16] jquerylib_0.1.4
## [17] yaml_2.3.10
## [18] plotrix_3.8-4
## [19] httpuv_1.6.15
## [20] ggtangle_0.0.4
## [21] cowplot_1.1.3
## [22] DBI_1.2.3
## [23] buildtools_1.0.0
## [24] RColorBrewer_1.1-3
## [25] ShortRead_1.64.0
## [26] abind_1.4-8
## [27] zlibbioc_1.52.0
## [28] purrr_1.0.2
## [29] R.utils_2.12.3
## [30] RCurl_1.98-1.16
## [31] yulab.utils_0.1.7
## [32] circlize_0.4.16
## [33] GenomeInfoDbData_1.2.13
## [34] enrichplot_1.27.0
## [35] ggrepel_0.9.6
## [36] tidytree_0.4.6
## [37] maketools_1.3.1
## [38] pheatmap_1.0.12
## [39] ChIPseeker_1.43.0
## [40] codetools_0.2-20
## [41] DelayedArray_0.33.1
## [42] DOSE_4.1.0
## [43] DT_0.33
## [44] tidyselect_1.2.1
## [45] shape_1.4.6.1
## [46] aplot_0.2.3
## [47] UCSC.utils_1.2.0
## [48] farver_2.1.2
## [49] GenomicAlignments_1.43.0
## [50] jsonlite_1.8.9
## [51] GetoptLong_1.0.5
## [52] iterators_1.0.14
## [53] foreach_1.5.2
## [54] tools_4.4.1
## [55] progress_1.2.3
## [56] treeio_1.30.0
## [57] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
## [58] Rcpp_1.0.13
## [59] glue_1.8.0
## [60] SparseArray_1.6.0
## [61] xfun_0.48
## [62] qvalue_2.38.0
## [63] dplyr_1.1.4
## [64] withr_3.0.2
## [65] BiocManager_1.30.25
## [66] fastmap_1.2.0
## [67] latticeExtra_0.6-30
## [68] boot_1.3-31
## [69] fansi_1.0.6
## [70] caTools_1.18.3
## [71] digest_0.6.37
## [72] R6_2.5.1
## [73] mime_0.12
## [74] gridGraphics_0.5-1
## [75] colorspace_2.1-1
## [76] Cairo_1.6-2
## [77] GO.db_3.20.0
## [78] gtools_3.9.5
## [79] jpeg_0.1-10
## [80] RSQLite_2.3.7
## [81] R.methodsS3_1.8.2
## [82] utf8_1.2.4
## [83] tidyr_1.3.1
## [84] generics_0.1.3
## [85] data.table_1.16.2
## [86] rtracklayer_1.66.0
## [87] prettyunits_1.2.0
## [88] httr_1.4.7
## [89] htmlwidgets_1.6.4
## [90] S4Arrays_1.6.0
## [91] pkgconfig_2.0.3
## [92] gtable_0.3.6
## [93] blob_1.2.4
## [94] hwriter_1.3.2.1
## [95] ComplexHeatmap_2.23.0
## [96] sys_3.4.3
## [97] chipseq_1.57.0
## [98] htmltools_0.5.8.1
## [99] fgsea_1.33.0
## [100] clue_0.3-65
## [101] scales_1.3.0
## [102] TxDb.Hsapiens.UCSC.hg38.knownGene_3.20.0
## [103] png_0.1-8
## [104] ggfun_0.1.7
## [105] knitr_1.48
## [106] reshape2_1.4.4
## [107] rjson_0.2.23
## [108] nlme_3.1-166
## [109] curl_5.2.3
## [110] org.Hs.eg.db_3.20.0
## [111] cachem_1.1.0
## [112] GlobalOptions_0.1.2
## [113] stringr_1.5.1
## [114] KernSmooth_2.23-24
## [115] parallel_4.4.1
## [116] restfulr_0.0.15
## [117] pillar_1.9.0
## [118] grid_4.4.1
## [119] vctrs_0.6.5
## [120] gplots_3.2.0
## [121] promises_1.3.0
## [122] cluster_2.1.6
## [123] xtable_1.8-4
## [124] evaluate_1.0.1
## [125] rGREAT_2.8.0
## [126] EnrichedHeatmap_1.37.0
## [127] GenomicFeatures_1.59.0
## [128] locfit_1.5-9.10
## [129] cli_3.6.3
## [130] compiler_4.4.1
## [131] rlang_1.1.4
## [132] crayon_1.5.3
## [133] labeling_0.4.3
## [134] interp_1.1-6
## [135] plyr_1.8.9
## [136] fs_1.6.5
## [137] stringi_1.8.4
## [138] TxDb.Mmusculus.UCSC.mm10.knownGene_3.10.0
## [139] deldir_2.0-4
## [140] BiocParallel_1.41.0
## [141] munsell_0.5.1
## [142] lazyeval_0.2.2
## [143] tiff_0.1-12
## [144] GOSemSim_2.33.0
## [145] Matrix_1.7-1
## [146] hms_1.1.3
## [147] patchwork_1.3.0
## [148] bit64_4.5.2
## [149] KEGGREST_1.47.0
## [150] shiny_1.9.1
## [151] highr_0.11
## [152] TxDb.Mmusculus.UCSC.mm9.knownGene_3.2.2
## [153] igraph_2.1.1
## [154] memoise_2.0.1
## [155] bslib_0.8.0
## [156] ggtree_3.15.0
## [157] fastmatch_1.1-4
## [158] soGGi_1.38.0
## [159] bit_4.5.0
## [160] ape_5.8- Introduction
- Import signal quantification from deepTools or soGGi
- The profileplyr object
- Export/Conversion of profileplyr object for heatmap visualization of ranges
- Summarize signal for ggplot or heatmap visualization
- Annotating of genomic ranges with clusters, genomic regions, and genes
- Grouping ranges by range metadata, gene list, or additional GRanges
- Combining multiple profileplyr functions for heatmap annotation
- Acknowledgements
- Session info