MURASAKI(1)

๐Ÿ’– MURASAKI ๐Ÿ’–

okay bestie, welcome to the murasaki manual page (that's the MURASAKI(1) one, user contributed Perl documentation, iykyk) โ€” it computes anchors between multiple sequences and honestly? it's giving main character energy. buckle up, this is part one of three, so don't go anywhere. ๐Ÿฆ‹


โœจ NAME โœจ

murasaki โ€” compute anchors between multiple sequences


๐ŸŒธ SYNOPSIS ๐ŸŒธ

murasaki [OPTIONS] -p[pattern] seq1.fa seq2.gbk [seq3.raw ...] #compute anchors between seq1.fa and seq2.gbk using [pattern]
mpirun murasaki [OPTIONS] -p[pattern] seq1.fa seq2.gbk [seq3.raw ...] #compute anchors between seq1.fa and seq2.gbk using [pattern] in parallel via MPI

๐ŸŽ€ DESCRIPTION ๐ŸŽ€

okay so here's the tea โ˜• โ€” Murasaki generates anchors based on all the sequences you hand it, using the pattern you supply plus hash tables. Basically, every single base of every single sequence gets masked by the pattern, which forms a seed, and that seed gets used to generate a hash. The location of the seed? Stored in the hash table, obviously. ๐Ÿ’…

Once all the seeds have been hashed and stored, Murasaki scans the hash table and generates anchors for all the matching seeds. An anchor is a set of intervals across a subset of the input sequences โ€” no cap. These live in name.anchors files, and they're described in "FILE FORMATS".

By default, anchors are maximally extended until their minimum pairwise ungapped alignment score drops below a threshold โ€” same fashion as the X-drop parameter in BLAST and BLAST-like searches. literally the same vibe. ๐Ÿ”ฅ


๐Ÿงฌ PATTERNS ๐Ÿงฌ

Murasaki uses spaced seed patterns when it's considering seeds. A spaced seed pattern is typically written as a string of 1s and 0s, and it necessarily starts and ends with a 1 (that part is not optional bestie, it's the rules).

  • 1s mean this base is considered part of the seed.
  • bases sitting at 0 positions are not.

For example, with a pattern "1011", the sequence "ACGT" would match "AGGT" and "ATGT" but not "ACTT". slay.

The number of 1s in the pattern is called the "weight" of the pattern, and the number of 1s and 0s combined is the "length" of the pattern. Murasaki lets you use any arbitrary pattern expressed as a string of 1s and 0s, and it also interprets patterns of the form "x:y" to mean a "random pattern of weight x and length y."

The choice of pattern obviously has an impact on sensitivity and specificity, but whether one pattern is "better" than another depends on the application and the input sequences under consideration. Calculating "maximally sensitive spaced seed patterns" is a computationally difficult problem, and there are a number of research papers describing various methods for approximation (peep "RELATED READING"). In general though? "Heavier" spaced seed patterns are less sensitive but more specific than lighter seeds. Anecdotally we find that seeds with weights approximately 60% to 75% (with lengths around 24 for bacteria, and 36 to 48 for mammals) are good for most applications. Extremely similar species (for example human and chimp) benefit from longer, heavier seeds. โœจ


๐Ÿ” HASH FUNCTIONS ๐Ÿ”

Hash functions (and hash parameters too) are generated automatically, based on the system environment and the input sequences. There are essentially two types of hash functions available in Murasaki: adaptive and cryptographic hashes.

  • The adaptive hashes are XOR combinations of various bitwise shifts of the seed, designed by analyzing the spaced seed pattern to maximize the entropy of the resulting hash.
  • The cryptographic hashes are available via the CryptoPP library and use the entire spaced seed pattern to generate a hash using one of the common cryptographic hashes like MD5 or SHA-1.

The adaptive hash functions are almost always faster and more efficient than MD5 and SHA-1, but the cryptographic functions are available for reference and may be useful as an alternative in the unlikely event you're dealing with an environment where the adaptive hasher is unsuitable (for example a sequence consisting of only A and T (leaving 1 out of every 2 bits uninitialized)). ๐Ÿ’€


๐Ÿง  MEMORY SCALING ๐Ÿง 

Murasaki can take a lot of memory. Storing the location of each seed in the hash table is the most costly part of the operation, requiring approximately ceil(log_2(N)) bits per seed where N is the total sequence length. Locations are, by default, stored in a bitpacked format to approach theoretical minimum. so efficient, we love her. ๐Ÿซถ

The second most costly element is the hash table structure, where each bucket carries a small overhead and unused [buckets] are simply wasted space. More hash table buckets (i.e. a longer hash table) decreases the expected number of collisions, leading to faster execution time. Therefore Murasaki tries to use as many buckets as possible by inspecting the available system memory and using as much as it can while still storing all the seed locations.

If this automatic scaling is ineffective, setting the hash table size directly via the --hashbits|-b options can force a specific hash table size. And if the memory of one computer is insufficient to store the desired hash table, PARALLELIZATION can be used to distribute the hash table across multiple computers. ๐Ÿ‘‘

๐ŸŒˆ PARALLELIZATION ๐ŸŒˆ

Murasaki is designed to run in parallel using MPI. Consult the documentation for the specific variations of your MPI implementation, however in general the execution method looks like:

mpirun [MPI options] murasaki [murasaki options] -p[pattern] [seq1 ...]

Murasaki in parallel divides the number of processors available (NP) into two groups: hasher nodes and storage nodes.

  • The storage nodes divide the hash table between themselves, each one being responsible for a different part of the table.
  • The hasher nodes divide the input sequence between themselves, each hashing a separate portion of the input sequence, and passing the seed location to the appropriate storage node for storage.

When all the hasher nodes are finished hashing, the storage nodes scan their portion of the hash table and pass matching sets of seeds to a hasher node, where they're assembled into anchors and extended. Finally, all the hasher nodes combine their independent anchor sets into one final set in ceil(log_2(H)) iterations (where H is the number of hasher nodes), with each hasher node number 2h passing its anchors to hasher number 2h-1 at each iteration. ๐ŸŒŸ

Because almost none of the parallelization steps require communication between all nodes, and each seed and each anchor can be processed in parallel, Murasaki scales very well in parallel โ€” running approximately twice as fast when twice as many nodes are available. fr fr. Furthermore, the hash table is automatically grown to take advantage of the combined memory from multiple machines. iconic behavior honestly. ๐Ÿ’œ


โœจ OPTIONS โœจ

okay so here's the thing bestie: most options can be written out in their long form โ€” like --directory out or --directory=out โ€” or you can go short, like -dout. both work, pick your vibe ๐Ÿ’…

and there's a little alphabet soup for the value types, so let's decode it:

  • <S> โ€” wants a string
  • <D> โ€” wants an integer
  • <F> โ€” wants a float
  • <B> โ€” wants a boolean: yes/on/true/1 for true, no/off/false/0 for false

and the booleans? most of them can just omit the value entirely and they'll flip from whatever they were to the opposite. toggling queen ๐Ÿ‘‘


murasaki hasโ€ฆ a lot of options. like, a lot a lot. so they're sorted into categories here to help you tell what each one's actually for. but heads up: in certain situations certain option choices may have unforeseen consequences, and obviously if the output you asked for is huge, the runtime is gonna be long. it is a mistake to think that everything outside the "tuning options" in the Performance section has no bearing on performance. no cap. ๐Ÿซถ

these ones shape what even counts as an "anchor".

  • --pattern|-p <S> โ€” sets the seed pattern (eg. 11101001010011011). if you use the format [<w>:<l>] it automatically generates a random pattern of weight <w> and length <l> โœจ
  • --repeatmask|-r <B> โ€” skip repeat masked data (ie: lowercase atgc). just be aware that some sequence files are distributed purely in lower case, soโ€ฆ vibe check your input first ๐Ÿ˜ญ
  • --seedfilter|-f <D> โ€” skip seeds that occur more than N times. exceptionally slow. see --hashfilter for a faster approximation.
  • --hashfilter|-m <D> โ€” like --seedfilter but it works on hash keys instead of seeds. may cause some collateral damage to otherwise unique seeds, but it's faster. tradeoffs, bestie ๐Ÿ’€
  • --skipfwd|-F <B> โ€” don't hash/match the forward strands.
  • --skiprev|-R <B> โ€” don't hash/match the reverse complement strands.
  • --skip1to1|-1 <B> โ€” skip matches along the 1:1 line (good for comparing to self ๐Ÿ‘€).
  • --hashonly|-Q <B> โ€” Hash Only. no anchor output, just statistics.
  • --hashskip|-S <D> โ€” hashes every n bases. default is 1 (i.e. hashing all positions). not supplying any argument increments the skip amount by 1.
  • --join|-j <D> โ€” join anchors within n bases of each other (default: 0). specifying a negative D implies -D*patternLength.
  • --bitscore|-B <B> โ€” toggles computation of a bitscore for all anchors (default is on).
  • --seedterms|-T <B> โ€” toggles retention of seed terms (defaults to off). these are necessary for computing TF-IDF scores).
  • --sectime|-e <B> โ€” always display times in seconds as opposed to human readable "1d 3h 45m 5s" style times.
  • --mergefilter|-Y <D> โ€” filter out matches which would cause more than D many anchors to be generated from 1 seed (default -Y100). use -Y0 to disable.
  • --scorefilter <D> โ€” set a minimum ungapped score for seeds.
  • --rifts|-/ <D> โ€” allow anchors to skip D sequences (default 0).
  • --islands|-% <D> โ€” same as --rifts=S-D (where S is number of input seqs).
  • --fuzzyextend|-z <B> โ€” enable (default) or disable fuzzy extension (i.e. ungapped alignment) of anchors.
  • --fuzzyextendlosslimit|-Z <D> โ€” set the cutoff at which to stop extending fuzzy hits (ie. the BLAST X parameter).
  • --gappedanchors <B> โ€” use gapped (true) or ungapped (false (default)) anchors.
  • --scorebyminimumpair <B> โ€” do anchor scoring by minimum pair when appropriate (default). the alternative is arithmetic mean (seldom useful, but theoretically faster).

๐ŸŒธ Output options

these are the ones that mostly decide what data goes where.

  • --directory|-d <S> โ€” output directory (default: output)
  • --name|-n <S> โ€” alignment name (default: test)
  • --repeatmap|-i <B> โ€” toggles keeping of a repeat map when --mergefilter is used (defaults to yes).
  • --histogram|-H <D> โ€” histogram computation level (-H alone implies -H1):
    • 0 - no histogram (default)
    • 1 - basic bucketsize/bucketcount histogram data
    • 2 - bucket-based scores to anchors.detils
    • 3 - perbucket count data
    • 4 - perbucket + perpattern count data any values above 2 are purely explorartory and can result in massive output files. like, unhinged file sizes. be careful ๐Ÿ’€
  • --tfidf|-k <B> โ€” perform accurate tfidf scoring from within murasaki (requires extra memory at anchor generation time). default is no.

๐Ÿ”ฅ Performance/tuning options

these mostly affect performance, and don't (in general) impact output. in general, okay, don't @ me.

  • --quickhash|-q <D> โ€” specify a hashing function:
    • 0 - adaptive with S-boxes (default when there's plenty of hash table to spare)
    • 1 - don't pack bits to make hash (use first word only)
    • 2 - naively use the first hashbits worth of pattern
    • 3 - adaptivevely find a good hash (default)
    • experimental CryptoPP hashes
    • 4 - MD5
    • 5 - SHA1
    • 6 - Whirlpool
    • 7 - CRC-32
    • 8 - Adler-32
Note: 3 and 0 are the only "recommended" hash functions, and the only ones automatically selected. the others are provided merely for reference. 1, 7, and 8 aren't even expected to utilize the entire hash space. ๐Ÿ™ˆ
  • --hashbits|-b <D> โ€” use D bit hashes (for n's of 1 to WORDSIZE. default 26)
  • --hashtype|-t <S> โ€” select hash table data structure to use:
    • OpenHash - open sub-word packing of hashbits (default when there's plenty of hash table to spare)
    • EcoHash - chained sub-word packing of hashbits (default)
    • ArrayHash - malloc/realloc (fast but fragmentation-prone)
    • MSetHash - memory exorbanant, almost pointless.
  • --probing <D> โ€” 0 - linear, 1 - quadratic (default). only applicable for --hashtype=OpenHash.
  • --hitfilter|-h <D> โ€” minimum number of hits to be outputted as an anchor (default 1). in PatternHunter this is 2.
  • --rseed|-s <D> โ€” random number seed for non-deterministic algorithms (ie: adative hash function generation). if you're doing any performance comparisons, it's probably imperative that you use the same seed for each run of the same settings. default is obtained from time() (ie: seconds since 1970).
  • --memory|-M [<F>|<S>] โ€” set the target amount of total memory (either in gb or as % total memory).
  • --reverseotf|-o <B> โ€” generate reverse complement on the fly (defaults to on). turning this off precomputes all the reverse complement strands and stores them in memory, which rarely provides a measurable performance improvement.
  • --binaryseq <B> โ€” enable (default) or disable binary sequence read/write

performance options related to adaptive hash function generation.

  • --hasherFairEntropy <B> โ€” use more balanced entropy estimation (default: yes).
  • --hasherCorrelationAdjust <B> โ€” adjust entropy estimates for nearby sources assuming some correlation (default: yes).
  • --hasherTargetGACycles <D> โ€” adaptive hash function generation genetic algorithm cycle cutoff.
  • --hasherEntropyAgro <F> โ€” how aggressive to be about pursuing maximum entropy hash functions (takes a real. default is 1).

๐Ÿ’œ MPI Specific

  • --hashers|-A [<F>|<D>] โ€” specify the number of processes to be used as hashers (only applies to MPI. if a number between 0 and 1 it refers to a ratio of np).
  • --localhash|-K <B> โ€” perform hashing locally on each storage node rather than sending it over the network (helpful for slow networks).
  • --mpidistro|-L <B> โ€” toggles use of MPI to distribute sequence data over (if the sequence is available on local disk on each node then turning this off may potentially accerlate the initial sequence loading).
  • --waittoanchor|-w <B> โ€” postpone actual anchor computation until all location sets have been received (as opposed to trying to work between receiving seed packets).
  • --buffers|-u <D> โ€” maximum number of unfinished buffers to allow while message passing (0 means unlimited). default is set based on the number of nodes participating.
MPI can crash or perform very poorly if this value is too high. โš ๏ธ
  • --nobuffers|-U <B> โ€” same as --buffers=1.
  • --bigfirst|-I <B> โ€” assign hashers to large memory nodes first.
  • --hostbalance|-l <B> โ€” if yes (default): spread out hashers evenly among all nodes. if no: ignore host name when assigning jobs.
  • --memorybalance|-a <B> โ€” if yes (default): balance hash storage between nodes based on the amount of available ram. if no: distribute storage evently. this more likely to achieve optimal run times, but might not utilize memory as efficiently.
  • --distmerge|-< <B> โ€” if yes (default): during the merge step, storage nodes send seeds to any available hasher. if no: send all seeds to one node only.
  • --distcollect|-> <B> โ€” if yes (default): collect anchor data from all hashers. if no: send all seeds to the final assembly node only.
  • --mpiredirectoutput <B> โ€” if yes (default): each rank redirects its stdout/stderr to a separate file (murasaki-mpiout-N). if no: do what comes naturally (ie: managed by mpirun (for OpenMPI see --output-filename and --tag-output in mpirun(1))).
  • --keepstdoe <B> โ€” don't erase the murasaki-mpiout files on success.
  • --sysvipc|-V <B> โ€” use System V IPC to negotiate shared memory regions (saves memory when one host runs multiple nodes). default is true.

โญ Universal options

  • --verbose|-v โ€” increases verbosity.
  • --version|-V โ€” prints version information and quits.
  • --help|-? โ€” prints a help message and quits.

๐Ÿ“‚ FILE FORMATS ๐Ÿ“‚

Murasaki drops a whole array of output files, and honestly bestie, most of the formats are designed to just ~make sense~ on sight. Every output file gets prefixed with whatever you gave to the --name parameter, okay? The big ones are described right here. Unless it says otherwise, files are line based and tab delimited. ๐Ÿ’–


๐Ÿงฌ .seqs

The .seqs file shows which sequences got used as input โ€” one per line, periodt. Various programs use this file together with the .anchors file, so it's generally important that the contents actually reflect the correct sequence files. ๐Ÿ‘‘

Heads up: if you move anchor results between computers, the paths might change and you'll have to update the .seqs file yourself. An alternative that lowkey fixes this: always use relative paths. ๐Ÿซถ

โš“ .anchors files

These are 1 anchor per line, with a 3-tuple per sequence. Each tuple = the start and stop coordinates plus the strand of the anchored interval on that sequence. The sequence order matches the order in the .seqs file. ๐ŸŒธ

Coordinates are structured so that 1 is the first base in the sequence, 2 is the second, etc. Negative values refer to the reverse complement sequence, where -1 is the last base of the reverse complement sequence (ie: the complement first base in the forward sequence). The "strand" element is a '+' or '-' that merely matches the sign of the coordinates โ€” yes, that's redundant info, but it's kept around to make parsing or filtering simpler. ๐Ÿ’…

For examle:

 1       18     +       -1      -18       -

This line describes an anchor where the first 18 bases of the first sequence match the first 18 bases of the reverse complement of the second sequence. Slay. โœจ

๐Ÿ“œ .anchors.details

This is an antiquated file format, ngl, but GMV uses it to calculate statistics like TF-IDF scores, so it's been kept around for that reason. The .anchors.details file has the same format and information as the .anchors file โ€” except after the anchor tuples come two more terms: a score, and a comma (,) delimited list of term and count pairs (written "term:count"). The score and count data might be varied depending on the --histogram option choices. ๐ŸŽ€

๐Ÿงฎ .anchors.bitscore

The term "bitscore" here is a misnomer, but it's maintained for historical reasons. In reality this file contains the mean number of matching bases and the length of each anchor (corresponding line by line to the .anchors file). ๐Ÿ’€

๐Ÿ“Š .stats.tfidf

Contains anchor TF-IDF scores (corresponding line by line to the .anchors file). Cute, simple, iconic. ๐Ÿ’œ

๐Ÿ“ˆ .histogram

Contains a simple histogram of the hash table usage. First field = the bucket size, second field = the frequency. So for example, a .histogram file like this:

 1  24
 2  1

Would indicate that there were 24 hash buckets that stored only 1 location (i.e. 24 unique seeds), and 1 hash bucket stored 2 locations (i.e. 1 seed that matched 2 locations (or 2 non-matching seeds that resulted in a hash collision)). ๐ŸŒˆ

โš™๏ธ .options

Maintains a record of the options used when running Murasaki. Receipts, basically. ๐Ÿงพ

๐Ÿ” .repeats

The .repeats file stores a record of "repeats" as defined by the --mergefilter option (i.e. seeds that would have have induced more anchors than permitted). In this file, each repeat record is separated by a blank line. A repeat record looks like this:

 R: G.GCCTTT.T.ACT.CACAA..AT
 0: 2145540494 -425039256 -113794380 1998323403
 1: 2480929222 -1874514626 2543723555 -2550045172

The first line (always prefixed "R:") shows the repeating seed itself (where the . are the bases masked by the pattern). The subsequent lines show where these seeds occurred in the input sequences (in the first (0) and second (1) sequences). Note that if there are no hits in a particular sequence, it doesn't include a blank line for that sequence. For example:

 R: G.GCCTTT.T.ACT.CACAA..AT
 0: 2145540494 -425039256 -113794380 1998323403
 2: 2480929222 -1874514626 2543723555 -2550045172

is also a valid .repeats file. ๐Ÿฆ‹


๐Ÿ“„ LICENSE ๐Ÿ“„

GNU General Public License, version 3 (GPLv3). Free as in bestie. ๐Ÿ’–

๐ŸŒ AVAILABILITY ๐ŸŒ

<http://murasaki.sourceforge.net>

๐Ÿ‘‘ AUTHOR ๐Ÿ‘‘

Kris Popendorf <krisp@dna.bio.keio.ac.jp>

๐Ÿ‘€ SEE ALSO ๐Ÿ‘€

mbfa(1), geneparse(1)

  • M. Csuros and B. Ma, "Rapid Homology Search with Two-Stage Extension and Daughter Seeds" (2005).
  • F. P. Preparata and L. Zhang and K. W. Choi, "Quick, practical selection of effective seeds for homology search" (2005).
  • KP Choi, et. al., "Good spaced seeds for homology search" (2004).

And that's a wrap, bestie!! ๐ŸŽ€ This page is perl v5.10.1, dated 2010-05-31, and it's the man page for MURASAKI(1) โ€” go forth and align some sequences, you absolute icon. ๐Ÿ’…โœจ