Beginning Python for Bioinformatics

Packt Publishing invited me to review Django 1.1 Testing and Debuggingby Karen M. Tracey. They also kindly provided a free chapter that you can download from the link below. A full review will be posted as soon as I finish the book.

preview chapter – Chapter No.3 “Testing 1, 2, 3: Basic Unit Testing”

Last week I made a 5 book order at Amazon and one of them was Bioinformatics Programming Using Python: Practical Programming for Biological Data (Animal Guide) by Mitchell L Model.

I started reading the book late Friday night, and I’m on the third chapter, where there’s an introduction to sequences. So far, I found the book very confusing, especially as it claims to be a book for people with no programming background. The examples are OK, but there’s a very messy mixture of Python interpreter and standalone script usage, as the author jumps back and forth them. Another thing is that some examples are explained in detail including the line number, while others you depend on the code’s docstring to understand it.

So far, I’m not impressed. The initial Python sequence example is a set and in this chapter there already are some functional programming concepts, what can be quite challenging to someone that has never programmed in their life. And in the second chapter the reader sees a ternary operator. Another criticism, is that in the preface the author suggests using Python 3, instead of 2, what might add to the frustration of the beginner when a module cannot be installed.

I will continue reading it and post whenever I have a more complete overview of the book.

But, I guess, it can be easily refactored/enhanced/despised by the audience that read or have access to this blog via Planet Python. Anyway, for someone like me, whose main task now is not to generate tons of code and lines, I think the code (or part of it) that I will present below is quite good. Feel free to comment, criticize and say bad and good things about it.

We needed a script that would take files coming out from protein search engines that would be able to compare the peptides and protein sequences, their abundance and some other characteristics. We had a combination of protein and peptide files, with a list of proteins (one protein per line in a tab delimited file) that was related to a list of peptides in another file (one peptide per line, with multiple peptides/lines related to one protein in the original list). Also each line in both files had more than 50 columns, and 8 or 16 of them were the values we wanted to extract. I say 8 or 16 of them because we didn’t know how many will output each time, as it would depend on the number of samples per run (4 to 8 samples) .

So, we had a couple of issues: we didn’t know how many proteins would be output (actually found) in each file, we wouldn’t know how many peptides for each protein would be found and we didn’t know before hand how many samples would be run at once. One good thing is that the 8-16 columns of values were fixed, always in the same position and with empty cells if no value was registered there. And we had a fourth problem: usually the samples attributions would be random, meaning a control could come in the first value column or could come in the last. And a fifth as we didn’t know before hand (the tech knew) how many treatments would be run each time. A treatment could be a different experimental condition, a sample grouping or some other extraneous factor. An extra issue is that we would need to compare multiple files, and get protein and peptide abundances in all of them at the same time and finally compare each treatment.

Basically, in order to create an universal script we needed something flexible enough that whatever the experiments threw at use we would be able to handle. First step we decided to use a YAML file that could be filled by the experimental researchers with sample assignments, treatments, etc. The YAML would look like this

B0:
– 114: A
– 115: D
– 116: B
– 117: C

B1:
– 114: C
– 115: A
– 116: D
– 117: B

In this file B0 and B1 would be the result file names, 114 is the column/channel where the sample was run and and A, B, C and D the treatment. With this set, out objective was to get all proteins and their peptides for treatment A in files B0 and B1, do some calculation and them compare to all proteins and peptides from treatment B, C and D extracted also from files B0 and B1.

First step was to get the names of the treatments from the YAML file

]def get_treatments(mapping):
    treats = set([])
    for entry in mapping:
        [treats.add(list(t.values())[0]) for t in mapping[entry]]

    return treats

where [code]mapping[/code] is the YAML file name. We used a set to store and sets have unique items, and treatment names can vary from file to file. In the code above we basically read the YAML and the value for each entry.

We then needed a class to store protein information, and there was the story got hairy. With all my (lack of ) experience, I decided to use [code]exec[/code] statements to fix all the uncertainty of the experimental data details. I didn’t have the treatment names before hand (or in a fixed immutable list), and didn’t have the columns (channels) that were being used at the time and I have to correctly assign each protein abundance (area) to its place. In the end our class look like this

class Protein():
    """Class Protein, stores all the information about channels and areas, name and accession"""
    def __init__(self, accession, name, treatments):
        self.accession = accession
        self.name = name

        #ratio channels are called based on their name
        for i in treatments:
            exec('self.%s = []' %i)
            exec('self.area%s = []' %i)

    def add_to_channel(self, channel, peptide):
        exec('self.%s.append(peptide)' % channel)

    def add_to_area(self, channel, area):
        exec('self.area%s.append(area)' % channel)

In order to be faithful to this blog’s name, I will explain how the code above is supposed to work. First, [code]exec[/code] is a Python statement that support dynamic execution of code. In our case above it was used to name the objects, so we would be able to access them by name in subsequent functions. Let’s take this for example

for i in treatments:
    exec('self.%s = []' %i)
    exec('self.area%s = []' %i)

In this snippet we were trying to create lists called (for the YAML file above) A, B, C and D, and another set of lists called areaA, areaB, areaC and areaD. Let’s say for another experiment we would have treatments “Control”, “Low” and “High” and so on.

The next two functions use the same approach, with exec, this time appending to the freshly created lists. This way it’s easy to control what the user is throwing at us.

I don’f know if this the best approach possible, or if it is or not harmful. Maybe experts reading this might have better ideas, and I appreciate them. We check the rest of the script next time.

I was invited by Packt Publishing to review Python Testing Beginner’s Guide by Daniel Arbuckle. This is a book on one of the most important aspects of scientific programming (even though the majority of scientific software don’t have any testing routines): code testing, checking if your code actually does what is intended to do. I can say I’m not really an expert on testing so I guess I’m the right audience for it:

You’ll learn about several of Python’s automated testing tools, and you’ll learn about the philosophies and methodologies that they were designed to support, like unit testing and test-driven development. When you’re done, you’ll be able to produce thoroughly tested code faster and more easily than ever before, and you’ll be able to do it in a way that doesn’t distract you from your “real” programming.

Packt also supplied a preview/sample chapter that you can download here.

I hope to get a review ready by the end of the week. before the Ontario Institute of Cancer Research retreat, otherwise I will try to post a full review next week.

Let me start by saying that Python for Bioinformatics (Chapman & Hall/Crc Mathematical & Computational Biology) is a massive book, massive in a way that it contains a lot of material. I still didn’t have enough time to check everything, but I’m well into the first section of the book that gives an initial view of Python and how to program it.

The initial section of the book is well written (I’m not going criticize the book in terms of good/poor English, as I’m not well qualified to do that), and gives a clear perspective on how to program Python for scientists, who are the main target demographic of the book. Of course, it always help to have some basic knowledge of command line shells, but the book also includes some explanations of IDLE and other Python-capable IDEs. I cannot say that I read this section with the enough care and attention, but what I can say is that you won’t miss a beat with PfB, as it has more material than I expected. I still have to start with the more advanced topics, like BioPython and so forth, what I plan to do in the coming month, and as I don’t have a lot of experience with BioPython, I’m looking forward to it.

On the other hand I have a small-ish complaint, that maybe is more about style than substance. I don’t like the design of the book, the way the code interleaves with the text and the way the code explanations are presented. Most of the code blocks are followed by a careful explanation, but this explanation works as a figure label for the code block. That is quite annoying because there are too many stops in the text fluidity as one tends to lose attention to it (my case, not exactly everyone’s).

Another minor detail is the use of “he” every time scientists are referred (one example is on page 3 on the second phrase of the introduction). The (politically) correct would be to use “he or she” or “she or he” (but that’s OK with me).

I will try to post more complete reviews of the sections that I don’t master. I would also like to thank Sebastian for sending me a copy of the book.

I’m slowly moving the posts from the blog to a wiki. It makes easier to display post series and allows people to modify/enhance/discuss.

The wiki address is http://wiki.genedrift.org.

Thanks to the comments and suggestions to the last post, it’s possible to make now a more pythonic and clearly generic database update class. Let’s check how the “generic” update/edit entry function is currently:

def update_data(self, values_list):
    '''edits and updates fields'''

    if sys.platform == 'darwin':
        (cursor, database) = link_db(self.db_path)
    else:
        (cursor, database) = link_db()

    cursor.execute("UPDATE bac SET  projects = ?, comments = ?, temperature = ?, cell = ?, box = ?, tubes = ?, chromosome = ?, sdate = ?, clone = ?, source
	= ?, location1 = ?, startpos = ?, endpos = ?,
	gene = ?, genelink = ?, dnaex = ?, validation = ?, pcr = ?, refs = ?, antibiotic = ? WHERE idbac = ?",
    values_list['projects'], values_list['comments'], values_list['temperature'], values_list['cell'], values_list['box'], values_list['tubes'],
    values_list['chromo'], values_list['date'], values_list['clone'], values_list['source'], values_list['location'], values_list['start']
    values_list['end'], values_list['gene'], values_list['genelink'], values_list['dna'], values_list['validation'], values_list['pcr'],
    values_list['refs'], values_list['antibiotic'], values_list['idbac']))

    database.commit()
    database.close()

which is really ugly and, although it works, is not really useful outside this small project. Based on the comments the best option was to use placeholders and a dictionary, similar to the approach used on the insert data function. Pre-formatting a string to have both the field name to be updated and a placeholder (for instance :idbac) that will receive the values

update = ','.join(['%s=:%s' % (y, y) for y in values_list])

where update is the string we want and values_list is the dictionary with all the key-value pairs. I tried this approach, using this structure in the generic function, but then I decided that the best alternative was to put this join in the derived class function and pre-populate the string with the values and then send this string directly to the update function. In the end I opted to use this

update = ','.join(['%s=\"%s\"' % (y, values_list[y]) for y in values_list])

The latter is slightly different to what was suggested. The original one would create a tuple with the keys from the dictionary, making for instance sdate:sdate. With all these place holders just pass the dictionary and you have all the values inserted. This would be handy if the insert string was being created on the “generic” function. If we move this to the derived class, we can use the the alternative, keeping in mind that the values parsed should be surrounded by quotes, otherwise the SQL UPDATE statement will have problems with spaces and other foreign characters that should not be there. So instead of placeholders we will have gene:"PTEN" and we can attache this joined string to the actual commands. We then can move all the machinery from the “generic” function that can be written as

def update_data(self, update_string):
    '''edits and updates fields'''

    if sys.platform == 'darwin':
        (cursor, database) = link_db(self.db_path)
    else:
        (cursor, database) = link_db()
    cursor.execute(update_string)

    database.commit()
    database.close()

That’s it, very elegant (we will see the derived class in the next post). And finishing our generic class, we would need a delete function, so the user can eliminate entries that he/she doesn’t want anymore. It’s also a very simple function

def delete_data(self, delete_string):
    '''deletes one field'''

    if sys.platform == 'darwin':
        (cursor, database) = link_db(self.db_path)
    else:
        (cursor, database) = link_db()
    cursor.execute(delete_string)

    database.commit()
    database.close()

We will check the delete string next time. Again, I would like to thank for all the comments, it has been really helpful for me.

Previously in the series:
Part 1
Part 2
Part 3
Part 4
Part 5
Part 6
Part 7

And we’re back. After a couple of weeks of inactivity we will get back to our small soap-opera pf Python, wxPython and SQLite. Continuing in our database management code let’s check two other functions that changed since our first inception of the code. The first one is the insert_data function that looks like this now

def insert_data(self, values_list, insert_string):
    '''inserts data in the database'''

    if sys.platform == 'darwin':
        (cursor, database) = link_db(self.db_path)
    else:
        (cursor, database) = link_db()

    cursor.execute(insert_string % self.table_name, values_list)

    database.commit()
    database.close()

Basically no changes, apart from the obvious check for the current running operating system, which was explained in the last post. The other function to check is the update_data. This function is new and it wasn’t in the first version, but as it can be seen it has a problem being a “generic” function, because it contains information pertained to the table and database being used in the interface. This function basically received information that needs to be updated in the table’s fields and by using the SQL UPDATE ... SET edits and updates data in the changed fields. I have tried several different syntaxes to make the execute generic, mainly trying to pre-format the string without success. IF anyone reading this can help, I’d appreciate.

def update_data(self, values_list):
    '''edits and updates fields'''

    if sys.platform == 'darwin':
        (cursor, database) = link_db(self.db_path)
    else:
        (cursor, database) = link_db()

    cursor.execute("UPDATE bac SET  projects = ?, comments = ?, temperature = ?, cell = ?, box = ?, tubes = ?, chromosome = ?, sdate = ?, clone = ?, source = ?, location1 = ?, startpos = ?, endpos = ?,
	gene = ?, genelink = ?, dnaex = ?, validation = ?, pcr = ?, refs = ?, antibiotic = ? WHERE idbac = ?",
    values_list['projects'], values_list['comments'], values_list['temperature'], values_list['cell'], values_list['box'], values_list['tubes'],
    values_list['chromo'], values_list['date'], values_list['clone'], values_list['source'], values_list['location'], values_list['start'],  values_list['end'],
    values_list['gene'], values_list['genelink'], values_list['dna'], values_list['validation'], values_list['pcr'],
    values_list['refs'], values_list['antibiotic'], values_list['idbac']))

    database.commit()
    database.close()

Anyway, I will explain the logic of the command (OK for a stop gap, but not as a definite solution). values_list is a dictionary that is passed to the function and contains the field names as keys and the new/changed information as values. The execute method simply parses the values from each key in the update string which is then sent to the database and table to be changed. Everything is committed and the database is closed.

As this is a “generic” function from a “generic” class the ideal scenario would be to the function to receive a pre-formatted string with all the information, as in the insert data function, and update the information in the database.

I would like to thank in advance anyone that can comment on this. Next time we will continue checking the generic class and finalize this part in order to start with the interface build process.

Previously in the series:
Part 1
Part 2
Part 3
Part 4
Part 5
Part 6

Let’s get back to our SQLite and wxPython project. We haven’t seen anything on wxPython yet, and we will check the interface only on the next post. For now, let’s see some extra code added to the SQLite access class. Remember that we have a generic class and one class derived from it that would work on accessing specific tables in our database file.

When we last covered the db access routines, there was no search for an entry (the function returned everything in the table no matter what), there was no update function in case someone would want to modify an entry and there was no delete method if you wanted to delete something. In the meantime, I added all of this functionality (and some other) to the generic class and extended it to the class derived from it. Let’s check how the generic class is now (you will notice that there is an issue in one of the methods, if someone can help me I’d appreciate. More details later.)

class DB_Generic():
    '''generic class to add DB functionality'''
    def __init__(self, table_name, db_path = ''):
        #par= name of the table to be used
        self.table_name = table_name
        if len(db_path) > 0:
            self.db_path = db_path
            print db_path

    def get_data_generic(self, range = 1, bac_to_get = 0):
        '''gets the data from the database'''       

        if sys.platform == 'darwin':
            (cursor, database) = link_db(self.db_path)
        else:
            (cursor, database) = link_db()

        if range == 1:
            cursor.execute("""SELECT * FROM %s""" % self.table_name)
        elif range == 2:
            cursor.execute("""SELECT * FROM %s where idbac = %d""" % (self.table_name, bac_to_get))

        table_data = cursor.fetchall()
        raw_data = []
        for i in table_data:
            raw_data.append(list(i))

        self.table_data = raw_data
        database.close()

    def insert_data(self, values_list, insert_string):
        '''inserts data in the database'''

        if sys.platform == 'darwin':
            (cursor, database) = link_db(self.db_path)
        else:
            (cursor, database) = link_db()

        cursor.execute(insert_string % self.table_name, values_list)

        database.commit()
        database.close()

    def update_data(self, values_list):
        '''edits and updates fields'''

        if sys.platform == 'darwin':
            (cursor, database) = link_db(self.db_path)
        else:
            (cursor, database) = link_db()

        #change this to generic!!!!!!!!!!!!
        cursor.execute("UPDATE bac SET  projects = ?, comments = ?, temperature = ?, cell = ?, box = ?, tubes = ?, chromosome = ?, sdate = ?, clone = ?, source = ?, location1 = ?, startpos = ?, endpos = ?,
		gene = ?, genelink = ?, dnaex = ?, validation = ?, pcr = ?, refs = ?, antibiotic = ? WHERE idbac = ?",
        (values_list['projects'], values_list['comments'], values_list['temperature'], values_list['cell'], values_list['box'], values_list['tubes'],
         values_list['chromo'], values_list['date'], values_list['clone'], values_list['source'], values_list['location'], values_list['start'], values_list['end'],
         values_list['gene'], values_list['genelink'], values_list['dna'], values_list['validation'], values_list['pcr'],
         values_list['refs'], values_list['antibiotic'], values_list['idbac']))

        database.commit()
        database.close()

    def delete_data(self, delete_string):
        '''deletes one field'''

        if sys.platform == 'darwin':
            (cursor, database) = link_db(self.db_path)
        else:
            (cursor, database) = link_db()
        cursor.execute(delete_string)

        database.commit()
        database.close()

In the next couple of posts we’ll dissect each function and see what’s going on. The class definition wasn’t changed, so we start with get_data_generic

def get_data_generic(self, range = 1, bac_to_get = 0):
	'''gets the data from the database'''       

	if sys.platform == 'darwin':
		(cursor, database) = link_db(self.db_path)
	else:
		(cursor, database) = link_db()

	if range == 1:
		cursor.execute("""SELECT * FROM %s""" % self.table_name)
	elif range == 2:
		cursor.execute("""SELECT * FROM %s where idbac = %d""" % (self.table_name, bac_to_get))

	table_data = cursor.fetchall()
	raw_data = []
	for i in table_data:
		raw_data.append(list(i))

	self.table_data = raw_data
	database.close()

The first difference we notice here is the sys.platform usage. This is required if we intend to package our application as an OS X app, using py2app. When a Python/wxPython application is packaged in OS X, the actual application executable is inside the a directory named after the application (or whatever you set up). In our case here we don’t provide a way for the Python script to receive the path and name for the database on a command line, as we expect it to be in the executable’s current directory. Because of that we need to provide a “config” file (in our case here a one-line text file with the database path) inside the application wrapper, something we will see in the end of the series.

Another modification here is the range parameter and the addition of the bac_to_get parameter. Notice that both parameters have a value assigned to it. This means that they are optional, the function’s call can pass them or not. If it doesn’t pass, their value will be the one assigned on the function declaration. So, here if we are interested in getting all bacs, range will have the value of 1 and we don’t need to worry about it. If we want an specific bac we will pass range as 2 and then pass the bac_to_get ID to be returned.

A final change/addition is that we added a new select statement for the cases when range equals 2. This time we are adding the bac ID to be returned.

Previously in the series:
Part 1
Part 2
Part 3
Part 4
Part 5

Just illustrating my point (or lack of), an animation about the commits of RoR to its repository. Notice the jump after it was migrated to Github

Ruby on Rails from Ilya Grigorik on Vimeo.

Sorry for the non-Python post.