Example usage of seqArchR on simulated DNA sequences
2023-04-25
Package
seqArchR 1.4.0
Contents
1 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.
2 seqArchR’s algorithm
seqArchR implements a chunking-based iterative procedure. Below is a schematic of seqArchR’s algorithm.
Further details to follow.
3 Installation
3.1 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.
4 Working with seqArchR
# Load seqArchR
library(seqArchR)
library(Biostrings, quietly = TRUE)
#>
#> 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, intersect, is.unsorted,
#> lapply, mapply, match, mget, order, paste, pmax, pmax.int, pmin,
#> pmin.int, rank, rbind, rownames, sapply, setdiff, sort, table,
#> tapply, union, 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)4.1 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.
4.2 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.
4.2.1 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^{2}\)) 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]))
4.2.2 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] TRUE4.3 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)
4.4 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))
4.5 Understanding the result object from seqArchR
In the version 1.4.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 elementsnBasisVectorsandbasisVectors; - 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].
4.5.1 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.10543655The NMF basis vectors can be visualized as a heatmap and/or sequence logo using viz_bas_vec_heat_seqlogo function.
4.5.2 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]]
4.5.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]]
4.5.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)
Figure 1: Clusters at iteration 1
seqArchR::viz_seqs_acgt_mat(seqs_str(seqArchRresult, iter = 2, ord = TRUE),
pos_lab = positions)
Figure 2: Clusters at iteration 2
5 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.
6 Session Info
sessionInfo()
#> R version 4.3.0 RC (2023-04-13 r84269)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 22.04.2 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.17-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB 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: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] Biostrings_2.68.0 GenomeInfoDb_1.36.0 XVector_0.40.0
#> [4] IRanges_2.34.0 S4Vectors_0.38.0 BiocGenerics_0.46.0
#> [7] seqArchR_1.4.0 BiocStyle_2.28.0
#>
#> loaded via a namespace (and not attached):
#> [1] gtable_0.3.3 xfun_0.39 bslib_0.4.2
#> [4] ggplot2_3.4.2 lattice_0.21-8 vctrs_0.6.2
#> [7] tools_4.3.0 bitops_1.0-7 generics_0.1.3
#> [10] parallel_4.3.0 tibble_3.2.1 fansi_1.0.4
#> [13] highr_0.10 pkgconfig_2.0.3 Matrix_1.5-4
#> [16] lifecycle_1.0.3 GenomeInfoDbData_1.2.10 compiler_4.3.0
#> [19] farver_2.1.1 stringr_1.5.0 munsell_0.5.0
#> [22] codetools_0.2-19 htmltools_0.5.5 sass_0.4.5
#> [25] RCurl_1.98-1.12 yaml_2.3.7 pillar_1.9.0
#> [28] crayon_1.5.2 jquerylib_0.1.4 BiocParallel_1.34.0
#> [31] cachem_1.0.7 magick_2.7.4 tidyselect_1.2.0
#> [34] digest_0.6.31 stringi_1.7.12 dplyr_1.1.2
#> [37] reshape2_1.4.4 bookdown_0.33 labeling_0.4.2
#> [40] cowplot_1.1.1 fastmap_1.1.1 grid_4.3.0
#> [43] colorspace_2.1-0 cli_3.6.1 magrittr_2.0.3
#> [46] utf8_1.2.3 withr_2.5.0 prettyunits_1.1.1
#> [49] scales_1.2.1 rmarkdown_2.21 evaluate_0.20
#> [52] knitr_1.42 rlang_1.1.0 Rcpp_1.0.10
#> [55] glue_1.6.2 BiocManager_1.30.20 jsonlite_1.8.4
#> [58] R6_2.5.1 plyr_1.8.8 zlibbioc_1.46.0
#> [61] ggseqlogo_0.1