Script to extract FASTA sequences

＃!/usr/bin/env python

import sys
import mmap
import re

indexes = {}

# read the index file first
infile = open("%s.fai" % sys.argv[1])
for line in infile:
    line = line.strip()
    if line != "": 
        fields = line.split()
        if len(fields) != 5:
            print "Incorrect index line: %s" % line
            sys.exit(1)
        indexes[fields[0]] = {
            "length" : int(fields[1]),
            "offset" : int(fields[2]),
            "base_per_row" : int(fields[3]),
            "chr_per_row" : int(fields[4]), 
        }
infile.close()

def posToOffset(chr, pos):
    pos -= 1 # change to 0 based index
    offset = indexes[chr]["offset"] # use the chromosome offset first
    offset += pos # add the number of based
    offset += (pos / indexes[chr]["base_per_row"]) * (indexes[chr]["chr_per_row"] - indexes[chr]["base_per_row"]) # account for newlines
    ＃print chr, pos, indexes[chr]["offset"], offset
    return offset

with open(sys.argv[1], "r+b") as f:
    mm = mmap.mmap(f.fileno(), 0)
    for region in sys.argv[2:]:
    # range
    mo = re.search("^(\\w+)\\:(\\d+)\\-(\\d+)$", region)
    if mo != None:
        chr = mo.group(1)
        if chr in indexes:
            start = int(mo.group(2))
            end = int(mo.group(3))
            spos = posToOffset(chr, start)
            epos = posToOffset(chr, end+1)
            seq = mm[spos:epos].replace("\n", "")
            if len(seq) != (end-start+1):
                sys.stderr.write("Incorrect sequence length found for %s, %s\n" % (region, seq))
            else:
                print seq else: sys.stderr.write("Unknown chromosome %s\n" % chr)
    else:
        # single position
        mo = re.search("^(\\w+)\\:(\\d+)$", region)
        if mo != None:
            chr = mo.group(1)
            if chr in indexes:
                spos = posToOffset(chr, int(mo.group(2)))
                print mm[spos].replace("\n", "")
            else:
                sys.stderr.write("Unknown chromosome %s\n" % chr)
        else:
            sys.stderr.write("Unknown region %s\n" % region)