python – 用于查找总质量m的可能氨基酸序列的算法优化

因此,这可以简化为线性规划问题 – 找到系数的值,使得M = a * A b * C c * D d * E e * G … r * W

一旦你有了解决方案,你就可以生成给定氨基酸组的所有可能的排列 – 或者如果你只需要排列的数量,你可以直接计算它.

编辑：

正如@Hooked指出的那样,这不是线性规划,原因有二：首先,我们需要整数系数,其次,我们正在寻找所有组合,而不是找到一个单一的最优解.

我已经制定了一个递归生成器,如下所示：

from math import floor,ceil
import profile

amino_weight = {
    'A':  71.038,'C': 103.009,'D': 115.027,'E': 129.043,'F': 147.068,'G':  57.021,'H': 137.059,'I': 113.084,'K': 128.095,'L': 113.084,# you omitted leutine?
    'M': 131.040,'N': 114.043,'P':  97.053,'Q': 128.059,# you omitted glutamine?
    'R': 156.101,'S':  87.032,'T': 101.048,'V':  99.068,'W': 186.079,'Y': 163.063
}

def get_float(prompt):
    while True:
        try:
            return float(raw_input(prompt))
        except ValueError:
            pass

# This is where the fun happens!
def get_mass_combos(aminos,pos,lo,hi,cutoff):
    this = aminos[pos]         # use a pointer into the string,to avoid copying 8 million partial strings around
    wt = amino_weight[this]
    kmax = int(floor(hi / wt))
    npos = pos - 1
    if npos:                   # more aminos to consider recursively
        for k in xrange(0,kmax + 1):
            mass    = k * wt
            nlo     = lo - mass
            nhi     = hi - mass
            ncutoff = cutoff - mass
            if nlo <= 0. and nhi >= 0.:
                # we found a winner!
                yield {this: k}
            elif ncutoff < 0.:
                # no further solution is possible
                break
            else:
                # recurse
                for cc in get_mass_combos(aminos,npos,nlo,nhi,ncutoff):
                    if k > 0: cc[this] = k
                    yield cc
    else:                      # last amino - it's this or nothing
        kmin = int(ceil(lo / wt))
        for k in xrange(kmin,kmax+1):
            yield {this: k}

def to_string(combo):
    keys = sorted(combo)
    return ''.join(k*combo[k] for k in keys)

def total_mass(combo):
    return sum(amino_weight[a]*n for a,n in combo.items())

def fact(n):
    num = 1
    for i in xrange(2,n+1):
        num *= i
    return num

def permutations(combo):
    num = 0
    div = 1
    for v in combo.values():
        num += v
        div *= fact(v)
    return fact(num) / div

def find_combos(lo,hi):
    total = 0
    bases = []
    aminos = ''.join(sorted(amino_weight,key = lambda x: amino_weight[x]))
    for combo in get_mass_combos(aminos,len(aminos)-1,hi - amino_weight[aminos[0]]):
        base = to_string(combo)
        bases.append(base)
        mass = total_mass(combo)
        cc = permutations(combo)
        total += cc
        print("{} (mass {},{} permutations)".format(base,mass,cc))
    print('Total: {} bases,{} permutations'.format(len(bases),total))

def main():
    lo = get_float('Bottom of target mass range? ')
    hi = get_float('Top of target mass range? ')

    prof = profile.Profile()
    prof.run('find_combos({},{})'.format(lo,hi))
    prof.print_stats()

if __name__=="__main__":
    main()

它还使用浮点氨基质量来寻找质量范围.在我的机器(i5-870)上搜索748.0和752.0之间的质量,返回7,505个碱基,总共9,400,528个排列,在3.82秒内.