Seq – A Language for Bioinformatics

As a former bioinformatician (if that’s a word), I’m not sure there’s much value in this. There’s high dispersion in the performance requirements of bioinformatics tools. The processes that need to be fast (alignment, BLAST, whatever, tree creation, etc) are already super fucking optimized (though, unfortunately, still slow). The things that don’t need to be fast can use whatever you want (I used Haskell and Racket for my own tools at the time). Python is.. not the greatest. The major value add is the multitude of scientific libraries. If you’re gonna throw that all away, why not just use something better? Things like Julia, OCaml, Haskell, etc. I personally think Julia is pretty dope and is what I would use today for bioinformatics research. Or maybe if I was feeling a little subversive, K/Q or J. Q’s time-series database KDB+ could probably be used for sequences. And maybe even for great effect. And the performance would be off the charts!

It seems like the purpose of this is Python without the performance penalty, which doesn’t make much sense to me. I’ve found Haskell absolutely perfect for bioinformatics as most operations you are doing are functional data transformations. Moreover, it’s pretty damn fast if you need it to be.

I’ve been out of the field a long time though (Roche-454 was still the main workhorse at the time). But let me tell you, bioinformatics is/was a fucking shit-show. The tools and ecosystem are/were like Linux in the mid 90s: fucking terrible. And another language is just gonna make it worse.

Hi guys,

one of the authors (Ibrahim) here! Thanks a lot for the comments--- we definitely appreciate them!

A quick explanation why we built Seq:

- We were not happy with the existing bioinformatics libraries for various reasons. And honestly, while Julia is amazing project (and we do talk to the Julia team from time to time as they are located two floors above our office at MIT), it never 'clicked' with us or many other people in the field.

- While the main application domain is bioinformatics (that's where we came from), Seq is pretty much a strongly typed statically compiled Python. One of the main goals we had was to push the boundaries of how much stuff in Python can be deduced by compiler, as we loved the Python's syntax (Seq is to Python as Cystal is to Ruby--- or at least that is what we are aiming for).

- We do not want people to learn a new language--- Seq should be pretty much a drop-in replacement for Python, at least for most scientific/bioinformatics software. There still remains a small gap, but we are actively working to close it.

- At some level, Python libraries cannot cut it, especially when dealing with next-gen sequencing data. Also, owning the whole stack gives us the control to perform low-level pipeline optimizations. Chief example is out prefetch statement that is rather hard to implement in other languages.

Also, check out the paper (https://dl.acm.org/doi/10.1145/3360551) for more information. Let me know if you cannot access it for various reasons.

Fascinating, but there exist already htslib [0] bindings for Python (and many other languages). htslib truly is the standard library with respect to high-throughput sequencing data file access, and with high level bindings, we can already write something like:

``` for seq in bamfile: print(seq.pos) ``` or whatever.

[0] https://github.com/samtools/htslib

Response from the BioJulia developers: https://biojulia.net/post/seq-lang/. High-level points:

* they were able reproduce the performance comparison between Seq and BioJulia

* BioJulia spends most of its time on these benchmarks validating and transcoding data into a more compact, efficient representation of gene sequences

* Seq, on the other hand, operates on the raw ASCII input data and does no validation

* BioJulia devs implemented Julia types that representing gene sequences the same way as Seq does in less than 100 lines of Julia code [1]

* when using the same representation in as Seq, BioJulia was significantly faster than Seq

* BioJulia devs were able to further optimize transcoding of gene sequences to get a 10x performance improvement [2]

* with these improvements BioJulia reaches similar performance to Seq while still doing validation and using less memory

The full post is well worth reading. For me the main takeaway is that there's no real need for a domain-specific language at least not based on these results. Julia is already a great language for this kind of work and you get C-like speed and JIT compilation for free.

[1] https://github.com/jakobnissen/SeqLangBenchmarks/blob/master...

[2] https://github.com/BioJulia/BioSequences.jl/issues/86

Why not just contribute to BioJulia[1][2] instead? It is a more mainstream language but the one that fits computational tasks better than Python.

[1] https://biojulia.net/

[2] https://github.com/BioJulia

I'm especially interested in how sequence searching and matching work in libraries like this. Seq has a "match" statement for this task, which implements ACGT characters and _ for a single wildcard base and "..." for multiple wildcard bases, and a recursive matching system I haven't quite grokked yet.

Personally I'm more comfortable with a regular expression syntax, so would prefer "." and ".*". Actually, even better than "." is "N" from the IUPAC notation: https://en.wikipedia.org/wiki/Nucleic_acid_notation

The IUPAC notation is nice because it standardises "character classes" for working with nucleic acid sequences. For example, "B" is "[CGT]".

I wrote a module a while ago for searching nucleic acid sequences with regexps: https://metacpan.org/pod/Bio::Regexp

When working with sequences there are a bunch of things to think about that aren't really obvious from other types of data (or at least weren't obvious to me!)

Exhaustive search: Often a regexp can match in many ways, and most regexp systems don't provide a way to get a complete list of them all. Fortunately there is the Regexp::Exhaustive perl module which is what I used. The way this module works is pretty awesome. It adds a special "FAIL" directive to the end of a pattern, so that the match can be recorded and artificially failed, triggering the regexp engine's backtracking mechanism to back up and find the next match (if any).

Reverse complements: Because DNA is double-stranded (well, usually... this is biology after all) there is a "complementary" pattern on the first strand that corresponds to the pattern you are interested in on the other strand. You almost always need to search for both. And what's more, DNA is directional so you actually need to search in the reverse direction for this complementary pattern. You can either reverse and complement your sequence (which seq has a special ~ operator for, neat!) and search again, or reverse and complement the search pattern itself, assemble a single combined regexp (Regexp::Assemble module), and do a single scan over the data, which is what Bio::Regexp does.

Circular DNA: Some DNA (plasmids) are actually circular in shape, meaning the start is connected to the end. So a comprehensive search needs to check for cases where the desired patterns span the arbitrary location selected as the "start" in your sequence.

Like a few others have mentioned, I do wonder how big the marginal return in baking in these features is.

From a quick look, stuff like FASTA/BAM parsing, translation, etc can be implemented in C-land a la numpy, and called from Python, right?

A language like Swift would also support the addition of powerful user-defined operators, as in the case of Swift for TensorFlow[1].

Language adoption is hard to drive, and I wonder if having domain-specific library calls built-in is worth the added effort for people in the field.

[1] https://www.tensorflow.org/swift

My personal impression isn't that bioinformatics needs a full language, but more tools in popular environments to lower the entry barrier for good software engineering practices.

Every bioinformatics codebase I've looked at has been a downright mess. Basically an ad hoc collection of scripts that transform data this way, maybe rendering some graphs or such, relying on 100 unstated assumptions. Nothing is maintainable, and often rely on messy approaches like loading your entire data set into memory (works fine for your 1-10GB data set, then not so much on a larger one) or what I would gently describe as mainframe compute in place of real software engineering.

I made a short explainer video on the language here:

https://youtu.be/5bk4Wc5Op2M

This should be marketed as a language for genomics not bioinformatics. There is more to bioinformatics than just genomics, but this language (at least the beginning of the docs that I looked at) seems to be marketing exclusively for genomics analysis.

Not a knock on the project at all, but it doesn't seem like someone doing analysis on cell images is going to get much out of this language.

As a working biological data scientist... This is useless until a tool I need, which I can't get anywhere else, is written in it.

My suggestion would be to rewrite some of the most popular tools, like bwa, in this language, and show the comparative performance etc.

Then, write a comprehensive open source package with great documentation and maintenance in this language, to demonstrate to others your investment.

Then, maybe, it will get some traction. But honestly, C, R and Java are so embedded it will be a hard road.

Since in bioinformatics, some processes are vastly more time-consuming than others it is not clear what the benefit of designing a new language is as opposed to adding an optimized, C library into Python itself.

What functionality would not work the same way if this were a python library?

That being said, it is a neat and impressive effort, I do fear though that it will have a whole lot less uptake as now requires bioinformaticians to learn a third language (Python, R and now seq)

Historically I note a similarity to the 'Mothur' and 'Qiime' split.

'Mothur' is a dynamic language similar to 'seq' whereas 'Qiime' is Python glue across many libraries. Frankly, I like 'mothur' better, but 'Qiime' is whole lot more popular.

I'm a bit surprised at all the negative comments here. I hope it isn't too discouraging for your team, because as author of a 50K LOC Python app (HashBackup), I could really use this! I love the Python language but sometimes the performance is a drag. For example, to plan a restore when the data isn't local, HashBackup has to traverse every block in every file to be restored and figure out when to load the block and when it can be released from the cache. This isn't particularly difficult, but for very large restores it requires long loops using large lists, arrays, and/or dicts. Parts are coded in Cython, and that works well for easily-isolated functions, but not so great for something like the restore plan that needs database access and is referenced during the restore.

I ran a small 10M entry {int:int} dict benchmark. In Python 2.7, the test used 1.1GB of RAM and about 8 seconds. In D (fully compiled) the same test used 881MB and 7.4 seconds. Here's the D version:

  import std.stdio : writeln;

  int main() {
    int[int] map;

    foreach (i; 0..10_000_000)
      map[i] = i;

    foreach (i; 0..3) {
       foreach (j; 0..10_000_000) {
         map[j] = map[j] + 1;
       }
    }
    return 0;
  }

In Seq it ran in 5 seconds and used 395MB. Here's the test program and Seq run:

  map = dict[int,int]()
  for i in range(0,10000000):
      map[i] = i
  for i in range(3):
    for j in range(10000000):
      map[j] = map[j] + 1  
  print map[12345]

  [root@hbseq ~]# /usr/bin/time -v ./map
  12348
    User time (seconds): 4.67
    System time (seconds): 0.55
    Percent of CPU this job got: 97%
    Elapsed (wall clock) time (h:mm:ss or m:ss): 0:05.36
    Maximum resident set size (kbytes): 395532

Looks pretty great to me, especially if I don't have to do a major rewrite! I guess I could have hit a case where Seq happened to have a higher HT load for 10M entries and D just did a resize, so it would be good to run the same kind of test at a lot of different hash table sizes. But the Python results are pretty terrible space-wise.

I don't really understand why to create a new language. But I like the improvements they made for the Python lang

As someone who noodles around in python, and knows nothing about bioinformatics, this looks very interesting for the 'nice' things it does to basic python like forcing single type of returns, array controls. Its almost like a 'safer-python' set of constructs that would be of great use to the general python programmer. And the pipe operator is very cool.....gotta try that out soon.

Hmm, it's interesting that you have DNA base sequences support built in. But (you know someone was going to ask ;-) ) I don't see similar support for Amino Acid sequences, or encoding/decoding between the two. Is this a deliberate design choice?

It is like a Julia-Nim hybrid.

I'm curious how this differentiates from Nim other than the builtin types for bfx stuff. Still cool either way and great to see something else joining this space.

How does it compare to Futhark/J? I'm especially curious how it compares to calling Futhark from Python. Are there special data structures?

Heterogeneous collections and no inheritance/polymorphism seems like a bad combination.

Is it possible to write a python module with seq?

This is similar to Crystal?

Obligatory “this name is already in use in this field” post: Seq is the name of a very popular structured logging sink that provides a web app with a query language interface for searching through and graphing log streams, often paired with Serilog when used in the .NET world.