Introduction

seqArchR is a non-negative matrix factorization (NMF)-based unsupervised learning approach for identifying different core promoter sequence architectures. seqArchR implements an algorithm based on chunking and iterative processing. While matrix factorization-based applications are known to scale poorly for large amounts of data, seqArchR’s algorithm enables scalable processing of large number of sequences. A notable advantage of seqArchR is that the sequence motifs – the lengths and positional specificities of individual motifs, and complex inter-relationships where multiple motifs are at play in tandem, all are simultaneously inferred from the data. To our knowledge, this is a novel application of NMF on biological sequence data capable of simultaneously discovering the sequence motifs and their positions. For a more detailed discussion, see preprint/publication.

This vignette demonstrates seqArchR’s usage with the help of a synthetic DNA sequences data set. Please refer to the paper for a detailed description of seqArchR’s algorithm. The paper also discusses the various parameters and their settings. For completeness, the following section gives a brief overview of the algorithm.

seqArchR’s algorithm

seqArchR implements a chunking-based iterative procedure. Below is a schematic of seqArchR’s algorithm.

Further details to follow.

Installation

Python scikit-learn dependency

seqArchR requires the Python module scikit-learn. Please see installation instructions here.

seqArchR is available on Bioconductor, and can be installed using:


if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("seqArchR")

In case of any errors, please consider looking up: https://github.com/snikumbh/seqArchR. If none of the already noted points with regards to troubleshooting seqArchR’s installation help, please file a new issue.

Working with seqArchR

# Load seqArchR
library(seqArchR)
library(Biostrings, quietly = TRUE)
#> 
#> Attaching package: 'generics'
#> The following objects are masked from 'package:base':
#> 
#>     as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
#>     setequal, union
#> 
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:stats':
#> 
#>     IQR, mad, sd, var, xtabs
#> The following objects are masked from 'package:base':
#> 
#>     Filter, Find, Map, Position, Reduce, anyDuplicated, aperm, append,
#>     as.data.frame, basename, cbind, colnames, dirname, do.call,
#>     duplicated, eval, evalq, get, grep, grepl, is.unsorted, lapply,
#>     mapply, match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
#>     rank, rbind, rownames, sapply, saveRDS, table, tapply, unique,
#>     unsplit, which.max, which.min
#> 
#> Attaching package: 'S4Vectors'
#> The following object is masked from 'package:utils':
#> 
#>     findMatches
#> The following objects are masked from 'package:base':
#> 
#>     I, expand.grid, unname
#> 
#> Attaching package: 'Biostrings'
#> The following object is masked from 'package:base':
#> 
#>     strsplit


# Set seed for reproducibility
set.seed(1234)

Synthetic data explained

In order to demonstrate the efficacy of seqArchR, we use seqArchR to cluster DNA sequences in a synthetic data set which was generated as follows. A set of 200 simulated DNA sequences was generated, each 100 nucleotides long and with uniform probability for all nucleotides. These sequences have four clusters in them, each with 50 sequences. The profiles of the four clusters are:

Cluster	Characteristic Motifs	Motif Occurrence Position	#Sequences
A	Dinucleotide repeat `AT`	every 10 nt	50
B	`GATTACA`	40	50
	`GAGAG`	60
C	`GAGAG`	60	50
D	`GAGAG`	80	50
	`TCAT`	40

All the motifs across the clusters were planted with a mutation rate of 0.

Input and feature representation

We use one-hot encoding to represent the dinucleotide profiles of each sequence in the data set. seqArchR provides functions to read input from (a) a FASTA file, and (b) Biostrings::DNAStringSet object.

Reading input as FASTA file

The function seqArchR::prepare_data_from_FASTA() enables one-hot-encoding the DNA sequences in the given FASTA file. The one-hot-encoded sequences are returned as a sparse matrix with as many columns as the number of sequences in the FASTA file and (sequence length x 4²) rows when dinucleotide profiles is selected. The number of rows will be (sequence length x 4) when mononucleotide profiles is selected. See the sinuc_or_dinuc argument.

Upon setting the logical argument rawSeq to TRUE, the function returns the raw sequences as a Biostrings::DNAStringSet object, with FALSE it returns the column-wise one-hot encoded representation as noted above. When raw_seq is TRUE, sinuc_or_dinuc argument is ignored.

# Creation of one-hot encoded data matrix from FASTA file
inputFname <- system.file("extdata", "example_data.fa.gz", 
                                  package = "seqArchR", 
                                  mustWork = TRUE)

# Specifying `dinuc` generates dinucleotide features
inputSeqsMat <- seqArchR::prepare_data_from_FASTA(fasta_fname = inputFname,
                                                  sinuc_or_dinuc = "dinuc")
#> Sequences OK,
#> Read 200 sequences
#> Generating dinucleotide profiles

inputSeqsRaw <- seqArchR::prepare_data_from_FASTA(fasta_fname = inputFname, 
                                               raw_seq = TRUE)

nSeqs <- length(inputSeqsRaw)
positions <- seq(1, Biostrings::width(inputSeqsRaw[1]))

Reading input as a DNAStringSet object

If you already have a Biostrings::DNAStringSet object, you can use the seqArchR::get_one_hot_encoded_seqs() function which directly accepts a Biostrings::DNAStringSet object.

# Creation of one-hot encoded data matrix from a DNAStringSet object
inputSeqs_direct <- seqArchR::get_one_hot_encoded_seqs(seqs = inputSeqsRaw, 
                                                  sinuc_or_dinuc = "dinuc")
#> Generating dinucleotide profiles

identical(inputSeqs_direct, inputSeqsMat)
#> [1] TRUE

Visualize input sequences as an image

# Visualize the sequences in a image matrix where the DNA bases are 
# assigned fixed colors

seqArchR::viz_seqs_acgt_mat(as.character(inputSeqsRaw), 
                          pos_lab = positions, save_fname = NULL)

Calling seqArchR

Setup seqArchR configuration as follows.

# Set seqArchR configuration
seqArchRconfig <- seqArchR::set_config(
        parallelize = TRUE,
        n_cores = 2,
        n_runs = 100,
        k_min = 1,
        k_max = 20,
        mod_sel_type = "stability",
        bound = 10^-6,
        chunk_size = 100,
        result_aggl = "ward.D", 
        result_dist = "euclid",
        flags = list(debug = FALSE, time = TRUE, verbose = TRUE,
                     plot = FALSE)
)

Once the configuration is setup, call the seqArchR::seqArchR() function with user-specified number of iterations.

# Call/Run seqArchR
seqArchRresult <- seqArchR::seqArchR(config = seqArchRconfig,
                            seqs_ohe_mat = inputSeqsMat,
                            seqs_raw = inputSeqsRaw,
                            seqs_pos = positions,
                            total_itr = 2,
                            set_ocollation = c(TRUE, FALSE))

Understanding the result object from seqArchR

In the version 1.11.0, seqArchR naively returns a result object which is a nested list of seven elements. These include:

the sequence cluster labels per iteration [seqsClustLabels];
the collection of NMF basis vectors per iteration [clustBasisVectors]: each is a list of two elements nBasisVectors and basisVectors;
the clustering solution, [clustSol], which is obtained upon combining raw clusters from the last iteration of seqArchR. This element stores the clustering of NMF basis vectors [basisVectorsClust] and the sequence clusters [clusters];
the raw sequences provided [rawSeqs];
if timeFlag is set, timing information (in minutes) per iteration [timeInfo];
the configuration setting [config]; and
the call itself [call].

NMF basis vectors

seqArchR stores the NMF basis vectors corresponding to each cluster in every iteration in the variable clustBasisVectors. clustBasisVectors is a numbered list corresponding to the number of iterations performed. This is then again a list holding two pieces of information: the number of basis vectors (nBasisVectors) and the basis vectors (basisVectors).


# Basis vectors at iteration 2
seqArchR::get_clBasVec_k(seqArchRresult, iter=2)
#> [1] 4

i2_bv <- seqArchR::get_clBasVec_m(seqArchRresult, iter=2)
dim(i2_bv)
#> [1] 1600    4
head(i2_bv)
#>             [,1]       [,2]       [,3]       [,4]
#> [1,] 0.120740995 0.07407753 0.05168607 0.08186580
#> [2,] 0.087066440 0.08679174 0.02486037 0.07917296
#> [3,] 0.080389797 0.03793642 0.05076204 0.05196331
#> [4,] 0.074610342 0.04276400 0.09945618 0.05334238
#> [5,] 0.080347189 0.11747636 0.07484865 0.10758754
#> [6,] 0.007201571 0.11872570 0.05051104 0.10543655

The NMF basis vectors can be visualized as a heatmap and/or sequence logo using viz_bas_vec_heat_seqlogo function.

Basis vectors at iteration 1

seqArchR::viz_bas_vec(feat_mat = get_clBasVec_m(seqArchRresult, 1), 
                      ptype = c("heatmap", "seqlogo"), method = "bits", 
                      sinuc_or_dinuc = "dinuc")
#> Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
#> of ggplot2 3.3.4.
#> ℹ The deprecated feature was likely used in the ggseqlogo package.
#>   Please report the issue at <https://github.com/omarwagih/ggseqlogo/issues>.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
#> 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.
#> [[1]]

#> 
#> [[2]]

#> 
#> [[3]]

Basis vectors at iteration 2

seqArchR::viz_bas_vec(feat_mat = get_clBasVec_m(seqArchRresult, 2), 
                      ptype = c("heatmap", "seqlogo"), method = "bits", 
                      sinuc_or_dinuc = "dinuc")
#> 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.
#> [[1]]

#> 
#> [[2]]

#> 
#> [[3]]

#> 
#> [[4]]

Visualize sequences by clusters

The clustered output from seqArchR can again be visualized as a matrix. Use the https://snikumbh.github.io/seqArchR/reference/seqs_str.html function to fetch sequences by clusters at any iteration and call seqArchR::viz_seqs_as_acgt_mat as shown.


seqArchR::viz_seqs_acgt_mat(seqs_str(seqArchRresult, iter = 1, ord = TRUE),
                                  pos_lab = positions)

Clusters at iteration 1


seqArchR::viz_seqs_acgt_mat(seqs_str(seqArchRresult, iter = 2, ord = TRUE),
                                  pos_lab = positions)

Clusters at iteration 2

Conclusion

seqArchR can detect de novo sequence features and simultaneously identify the complex interactions of different features together with their positional specificities.

Note that the sequence architectures identified by seqArchR have no limitations due to the size of the motifs or gaps in them, distance between motifs, compositional and positional variations in the individual motifs and their effects on the complex interactions, and number of motifs involved in any interaction.

Session Info

sessionInfo()
#> R version 4.4.2 (2024-10-31)
#> 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] Biostrings_2.75.3   GenomeInfoDb_1.43.2 XVector_0.47.1     
#> [4] IRanges_2.41.2      S4Vectors_0.45.2    BiocGenerics_0.53.3
#> [7] generics_0.1.3      seqArchR_1.11.0     BiocStyle_2.35.0   
#> 
#> loaded via a namespace (and not attached):
#>  [1] sass_0.4.9              stringi_1.8.4           lattice_0.22-6         
#>  [4] digest_0.6.37           magrittr_2.0.3          evaluate_1.0.1         
#>  [7] grid_4.4.2              fastmap_1.2.0           plyr_1.8.9             
#> [10] jsonlite_1.8.9          Matrix_1.7-1            BiocManager_1.30.25    
#> [13] httr_1.4.7              UCSC.utils_1.3.0        scales_1.3.0           
#> [16] codetools_0.2-20        jquerylib_0.1.4         cli_3.6.3              
#> [19] rlang_1.1.4             crayon_1.5.3            cowplot_1.1.3          
#> [22] munsell_0.5.1           withr_3.0.2             cachem_1.1.0           
#> [25] yaml_2.3.10             tools_4.4.2             ggseqlogo_0.2          
#> [28] parallel_4.4.2          reshape2_1.4.4          BiocParallel_1.41.0    
#> [31] colorspace_2.1-1        ggplot2_3.5.1           GenomeInfoDbData_1.2.13
#> [34] buildtools_1.0.0        vctrs_0.6.5             R6_2.5.1               
#> [37] lifecycle_1.0.4         stringr_1.5.1           pkgconfig_2.0.3        
#> [40] bslib_0.8.0             pillar_1.10.0           gtable_0.3.6           
#> [43] Rcpp_1.0.13-1           glue_1.8.0              xfun_0.49              
#> [46] tibble_3.2.1            sys_3.4.3               knitr_1.49             
#> [49] farver_2.1.2            htmltools_0.5.8.1       labeling_0.4.3         
#> [52] rmarkdown_2.29          maketools_1.3.1         compiler_4.4.2         
#> [55] prettyunits_1.2.0

Example usage of seqArchR on simulated DNA sequences