Just one technical comment - that parser will break on many FASTA files. They can also include '.' & * characters ('.' is often output by sequence aligners to indicate sections outside of the model, * is there for aa but not dna/rna and I have seen them in sequencer output for uncertain positions). I think this kind of edge case variance is possibly why the various Bio-* tools have not been very successful in all, you end up spending hours trying to work out why it's not working with your particular inputs.
On a personal viewpoint: "Take a look at the REPL dump below and see if you can work out where we are implicitly using $seq as a Str rather than Seq:" - this section pretty much sums up why I've always hated working with Other Peoples' Perl. To accurately follow someone else's perl requires either an encyclopedic knowledge of the language and it's edge cases or very heavy documentation and a set standard of coding practices. It's a shame that this looks true for Perl 6 as well.
Great for personal use, I guess. But by and large I'm glad to see more projects switch over to python.
Sure but the nice thing is you can just add an extra catch all to the Grammar. Or better yet create a new Grammar for that specific sequencing machine and override the token.
For example just adding:
rule sequence:sym<seq> { <-[\>]>+ }
Would mean you can just grab anything without worrying about if it's correct. Because a tonne of people abuse FASTA doesn't mean you shouldn't at first attempt to enforce something valid. The nice part of that is you know within the parser if you have clean data or something with a badly defined character like * in DNA when it should just be X or N or one of the other confusion characters.
Also it's an advent post not some formal public codebase ;P The focus was on showing some of the more obscure features of Perl 6 that aren't like other languages. Python can look as crazy as you like, it's just the majority of people using it fall into the camp of avoiding the crazy. For example here's (https://gist.github.com/MattOates/6195743) some crazy Python I've written. The six line list comprehension is especially readable :P
There are a couple of problems with starting with a strict approach. I predict the performance is going to be worse than something which allows (almost) all characters, because of validation overhead.
More seriously, I wrote a full grammar based parser engine for Biopython along these lines. The Achilles heel was error recovery. Suppose record #123456 in a database contains a syntax error (I can list dozens I found, like a copyright date of 19995 breaking my four-digit definition, or SWISS-PROT record with the copyright lines in the wrong spot.) How do I have the parser ignore that record, report the problem, and continue parsing with the next record?
Error recovery in context free grammars is hard, but doable.
Then how do I say that "if the ID is XYZ and the version number is 12 then use this work-around because that release had a bug which has since been fixed"? In practice and unlike a computer language, the format spec ends up being less important than the actual data, so you end up with lots of tweaks to handle vendor errors.
P.S.: Wouldn't
names = [SQLiteMagician.types.get( type( d ).__name__) for d in data[0]]
types = ["NONE" if t == "NoneType" else t for t in names]
have been less crazy, and to be preferred? Also, and I'm not sure what that code is trying to do, but sys._getframe() and linecache from the standard library should help, as might this example, which I think does what you want without using the ast module:
import sys
def show_me():
caller = sys._getframe(1)
f_code = caller.f_code
names = f_code.co_varnames[:f_code.co_argcount]
local_vars = caller.f_locals
for i, name in enumerate(names):
print repr(name), "=", repr(local_vars[name])
def spam(x, y, z=4):
a = "%s %s %s" % (x, y, z)
show_me()
spam("Hello", 2)
Great thanks! So I think there are two things about your code which are different. You are getting the names and values of the spam function with show_me() ? My version gets the next level up which is the names of the variables passed in to spam. I believe I tried to use some of the same standard functions at first. But this had problems when a function is called over more than one line of source text, and some other limitations. You don't necessarily get back what you expect without properly going into the source and parsing. I'll take another look though since I might have just missed something obvious. My Python isn't /that/ strong.
Completely agree with your list comprehension suggestion. The list comprehensions are crazy by design to troll monoglot Python people who hold a false belief about readability of code just being identical to language syntax. If I handed you that code without the comments, and without sensible variable names. Which is frequently what I do get handed to me by people I work with. I've found a myth persists about Python being inherently readable (the opposite of Perl!), and many scientists who don't have a programming background just go with this idea especially as it gets them off the hook of learning to program well.
With parsing and grammars I am still more happy having a parse fail and tell me what the garbage data was and where. It's very rare any code would be able to recover from completely arbitrary and unknown problems with a file. With exceptions being thrown you might at least start processing the bit of the file you already have and just set a state that shows the whole thing didn't run. Asking a language to do more for you feels a bit much. Your point about validation overhead I'm not convinced by since you have to check the characters for a record start/end state anyway. Plus it just means garbage data gets even deeper into your program! At some point you have to check it's valid. Fail fast. Checking a small set of match states isn't going to kill performance if the grammar/regex engine is any good. One of the ways I'm getting around this in my Perl 6 libraries is just having a thread that deals with file IO, parsing and object creation. Then a users main program can just deal with the stream of sequence objects as they come in.
Getting the variable names used for the caller is of course much harder, and meaningless for a call like f(1+g(b)). py.test, which made a more serious attempt at handling the f() case, gave up on that approach, preferring to rewrite the original code behind the scenes rather than infer it backwards.
So you're a troll. Thank you for letting me know. By trolling the uninformed you seemingly have no meaningful point other than to point your finger and laugh at the views of others. This action runs counter to your stated goal of wanting scientists to be better programmers. What good do you think it will do to the world to troll people using your novice Python skills?
You do realize that when someone searches for "Matt Oates", "bioinformatics" and "Python" they are likely to find this page in the future, and learn about your love of trolling, right?
I do not want to live in a troll world, just like I don't want to live in a world where people use false and undercutting statements disguised as joviality. Instead, I teach scientists how to be better at the programming they need to do good science.
Your last paragraph tells me that you have very little experience in writing parsers. In a record oriented format, and with hand-written parsers, resynchronization to the record boundary is easy. In Perl's it's almost trivial, because nearly all record-oriented formats have a boundary which can be handled with $/ ; FASTA is an exception but still easy to handle. On the other hand, my point is that it's much harder to resynchronize with a context-free grammar.
Here's an example of a failing case with my strict parser. HAS2_CHICK was the only record in SWISS-PROT 39 with a DT line of:
DT 30-MAY-2000 (REL. 39, Created)
This failed because the spec said that the "REL." was only supposed to be a "Rel." My parser failed. Every other parser basically said (by lack of full validation) that they weren't interested in that field, and ignored it.
You can certainly condemn the entire SWISS-PROT 39 release as "garbage" for having a single record which didn't match the spec. Effectively everyone who wants the data will disagree with you. You can also train your parser to ignore those details, and retrain it for every release, though you'll find it increasingly difficult to accept the format as it changes across 5 different database releases.
Here's another example. Some of the GENBANK records included a qualifier of "/primer", which wasn't in the specification. My parsers failed, because it validated the qualifier names. And you know what? No one really cares about that failure, so I had to tweak my parser to allow the non-standard name. The same with EMBL: the ID line is supposed to have the molecule type (one of "DNA|RNA|circular DNA|circular RNA"). In actuality? Record S40706 had a type of 'XXX'... and so again, my validating parser stopped, when in reality almost no one cares.
As these are further examples of simple differences causing a failure within a record, your bringing up "completely arbitrary and unknown problems with a file" suggests to me that you have little real desire to engage in my point, and rather prefer to troll me by diverting the topic. In case I am wrong, please assume that Perl's $/ or a similarly simple method is sufficient to determine the record boundary and resynchronize.
The underlying issue is that a database format spec changes after almost every release, so a validating parser ends up having to handle a family of formats, not one. Which is much easier for a line-oriented parser to do, since those formats are written with line-oriented hand-written parsers in mind, and usually updated in such a way as to not break those parsers.
When you find out that your nice, elegant, validating context-free parser breaks after nearly every release, due to details that are irrelevant to what you care about, then you might conclude that the rest of the world is a bunch of ignorant fools who can't program; or reject the idea of a validating parser for this sort of data and switch to a line-oriented one; or provide a recovery mechanism to determine what to do with the records which couldn't be parsed, and continue if requested. I do the latter now.
(It's also difficult to use a CFG-based parser that doesn't know about record boundaries to parse files that don't fit into memory even though each record does. I don't know if that's still an issue in bioinformatics these days, but it certainly is in cheminformatics. In any case, a mix of simple record extraction and a CFG-based grammar suffice, excepting for records containing an entire genome, which is outside the scope of this sort of task, and for PIR, which had header and footer components to the file.)
Regarding performance, don't trust your beliefs or mine. Measure it. Empiricism trumps faith. Try benchmarking a hand-written, line-oriented FASTA parser against your grammar. I wrote one just now in Python. The regex which used "[^>\n].*\n" was about 2% faster than checking the [specific sequence letters]. A 2% here and a 2% there and those validations quickly incur noticeable overhead. By comparison, my line-oriented parser took about 1/2 the time of the regex-based parser.
Of course, Perl has a different runtime environment and will have different conclusions.
Then again, there's already a history in Perl of variable stringency levels for improved performance. Swissknife (doi: 10.1093/bioinformatics/15.9.771) was written as a lazy parser because Perl5 wasn't fast enough for full validation and parsing of each record. Partial evaluation gave a 4x performance improvement when full validation wasn't needed.
> my point is that it's much harder to resynchronize with a context-free grammar.
So, don't use a context free grammar. [1]
Use something else, eg a PEG. [2]
Not only are Perl 6 grammars PEGs, not CFGs, but individual rules are closures (of arbitrary code) as well as input matchers so one can write arbitrary code to handle anything that happens during a parse. [3]
As a very trivial example, in a grammar I wrote for parsing a CPAN control file, I added a rule for handling unexpected field types. It prints a warning to stderr and then keeps on going. [4]
Imo Perl 6 is, among its many strengths, an innovative parsing tool with at least three big problems right now:
* The grammar engine is slow. [5]
* Cat is likely a post Perl 6.0.0 feature. (Among other things, Cat is for enabling practical parsing of arbitrarily long and arbitrarily buffered data streams without having to tweak userland code.) [6]
* Macro, slang and parse tweaking features are still evolving. [7]
"The Cat type allows you to have an infinitely extensible string. ... Then a Regex can be used against it as if it were an ordinary string ... The Cat interface allows the regex to match element boundaries with the <,> assertion ... Strings, arrays, lists, sequences, captures, and tree nodes can all be pattern matched by regexes or by signatures more or less interchangeably."
To add some more background for grammar-based bioinformatics parsers, take a look at SRS (Sequence Retrieval System) from the early 1990s, by Etzold and Argos. It was the last generation of purpose-built database systems, before off-the-shelf object-relational systems took over, and used a set of grammars in a declarative language to parse and index flat files.
Sorry, I was just pointing it out in case for anyone copying the code as they are both very common cases - especially dots as gap characters.
The Grammar system does look quite nice. Parsing is definitely a Perl forte - though I seriously hope bioinformatics in general will continue to abandon these crazy formats that can't be read easily by humans or machine.
I'm increasingly of the opinion, after having to deal with biologist data format hackery for many years (they have no shame!), that every biologist should be taught relational database normalisation. That would make them actually stop and think about the structure of the data they're recording and how to make it computationally tractable. Even more, they should not even start their experiment until they've cleared their data recording structure with a bioinformatician.
Then we ban them from using flat files. Related tables or nothing. I've wasted too much time, and they're not competent enough to work with JSON/XML/Excel.
The FASTA parsers for the various Bio* all support additional characters; if not then that would definitely be a bug.
Saying that, 'support' completely depends on what the FASTA is and where it comes from (assembly, sequence database, alignment tool, etc), which is something the parser/grammar can't define up front from such a simple format.
Much of the problem comes when validating the sequence for a specific alphabet or symbol set. If the alphabet isn't explicitly defined (again impossible to determine from the format without guessing) then you can certainly run into problems.
This is also an issue when using FASTA for both regular sequence data and for alignments (e.g. length of the sequence would have to take into account possible gap characters generated from various tools like '-?.', or stops in protein seqs like '*').
These are tools that have been in widespread use for ~15 years, so if you have run into problems you should be more explicit in what they are so they can be addressed.
BioPerl has probably got a lot better over the years, I haven't touched it in a long time (5 years+), so my views may be very out of date. I don't really want to get into reviewing Bio* libraries generally :)
There is a general problem I think though, and other people have experienced the same. I've been trying to reason why, but I think it's due to multiple factors. By the time you've managed to get that bit working reliably with BioPerl you could have written it yourself. And it would be faster and use less memory. BioPerl has the problem of trying to solve every case, whereas normally I'm working with a limited subset of possibilities.
So why? I think it's a combination of dodgy, manually hacked inputs (e.g PDB files!), the learning curve, poor backwards compatibility preventing upgrading, and, ah I don't know. Every time it's seemed like a good idea, and yet every time I've ended up abandoning BioPerl. Maybe it's too integrated with itself as a library of functions?
Take parsing a FASTA. By the time I've read the BioPerl documentation, I could have already written that one line split string statement (because I e.g. know in advance the sequence string is always on one line). It's hard to overcome that laziness and make the commitment to learn it and become fluent in it.
>>this section pretty much sums up why I've always hated working with Other Peoples' Perl.
I think that statement was really more of 'The remainder is left for the users as an assignment' kind of statement.
>>To accurately follow someone else's perl requires either an encyclopedic knowledge of the language and it's edge cases or very heavy documentation and a set standard of coding practices.
To accurately follow anything non-trivial written by anybody in any programming language will require you a good grip over the language and its corner cases. There is no magic here, if a language doesn't offer facilities to heavy lift a certain thing you will end up writing those yourself.
>> To accurately follow anything non-trivial written by anybody in any programming language will require you a good grip over the language and its corner cases. There is no magic here, if a language doesn't offer facilities to heavy lift a certain thing you will end up writing those yourself.
I'm going to have speak subjectively and relative to my own skills here. Without rigorous use of standards Perl 5 is a nightmare to read. With very clear and consistent style then it can be readable, but this has to be enforced extremely strictly for any code base beyond trivial. Perl's 100 ways to do things is not a positive in this context.
Case in point, I was given a bit of perl two weeks ago by a bioinformatics lecturer. Did not use strict or warnings; variables weren't declared in scope; no subroutines, let alone Moose objects. And this is typical of ~90% of bioinformatics perl that gets written. What's more, it was a trivial script that was simply parsing some BLAST matches from CSV, checking for overlaps and assigning a type based on the match gene combination - and was completely indecipherable. It took me two days of refactoring to pick out all the details. And he managed to use a couple of magic symbols I'd never seen before.
The first time I was given a bit of python, without learning the language at all, I was able to extend the script and get it do some new things (it was considerably more complex than the aforementioned perl script).
And even when it is well written and documented - I did work with some very good perl programmers - it's tended to end up as very fragile and difficult to extend or even substantially redevelop due to the heavy use of context for determining behaviour (and the actual context can be difficult to determine). Unit testing had to be very extensive since it was difficult to isolate changes.
I'm going to have to be precise here - I'm attacking Perl 5. It just seemed from that section that Perl 6 doesn't solve the major problem (difficult to deduce the context and work out the resulting effect) I had with Perl 5. I may be completely missing the mark and actually readability is much better. I really hope so.
>> Perl 6 isn't production ready yet. But even a cursory glance over the features tells, Python isn't even in the same league of tools Perl 6 will be.
Vapourware is vapourware. Perl 6 was meant to be production ready a decade ago. On the other hand, I'm interested in knowing what you think Perl 6 brings that will give it the edge over Python? Do you see it displacing Python as rapidly as Python replaced Perl 5 as the science "glue" language of choice? With the apparent rise of Julia as well, could it just be too late?
I'm not talking about Perl 5 per se. But I've worked with C, Java, Pig, Python, Perl and a bit of shell scripting so far. And I've never had a situation where I could simply glance over something and know all of it(This is for reasonably big applications). I'm currently working on a big Java application where logic is often deep buried below Factor classes, Interfaces, Enums and layers of classes simply passing control to level below. And when you find the code that does something finally, you will see it wrapped will kinds of exception handling code. Its simply too verbose, with pointless ceremonial code.
In pig, any reasonably large workflow has had me hooked for hours working with stub data trying to figure out how data flows and comes out.
In fact I see this nothing new to Perl. If you are working with any serious application. A good touch over the language plus experience with handling things with patience is a must. Python is no exception to this rule, neither is any language I know of.
>>I may be completely missing the mark and actually readability is much better. I really hope so.
It depends on how you define readability. I would say verbose code is equally unreadable especially if you work with large applications. If you simply don't like regular expression you have a different problem. In some cases use of regular expression is simply unavoidable. IMO Python works very well for applications that interact with standard data sources(Databases, API's, XML's). Once you get the data to some very easily ingestable form, and there is little work to do from there on, Python works just fine. But Perl works fine for other use cases, heavy lifting and other kinds of data work. Perl may not be the right language for all use cases, but neither is Python or Java.
>>Vapourware is vapourware. Perl 6 was meant to be production ready a decade ago.
No, I don't think they ever gave out a date when they wanted a finished version out. Besides Rakudo is feature ready already. They are just working out last bit of performance issues and good to have things. And plan to have it out by 2015: https://fosdem.org/2015/schedule/event/get_ready_to_party/
>>On the other hand, I'm interested in knowing what you think Perl 6 brings that will give it the edge over Python?
Grammars, Macros, Concurrency, Traits, Types, Improved regexes, Better OO, Functional programming, User defined operators, Mutable grammar, Junctions, Currying, Continuations, Multiple VM backends to name a few.
>>Do you see it displacing Python as rapidly as Python replaced Perl 5 as the science "glue" language of choice?
Perl and Python both co exist and will continue to, albeit the web frameworks market is dominated by Python and Ruby at this point in time, Perl is still pretty big in a lot of places. This isn't Perl or Python, people use tools which work for them within their project limitations.
But I guess people will use what they find useful. And that only individual projects can decide as to what works best for them.
>>With the apparent rise of Julia as well, could it just be too late?
Reading the Perl 6 spec will tell how different Perl 6 is with regards to everything else. Even with the delay Perl 6 has suffered, I don't see any other tool that has the feature set Perl 6 has to offer.
Absolutely, with any serious application, it takes time to break it down and understand it. But the majority of perl generated in bioinformatics are basically advanced shell scripts and should be immediately comprehensible. This is rarely true, mostly because there's so many ways to do things that everyone uses their own personal conventions. Learning perl 5 is a bit like collecting pokemon.
As for the Perl 6 features, I had already looked at the list and was hoping you may have some insight as to why these might be killer features for bioinformatics in particular. e.g I know what currying is, but I don't think most computational biologists will care. I actually can't see anything there that says, "Science, do it here".
Julia is a very promising language due to its focus on mathematics, readability & performance and seems to be gaining traction. R has an unbeatable set of statistical libraries. Python is now deeply embedded and has IPython, and links back end to web very nicely. C/C++ & Fortran are are miles out in front for number crunching. Java is excellent for reusable code and distributed development, and can be used in almost any layer and has proven very difficult to dislodge as the general purpose language. At the moment I don't see a niche being available to Perl 6, unless something really useful for scientists (like IPython notebooks) is brought to the table. If they really backed it as a parsing language, and provided some really sweet tools for handling biological data files - e.g. something more like an interactive IDE rather than having to write out a script in emacs. Actually give me that, I'd be very happy. But I honestly think it's the associated tools & ecosystem which will determine if Perl 6 can succeed, not a laundry list of features.
And will it ever come out? :) I remember going to a Perl learning course in 2002 where the instructor was very excited about the new version 6, which was going to be out by 2004.
Ermmm I wrote a whole blogpost on some of the killer features for Perl 6 in bioinformatics :P There are a lot of primitive types that are useful for scientific programming. A lot of my bioinformatics work at least gets reduces to bag models and set operations on big files quite often. On the computational side the stuff with native arrays is being built right now, so it remains to be seen how well that will pan out. But in theory you will get an array of native types that will be about as efficient as C. So that's like numpy included as a core part of the language. I left out the nice concurrency/parallel model and native code integration too which are really killer features in general. The ability to write first class domain specific languages that integrate fully into the compiler is kind of incredible. An example is the recent SQL Slang someone hacked together https://github.com/tony-o/perl6-slang-sql/ . You can just inline write a SQL statement and it acts like a for loop over the results returned. Kind of cool!
Thanks for the interesting reply. At the moment I guess it's a case of wait and see. I think it's fair enough to have some non-specific skepticism about Perl 6 given the history, however it's clear the designers are putting a lot of thought and love into it. And I can see how it will provide a powerful means of generating pipeline-type analysis tools.
However, it actually makes it seem to me that even more so Perl 6 makes the same "mistake" as Perl 5. The surface area of the language now seems fractal in complexity, and unless rigorous discipline & best practices are applied it's going to be very hard to build working complex projects. So great for individuals, very powerful for highly disciplined teams, but not good for the below average programmer - as 90% of biologists are. Perl 5 was full of magic, but in the hands of most biologists & less disciplined bioinformaticians that was a bad thing.
Take that SQL Slang example, great for quick scripts, but it's seems a long long way from a production ready library that can be used in a complex application. And maybe I just don't see the great advantage in being able to hide a couple of lines that make it exactly clear about what you're doing. If I'm maintaining code I want to see those explicit statements if possible. And really it's not exactly a big win over "&sql($statement)".
Anyway, it's interesting and luckily I can observe from a distance now.
I don't think they ever gave out a date when they wanted a finished version out.
"Larry Wall and other active Perl porters and Perl helpers met on Tuesday afternoon at Perl Conference 4.0 and mapped out a what is planned to become a complete rewrite of Perl that will become Perl 6 in 18 to 24 months."
-- Linux Today, byline July 19, 2000http://www.linuxtoday.com/developer/2000071901704OSSW
Perl is nothing but corner cases and magic context sensitive punctuation. If you think Python is as hard to understand as Perl, then you've never used Python.
Bioinformaticians are generally capable of writing truly awful code... that's why I like to call myself a Computational Biologist. I still write crap code, but no one's expecting any better from a biologist ...
On a personal viewpoint: "Take a look at the REPL dump below and see if you can work out where we are implicitly using $seq as a Str rather than Seq:" - this section pretty much sums up why I've always hated working with Other Peoples' Perl. To accurately follow someone else's perl requires either an encyclopedic knowledge of the language and it's edge cases or very heavy documentation and a set standard of coding practices. It's a shame that this looks true for Perl 6 as well.
Great for personal use, I guess. But by and large I'm glad to see more projects switch over to python.