Sorting genomic alignments using Python

Apr 3, 2013 3 min read

Originally published on my old Blogger blog on April 3, 2013. Read the original post and its comments.

http://bcbio.wordpress.com/2009/07/26/sorting-genomic-alignments-using-python/

Blue Collar Bioinformatics

Sorting genomic alignments using Python

Recent, I was asked to sort a MAF (Multiple Alignment Format) file containing genomic alignments. This was a large vertebrate alignment file and the researchers were interested in manually examining the longest matches. The tricky part of doing this is that the entire file cannot comfortably fit in memory on a standard machine.
The common solution to dealing with large memory intensive files is to pre-parse the file and provide a way to look up individual records based on an identifier. For instance, items could be pushed into a relational database table and retrieved based on primary keys. While moving to disk access is less efficient, it allows your code to scale to these type of real life problems without resorting to buying large memory machines.
We will use the bx-python library to read the MAF alignment files, prepare an index for accessing the alignments individually, and ultimately perform the sorting.
The first step is to define the index retrieval method. bx-python has index retrival functionality that allows querying and retrieving records based on genomic location. We don’t need that advanced query capability here, and can retrieve records based on position in the file as well. We build the index by looping through each record and recording the file position and genomic locations in the index.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

from bx.alignimport maf
from bximport interval_index_file

def build_index(in_file, index_file):
    indexes= interval_index_file.Indexes()
    withopen(in_file) as in_handle:
        reader= maf.Reader(in_handle)
        while 1:
            pos= reader.file.tell()
            rec= reader.next()
            if recis None:
                break
            for cin rec.components:
                indexes.add(c.src, c.forward_strand_start,
                        c.forward_strand_end, pos,max=c.src_size )

withopen(index_file,"w") as index_handle:
indexes.write(index_handle)

Now we have a unique index to retrieve a record — the position in the file. The second step is to loop through the file again and build a list of alignment sizes and positions. Alignment size is the first value in each list item, since this is the criteria to sort by. This list easily fits in memory, and we sort it with the largest alignments first.

1
2
3
4
5
6
7
8
9
10

rec_info= []
withopen(in_file) as in_handle:
    reader= maf.Reader(in_handle)
    while 1:
        pos= reader.file.tell()
        rec= reader.next()
        if recis None:
            break
        rec_info.append((rec.text_size, pos))
rec_info.sort(reverse=True)

Finally, we loop though the sorted list of sizes and use the associated positions to get records from our index. Each record is then written to an output file.

1
2
3
4
5
6

index= maf.Indexed(in_file, index_file)
withopen(out_file,"w") as out_handle:
    writer= maf.Writer(out_handle)
    for size, posin rec_info:
        rec= index.get_at_offset(pos)
        writer.write(rec)

The full script puts this together into a usable program. This could be adjusted to deal with other alignment formats supported by bx-python, like axt files.

genomics python

Sorting genomic alignments using Python

http://bcbio.wordpress.com/2009/07/26/sorting-genomic-alignments-using-python/

Blue Collar Bioinformatics

Sorting genomic alignments using Python

Related