Bioinformatics and the joy of Perl6
perl6advent.wordpress.com
perl6advent.wordpress.com
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.
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
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]
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)
'x' = 'Hello'
'y' = 2
'z' = 4Completely 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.
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)
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.
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]
[1] http://en.wikipedia.org/wiki/Context-free_grammar
[2] http://en.wikipedia.org/wiki/Parsing_expression_grammar
[3] http://en.wikipedia.org/wiki/Perl_6_rules
[4] https://gist.github.com/raiph/63508a5f4ade000acc7f
See the `data-item:dst_GENERIC-MATCH` rule; note the embedded `{ warn ... }` closure.
[5] http://www.reddit.com/r/perl6/comments/2ec2hd/rakudo_perl_6_...
[6] http://perlcabal.org/syn/S03.html which includes:
"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."
[7] http://perl6advent.wordpress.com/2014/12/17/day-17-three-str...
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.
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.
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.
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.
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.
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?
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/
http://perl6.org/compilers/features
>>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.
The full list is at : http://perlcabal.org/syn/
>>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.
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.
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.
"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, 2000 http://www.linuxtoday.com/developer/2000071901704OSSW
Cross-checked against an eyewitness report by Chris Nandor, published at http://use.perl.org/use.perl.org/articled5d3.html?sid=00/07/... on July 19, 2000.
Just getting into bioinformatics and focusing on python. Glad to hear some validation about that decision.
It's interesting that the example supports a comment in the FASTA record. I wonder why the author included it. I've never seen a comment used in the wild, but then again I haven't worked with FASTA files for some time. My experience was that nearly every tool would break if you give it a FASTA file with a comment. For example, the BioPerl FASTA parser is at https://github.com/bioperl/bioperl-live/blob/master/Bio/SeqI... , where you can see it doesn't support a comment.
Since the chance of breakage is so high, my belief 10 years ago was that this original feature from Pearson is dead, and will never be revived.
The Bio* refer to the additional stuff after the ID as the sequence 'description'. There is a subtle difference between that and a simple comment, though I have seen even NCBI use this as the 'junk drawer' for any additional seq information (such as the headers in the nr database).
Re: FASTA comment support, I haven't ever seen it's use in the wild.
>ENST00000517143 ncrna:snRNA
ATGC...
The problems comes to an ambiguity in the Pearson's original FASTA distribution, from http://faculty.virginia.edu/wrpearson/fasta/fasta3/ . In my copy of FASTA (fasta-35.1.5) in fasta20.doc is the following:
0 Pearson/FASTA (>SEQID - comment/sequence)
...
Standard library files. These are the same as plain sequence
files, each sequence is preceded by a comment line with a '>'
in the first column.
...
I have included several sample test files, *.AA. The first
line may begin with a '>' or ';' followed by a comment. The
text after ';' in other lines will be ignored. Spaces and
tabs (and anything else that is not an amino-acid code) are
ignored.
...
This is often referred to as "FASTA" or "Pearson" format. You
can build your own library by concatenating several sequence
files. Just be sure that each sequence is preceded by a line
beginning with a '>' with a sequence name.
>H10_HUMAN | 90538 | HISTONE H1' (H1.0) (H1(0)).
TENSTSAPAAKPKRAKASKKSTDHPKYSDMIVAAIQAEKNRAGSSRQSIQKYIKSHYKVGENADSQIKLSIKRLV
TTGVLKQTKGVGASGSFRLAKSDEPKKSVAFKKTKKEIKKVATPKKASKPKKAASKAPTKKPKATPVKKAKKKLA
ATPKKAKKPKTVKAKPVKASKPKKAKPVKPKAKSSAKRAGKKK
(Also, there are '.seq' records with DNA in them, so the comment about ignoring non-amino-acid codes is only referring to protein FASTA files.)
Obviously the text after '>' is important, while the ignorable text after a ';' is not ... unless it's the first line of the file. If so, what name should be used to distinguish between one and another? The code calls the first line a title, for example.
Some people implemented ';' as a generic comment field, to be ignored. (See that Talk page for a couple of examples; "read.fasta in the seqinR package and by the function readFASTA in the Biostrings package") Most others did not.
After Pearson came NCBI. They describe a backwards compatible subset of the original FASTA at http://blast.ncbi.nlm.nih.gov/blastcgihelp.shtml . It calls the '>' line a "description line (defline)". The NCBI toolkit further break up the record into "sequence, description, and identifiers" (see http://www.ncbi.nlm.nih.gov/books/NBK21097/ ).
BioPerl and Biopython, and I assume the other Bio* languages, follow NCBI's lead and use the same, or similar, names.
I am a member of the NCBI FASTA camp, not the Pearson FASTA camp, so when I see the term "comment" I think it unambiguously refers to text after a ';' line in a Peason file. I can see how someone from a different intellectual heritage would call the description a "comment", but as most of the world uses NCBI FASTA and not Pearson FASTA, I think it's a bit confusing to do so.
I've been in the bioinformatics space for about 13 years now and I've only come a small handful of perl scripts, all of which have been legacy programs. Even when I started in the field, perl was viewed as somewhat anachronistic. All of the end-user software I've written (read: For researchers) has been in R, Java and a wee bit of scala.
These days all of the researchers I know are using some combo of R, Python and to a lesser extent MATLAB. Many have useful C and/or Java (and/or Scala) skills as well.
My personal experience is obviously not exhaustive but I find it hard to believe that many people are using Perl these days, and IMO that's a matter of good riddance.
Perl is mostly used in the niches it's traditionally good at, munging data or cleaning it up with quick to write one-off scripts. The data is then fed into a higher performance C or Java program for the actual analysis.
They tried to use Python for that kind of task, but the language is stricter than Perl (you end up writing lots of exception handling code) and the text gobbling faculties are just ever so slightly not as good and it ends up just being easier and quicker to Perl the problem away. However, the newer folks are coming in not knowing Perl at all and it's likely it will simply go away in a few years and everything will be done in Python anyways.
Python and Java tend to be used on the back-end of the web reporting systems they build since there's better server-side support and tooling for that case. And I think some of the analytic tools are cross platform Java desktop clients.
As someone who knows the author of the post, that is mostly true ;)
In all seriousness, while I generally agree with your sentiments I don't think that Perl is dead within Bioinformatics. Not least since major data providers, e.g. Ensembl, write their tools in Perl.
Oh I know it still exists, but while this is admittedly anecdotal on my part (and likely heavily biased by the institutions I've been a part of) I haven't seen a single new person who defaults to Perl in a very, very, very long time.
And to be clear, I'm not trying to rip the author or the article, it just seems ... oddly timed, like writing a contemporary article on how to use COBOL in the business computing space.
I'm not sure I agree with that, however.
Perl 6 is probably best described as a dialect. The look and design of the language are very Perl-ish, but some of the key sore points from Perl 5 (OO, concurrency, etc) are addressed and it has added a number of killer features. This is not Python 2->32; it is a major overhaul of the language, with no backwards compat beyond a suggested perl 5-compatible layer (I believe this is called 'v5').
The article touches upon a few (Grammars for instance), but I personally think the concurrency work will also be a real draw.
The other key difference is that Perl 6 is actually a specification with an official test suite and Grammar (STD). I believe the specs indicate that anything that passes the test suite can be deemed to support 'Perl 6', which really opens up the use of various backends. The Rakudo Perl 6 implementation has support for three (MoarVM, Parrot, and JVM).
The neat things is that you can use any language you want.
I haven't used perl myself, so can anyone comment on the merits of perl vs awk here? They often compete for this use case.