1 Introduction
2 Installation
3 seqArchRplus for downstream analysis of promoter sequences
4 Conclusion
5 Session Info
References

1 Introduction

Analysis of the promoterome of any organism entails analyzing data from an experiment like CAGE (Cap analysis of gene expression) (Kodzius et al. 2006) that provides information on genome-wide transcription start sites at single nucleotide resolution. These promoters can then be further studied to identify different classes (Carninci et al. 2006) among them based on different attributes, for instance, their shape (broad or focused/sharp promoters), gene function (tissue specific vs. housekeeping) etc. These different promoter classes harbour a variety of promoter architectures orchestrated by different proteins together with sequence elements at near-fixed positions in the promoter that determine the position of the transcription start site. The different promoter architectures are known to be used differentially by genes either in different tissues or at different developmental timepoints. For instance, (Haberle et al. 2014) have shown a dynamic interchange of promoter architecture within the same genomic region between maternal and zygotic stages of zebrafish development. Identifying and studying these different promoter architectures further is thus vital.

While the R/Bioconductor package seqArchR enables de novo identification of clusters of promoter sequence architectures (Nikumbh and Lenhard 2023), this package, seqArchRplus, enables performing various steps in their downstream bioinformatic analyses and produce publication-ready plots (building on various other Bioconductor packages). The many steps in the downstream analyses of promoter sequence architectures enabled by seqArchRplus are:

Curate the final set of clusters from seqArchR. See the seqArchR vignette , or the bioRxiv preprint to understand in detail why this may be required
Order the sequence architectures by the interquantile widths (IQWs) of the tag clusters (aka promoter shape). See CAGEr vignette for more information on tag clusters and their IQWs
Visualize distributions of IQW, TPM (Tags per million) and conservation scores (or other when available) per cluster
Visualize the percentages of annotations for genome-wide CAGE-derived transcription tag clusters for each architecture-based clusters
Ease of comparison across samples/stages: Visualize the above plots as (publication ready) combined panels
Many per cluster visualizations including:
- sequence logos of cluster architectures
- strand-separated sequence logos of architectures
- distributions of promoters on different chromosomes and strands
- GO terms enriched for each cluster/architecture
Produce BED track files of seqArchR clusters for visualization in a genome browser or IGV

Examples of most of these capabilities are provided in this vignette.

We hope that using seqArchRplus (together with seqArchR) will be useful in swiftly, but comprehensively, analyzing promoters identified using CAGE for any organism, and enable insights and hypotheses generation with ease. So far, we have tested seqArchR and seqArchRplus on promoters in drosophila from (Schor et al. 2017) and modENCODE (Chen et al. 2014), zebrafish (Nepal et al. 2013), mice (Fantom (The FANTOM Consortium 2014)), humans (ENCODE), and also in plants (barley and maize).

We already have a plan for a feature in the future version of seqArchRplus: generate HTML reports that help you navigate this wealth of information with relative ease.

2 Installation

The latest stable version of seqArchjRplus can be installed from Bioconductor as shown below.

## When available on Bioconductor
if (!require("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}
BiocManager::install("seqArchRplus")

In case of any errors or help required, please consider looking up: https://github.com/snikumbh/seqArchRplus and file a new issue.

3 seqArchRplus for downstream analysis of promoter sequences

# Load seqArchRplus
library(seqArchRplus)
library(seqArchR)
library(Biostrings)
library(TxDb.Dmelanogaster.UCSC.dm6.ensGene)
library(org.Dm.eg.db)
library(ChIPseeker)

# Set seed for reproducibility
set.seed(1234)

3.1 Setting up

To begin using seqArchRplus, you require a set of promoters sequences (as a DNAStringSet object) and their clustering information (as a simple list storing the sequence IDs belonging to each cluster). With this, you can already visualize the cluster-wise sequence logos, distribution of chromosome and strand locations, GO term-enrichments, motif occurrence heatmaps. When you have additional information such as the IQWs (shape information), seqArchRplus can use these to order the clusters in the visualizations. This visualization can be supplemented with the distribution of per-cluster TPM values and/or conservation scores etc. when provided.

If your workflow involves CAGEr (for pre-processing raw CAGE data) and seqArchR (for clusters), they can be seamlessly used to utilize the full scope of seqArchRplus. In this case, it will be helpful if you have

the CAGEr object or information on the tag clusters from https://bioconductor.org/packages/release/bioc/html/CAGEr.html, specifically, the width of the tag clusters, total TPM value of the cluster and that of the dominant CTSS in the cluster
the seqArchR result object.

3.1.1 Input data

In this vignette, we use the example data provided with the package. The raw CAGE data for different developmental timepoints in Drosophila melanogaster (Schor et al. 2017) were pre-processed with CAGEr to obtain promoter sequences. These were then processed with seqArchR to cluster the promoter DNA sequences. Only a subset of the complete data is provided with the package and used here to enable demonstration of package utility. In particular, the data included is for the timepoint 10-12h after egg laying.

Specifically, the following files are provided:

promoter sequences stored as a gzipped FASTA file: two kinds, with 45 and 200 bp flanks around the TSS (filenames: example_promoters45.fa.gz and example_promoters200.fa.gz)
clustering information (filename: example_clust_info.rds)
GRanges object obtained from CAGEr which holds information on the tag clusters (promoters) with additional information such as the position of dominant TSS, TPM value for the dominant TSS, position of quantile boundaries (filename: example_tc_gr.rds)
BED file with same information as the GRanges object (filename: example_info_df.rds)
seqArchR result object (filename: seqArchR_result.rds)

## 1. Raw DNA sequences
raw_seqs <- Biostrings::readDNAStringSet(
    filepath = system.file("extdata",
        "example_promoters45.fa.gz",
        package = "seqArchRplus",
        mustWork = TRUE
    )
)

## 2. Clustering information (arbitrary order/unordered)
unord_clusts <- readRDS(system.file("extdata", "example_clust_info.rds",
                            package = "seqArchRplus", mustWork = TRUE
                        ))

## 3. GRanges object
tc_gr <- readRDS(system.file("extdata", "example_tc_gr.rds",
                    package = "seqArchRplus", mustWork = TRUE
                ))

## 4. BED file
bed_info_fname <- system.file("extdata", "example_info_df.bed.gz",
                                package = "seqArchRplus", mustWork = TRUE
                            )
info_df <- read.delim(file = bed_info_fname, sep = "\t", header = TRUE)

## 5. seqArchR result
seqArchR_result <- readRDS(system.file("extdata", "seqArchR_result.rds",
      package = "seqArchRplus", mustWork = TRUE))

## **NOTE** Only for the example seqArchR result object provided with this
## package:
## 
## Any seqArchR result object has the raw DNA sequences as a part of it.
## See details here: https://snikumbh.github.io/seqArchR/reference/seqArchR.html
## 
## But, in order to reduce the size of the example seqArchR result object 
## provided with this (seqArchRplus) package, the `rawSeqs` element of the 
## result is missing. This is reassigned below for the purposes of this 
## vignette.
## 
## The seqArchR result was obtained by processing the DNA sequences in 
## `example_promoters45.fa.gz`. Thus we reassign them to 
## `seqArchR_result$rawSeqs` here
##  
if(!("rawSeqs" %in% names(seqArchR_result)))
    seqArchR_result$rawSeqs <- raw_seqs

If you are not using the result from seqArchR or are already using a curated set of clusters, you can skip the following subsection and jump to section 3.2 for a demonstration of generating all plots.

3.1.2 Curating raw clusters from `seqArchR` result

Raw clusters from seqArchR result (say, the final iteration) often require curation.

The seqArchR result does have a final clustSol (the clustering solution) where the clusters from the final iteration are collated. seqArchR uses hierarchical clustering for this purpose. However, hierarchical clustering, with a chosen agglomeration and distance method, does not necessarily output the best suitable collated set of clusters. Therefore, some small amount of manual curation (re-assignments) may be required to reach the ideal collated set of clusters as the final solution. This is achieved with the help of the seqArchRplus::curate_clusters function.

The basic idea of the curate_clusters function is available in help(seqArchRplus::curate_clusters). A more elaborate description is as follows. It takes as input an agglomeration method and also a distance method (see more arguments in ). On the first call to the function, a plot with the associated dendrogram resulting from the hierarchical clustering is shown together with the per cluster sequence logos. This should help in identifying if the chosen agglomeration and distance methods worked well (i.e., are clusters with similar sequence logos also together in the dendrogram?). If this is satisfactory (i.e., chosen agglomeration and distance method have done well already and only few curations will be required), count the number of overall clusters that you can visually see. Generally, if (estimated) 3-4 (even 5-6) clusters out of 30-40 require curation, I would consider it as a satisfactory outcome from the hierarchical clustering output. Now, in the second call to curate_clusters, specify the number of clusters based on your count. View the resulting clustering in the new plot produced as output. This time the dendrogram shows this clustering with colors and grey-scaled boxes drawn around the clusters. Now you can exactly note which clusters need curation.

Imagine a hypothetical scenario where a total of 16 raw clusters are obtained from the final seqArchR iteration. Collating them into 4 clusters using hierarchical clustering results in:

collated cluster 1: 1, 7, 9
collated cluster 2: 2, 3, 4, 8, 10
collated cluster 3: 11, 12, 13
collated cluster 4: 5, 6, 14, 15, 16

But, due to very minor differences in the sequence logos, you may want to move raw cluster 9 from collated cluster 1 to collated cluster 3 (that of raw clusters 11-13). To do this, set the need_change argument of curate_clusters to c(9) and change_to argument to c(11). Any one of the destination cluster members can be specified in change_to. Similarly, any other such curations can be added. Specify all such curations together in one go like this:

need_change <- list(c(9), c(4, 7), c(16)) and
change_to <- list(11, 14, 0)

Here, (a) raw cluster 9 is moved to the collated cluster containing raw cluster 11; (b) 4 and 7 are moved to that of 14; and (c) 16 is moved to a totally new cluster of itself. The result is the re-assigned clusters.

Below, this procedure of curation is re-iterated with the help of actual function calls and resulting figures using a reduced example data available with the package.


sn <- "RAL28_10_to_12h"
use_aggl <- "complete"
use_dist <- "euclid"

## get seqArchR clusters custom curated
seqArchR_curated_clusts <- curate_clusters(sname = sn,
                                use_aggl = use_aggl, use_dist = use_dist,
                                seqArchR_result = seqArchR_result, iter = 5,
                                pos_lab = NULL, regularize = TRUE, topn = 50,
                                use_cutk = 2, final = FALSE, dir_path = NULL
                            )
#> 
#> ── seqArchR result clusters curation ───────────────────────────────────────────
#> 
#> ── Sample: RAL28_10_to_12h ──
#> 
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
seqArchR_curated_clusts$curation_plot

Figure 1: Figure showing a combined panel of dendrogram + sequence logos at the first step of curation

The motive of Figure 1 is to look at the clusters ordered by the hierarchical clustering step. This way, counting the tentative number of clusters (K) is easier.

In Figure 1, we can see that the first and the second clusters can be collated together. Therefore, to demonstrate curation, we will set K=5 clusters and perform this minor curation as shown below.

## Let us set K=5 for this example, and combine clusters 1 and 2 into one.
set_cutk <- 5

## Form the lists `need_change` and `change_to` for re-assignments
need_change <- list(c(2))
change_to <- list(c(1))

seqArchR_curated_clusts <- curate_clusters(sname = sn,
                                use_aggl = use_aggl, use_dist = use_dist,
                                seqArchR_result = seqArchR_result, iter = 5,
                                regularize = TRUE, topn = 50,
                                use_cutk = set_cutk, #***
                                need_change = need_change,
                                change_to = change_to,
                                final = FALSE, #***
                                dir_path = NULL
                            )
#> 
#> ── seqArchR result clusters curation ───────────────────────────────────────────
#> 
#> ── Sample: RAL28_10_to_12h ──
#> 
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
seqArchR_curated_clusts$curation_plot

Figure 2: Figure showing a combined panel of dendrogram + sequence logos at the second step of curation

Notice that each dendrogram branch (and leaf node) has been assigned a different color because we set K=5. We are now ready to make the final call to curate_clusters and obtain the curated clusters’ list.

## Satisfied with the re-assignments? now set final = TRUE
seqArchR_curated_clusts <- curate_clusters(sname = sn,
                                use_aggl = use_aggl, use_dist = use_dist,
                                seqArchR_result = seqArchR_result, iter = 5,
                                pos_lab = NULL, regularize = FALSE, topn = 50,
                                use_cutk = set_cutk, #***
                                need_change = need_change,
                                change_to = change_to,
                                final = TRUE, #***
                                dir_path = NULL
                            )
#> 
#> ── seqArchR result clusters curation ───────────────────────────────────────────
#> 
#> ── Sample: RAL28_10_to_12h ──
#> 
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
seqArchR_curated_clusts$curation_plot

Figure showing three panels combined: dendrogram + original sequence logos + collated cluster sequence logos at the final step of curation

Figure 3: Figure showing three panels combined: dendrogram + original sequence logos + collated cluster sequence logos at the final step of curation

Notice that the re-assigned clusters (1 and 2) now have different colors than the corresponding dendrogram branches.

str(seqArchR_curated_clusts$clusters_list)
#> List of 4
#>  $ : int [1:244] 1 2 3 4 5 6 7 8 9 10 ...
#>  $ : int [1:111] 245 246 247 248 249 250 251 252 253 254 ...
#>  $ : int [1:93] 356 357 358 359 360 361 362 363 364 365 ...
#>  $ : int [1:142] 449 450 451 452 453 454 455 456 457 458 ...

This fetches us clusters in arbitrary order, i.e., not related to the IQWs of the clusters. [See below for the function call order_clusters_iqw() that orders these clusters by their median/mean interquantile widths (section 3.2.2).] This curated set of clusters can now be used for further downstream analyses.

3.2 Generating individual plots via individual functions

The following subsections demonstrate how the seqArchRplus functions can be used to generate the various visualizations.

3.2.1 Visualize architectures for each promoter cluster as sequence logos

One can visualize the architecture for each cluster of promoters as a sequence logo using the function per_cluster_seqlogos. When the argument one_plot is set to TRUE, the function returns a single plot where sequence logos for all clusters are arranged one below the other. When set to FALSE, a list of sequence logo plots is returned instead. In the later case, each plot has a title with information on the number of sequences in the cluster.

##
seqlogos_oneplot_pl <- per_cluster_seqlogos(
                            sname = "RAL28_10_to_12h",
                            seqs = raw_seqs,
                            clusts = unord_clusts,
                            pos_lab = -45:45, bits_yax = "auto",
                            strand_sep = FALSE, one_plot = TRUE,
                            txt_size = 12, dir_path = NULL
                        )
#> 
#> ── All clusters' sequence logos ────────────────────────────────────────────────
#> 
#> ── Sample: RAL28_10_to_12h ──
#> 
#> Generating architectures for clusters of sequences...
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
#> Scale for y is already present.
#> Adding another scale for y, which will replace the existing scale.
#> Scale for y is already present.
#> Adding another scale for y, which will replace the existing scale.
#> Scale for y is already present.
#> Adding another scale for y, which will replace the existing scale.
#> Scale for y is already present.
#> Adding another scale for y, which will replace the existing scale.
#> Scale for y is already present.
#> Adding another scale for y, which will replace the existing scale.
seqlogos_oneplot_pl

The sequence logos can also be obtained as a list instead of already combined into one plot by simply setting the argument one_plot = FALSE. This gives you liberty to use them (either individually or otherwise) together with any other plots as is suitable.

## Obtain the sequence logos as a list for combining later by setting the 
## 'one_plot' argument to FALSE
seqlogos_list_pl <- per_cluster_seqlogos(
                        sname = "RAL28_10_to_12h",
                        seqs = raw_seqs,
                        clusts = unord_clusts,
                        pos_lab = -45:45, bits_yax = "auto",
                        strand_sep = FALSE, one_plot = FALSE,
                        dir_path = NULL,
                        txt_size = 12
                    )

Obtaining strand-separated sequence logos is demonstrated in section 3.2.7.

3.2.2 Visualize the distribution of IQW and TPM values for each cluster

The function iqw_tpm_plots orders the clusters by their IQW values (when quantitative information on promoter shape is available, for instance, the interquantile widths of the CAGEr tag clusters) and visualizes the distribution of the IQW and TPM values for each cluster as boxplots. The input clusters are ordered from sharp (top) to broad (bottom).

##
iqw_tpm_pl <- iqw_tpm_plots(sname = "RAL28_10_to_12h",
                            dir_path = NULL,
                            info_df = info_df,
                            iqw = TRUE, tpm = TRUE, cons = FALSE,
                            clusts = unord_clusts,
                            txt_size = 15
                            )
#> 
#> ── IQW-ordered boxplots ────────────────────────────────────────────────────────
#> 
#> ── Sample: RAL28_10_to_12h ──
#> 
#> ── Plotting: _IQW _TPM  plots ──
#> 
#> ── Order clusters by IQW ───────────────────────────────────────────────────────
#> 
#> ── Sample: RAL28_10_to_12h ──
#> 
#> ── Order clusters by IQW ───────────────────────────────────────────────────────
#> 
#> ── Sample: RAL28_10_to_12h ──
#> 
iqw_tpm_pl

3.2.3 Fetch clusters ordered by promoter shape

The cluster order visualized in the IQW_TPM plots above can be fetched using the function order_cluster_iqw.

## get clusters ordered by median IQW values
seqArchR_clusts_ord <- order_clusters_iqw(sname = "RAL28_10_to_12h", 
                                                        clusts = unord_clusts,
                                                        info_df = info_df, 
                                                        order_by_median = TRUE
                                    )
#> ── Order clusters by IQW ───────────────────────────────────────────────────────
#> 
#> ── Sample: RAL28_10_to_12h ──
#> 
str(unord_clusts)
#> List of 5
#>  $ : int [1:140] 1 2 3 4 5 6 7 8 9 10 ...
#>  $ : int [1:104] 141 142 143 144 145 146 147 148 149 150 ...
#>  $ : int [1:111] 245 246 247 248 249 250 251 252 253 254 ...
#>  $ : int [1:93] 356 357 358 359 360 361 362 363 364 365 ...
#>  $ : int [1:142] 449 450 451 452 453 454 455 456 457 458 ...
str(seqArchR_clusts_ord)
#> List of 5
#>  $ : int [1:140] 1 2 3 4 5 6 7 8 9 10 ...
#>  $ : int [1:104] 141 142 143 144 145 146 147 148 149 150 ...
#>  $ : int [1:93] 356 357 358 359 360 361 362 363 364 365 ...
#>  $ : int [1:142] 449 450 451 452 453 454 455 456 457 458 ...
#>  $ : int [1:111] 245 246 247 248 249 250 251 252 253 254 ...

3.2.4 Visualize sequence clusters as images colored by the nucleotides

The function seqs_acgt_image enables visualizing a set of sequences as an image colored by the nucleotides. Here, one can choose to order sequences by the clusters (see argument seqs_ord).

##
seqs_acgt_image(sname = "RAL28_10_to_12h",
                        seqs = raw_seqs,
                        seqs_ord = unlist(seqArchR_clusts_ord),
                        pos_lab = -45:45, dir_path = NULL
                    )

seqArchRplus facilitates downstream analysis of clusters of promoter sequence architectures

2024-05-01

Package

Contents

1 Introduction

2 Installation

3 seqArchRplus for downstream analysis of promoter sequences

3.1 Setting up

3.1.1 Input data

3.1.2 Curating raw clusters from `seqArchR` result

3.2 Generating individual plots via individual functions

3.2.1 Visualize architectures for each promoter cluster as sequence logos

3.2.2 Visualize the distribution of IQW and TPM values for each cluster

3.2.3 Fetch clusters ordered by promoter shape

3.2.4 Visualize sequence clusters as images colored by the nucleotides

seqArchRplus facilitates downstream analysis of clusters of promoter sequence architectures

2024-05-01

Package

Contents

1 Introduction

2 Installation

3 seqArchRplus for downstream analysis of promoter sequences

3.1 Setting up

3.1.1 Input data

3.1.2 Curating raw clusters from seqArchR result

3.2 Generating individual plots via individual functions

3.2.1 Visualize architectures for each promoter cluster as sequence logos

3.2.2 Visualize the distribution of IQW and TPM values for each cluster

3.2.3 Fetch clusters ordered by promoter shape

3.2.4 Visualize sequence clusters as images colored by the nucleotides

3.1.2 Curating raw clusters from `seqArchR` result