Hi,
I'm interested in seeing your favorite bit of code that you programmed yourself for bioinformatics purposes.
It can be a function, class, parser, you name it! Mention the language you coded it in.
Cheers to you bioinformaticians!
This almost sounds as bad as a bash script labeled "Unix Russian Roulette" I came across some time ago. The script was basically a short one liner that gave you a one in six chance (if you can call it that) of erasing your disk. I was linux-naive enough to try it and, luckylly, I didn't put "sudo" in front of it to gain admin rights... [sic]
My favorite codes are those that allow me to use powerful algorithms directly from my favorite language without any penalty. Here is such a code using the bx-python intersecter (a module written by Brent Pedersen). It queries 100,000 times a prebuilt data structure containing 100,000 reference intervals. Each interval is of a random size of up to 10,000 spanning a region of 10 million. And it only takes ... drumroll ... 7 seconds, thus performing at around 50 ms per query!
import time, random
from itertools import imap
from bx.intervals.intersection import Intersecter, Interval
def random_interval(x):
start = random.randint(1, 10**7)
end = start + random.randint(1, 10**4)
return start, end
# create 100K reference intervals
source_intervals = imap(random_interval, xrange(10**5))
print "Building:",
tick = time.time()
inter = Intersecter()
for i, ( start, end ) in enumerate( source_intervals ):
inter.add_interval( Interval( start, end, i ) )
print time.time( ) - tick
# create 100K query intervals
query_intervals = imap(random_interval, xrange(10**5))
print "Using:",
tick = time.time()
for start, end in query_intervals:
inter.find( start, end )
print time.time( ) - tick
The output is:
---------- Python 2.5 ----------
Building: 2.11999988556
Using: 5.25200009346
Output completed (7 sec consumed) - Normal Termination
chomp (@a = <cat list.txt>);
foreach $a(@a)
{
<whatever> system commands*;
}
Most of my project starts with a list of IDs, first step will be mostly a BLAST or hmmpfam search or to run a custom-program using files in list.txt
I think clear help function is really important for users. So, even if it is not a bioinformatic code, this is a useful code for me. It is a small code which reproduce the "man" aspect.
import sys
print """
test(1) User Commands test(1)
\033[1mNAME\033[0m
\t%s - small description
\033[1mSYNOPSIS\033[0m
\t\033[1mpython %s \033[0m[\033[4mOPTION\033[0m]... [\033[4mFILE\033[0m]...
\033[1mDESCRIPTION\033[0m
\tLong description
\twith differents paragraphs
\t\033[1m-f, --file\033[0m
\t\tThe file...
\t\033[1m-h, --help\033[0m
\t\tDisplay the man of this program
\033[1mAUTHOR\033[0m
\tWritten by Nico.
""" % (sys.argv[0],sys.argv[0])
A screenshot to show result:
It is definitely NOT a bioinformatic code... ^^'
Here is my own example. Although simplistic, I have used the following FASTA parser numerous times. Coded quick and dirty in Python, it produces a list of lists. Each of the inner lists contains two strings, one with the sequence name and the other with the sequence itself. Here it is:
def read_fasta(path):
"""
Read FASTA file into a list of lists.
Each inner list contains a name and a sequence.
"""
out = []
last_record = -1
with open(path) as file_:
for line in file_:
if line.startswith(">"):
contig_name = line.split()[0]
contig_seq = ""
out.append([contig_name, contig_seq])
else:
out[last_record][1] += line.rstrip()
return out
Cheers!
FYI: for a high performance python based Fasta library see PyFasta, http://pypi.python.org/pypi/pyfasta/ I for one wish I heard of this package sooner.
We can be fairly sure the code is tested - only that this is not the right place to show/distribute them. As for FASTA parsing, indeed there are many choices, though the vast majority of them have substantial limitations. Although a seemingly simple task, the problem of parsing a FASTA file not nearly as easy to solve as most people assume. Take your FASTA reader; use it to extract the last 100 bases of chrom1 of the human genome. Now use the same FASTA reader to count 50 million short reads with it. You will soon see that most implementation fail spectacularly at one or both of these tasks.
I see very justified critics concerning my FASTA script. Maybe I can offer an explanation as to why it figures here :) The main reason why this is my FAVORITE (which is the question) is that this specific script has been at the core of my Python learning process. Python has opened tons of interesting applications for me in my everyday work. It is also an example of how doing complex things (here at best 'lightly' complex) is rendered easy and even beautiful with some scripting in your favorite language. I do not even suggest you use this. As mentioned, there are dependable parsers out there :)
For batch processing I like to use GNU parallel, especially for single-thread tools. For instance:
parallel --results gene_counts -j 8 htseq-count -f bam -r pos {} my.gtf ::: *bam
It stores in gene_counts/1 for each bam-file the htseq-count sdtout and stderr output.
sed -i 's/>chr/>Chr/g' genome.fasta
Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
version 2.3.6
you should ask for any neat code not just Python
You are probably right. This is a bit narrow for a forum of bioinformaticiant regrouping 130 people :)
I thought this would be about preferred acronyms for bioinformatics; mine's bfx