Как я могу улучшить свой цикл сценариев Python, включая различные математические операции для разных условий в каждом цикле?

Question

Я пишу снова, потому что мне не повезло, и я попытался сделать следующий скрипт более эффективным. Для более подробной информации, посмотрите мой предыдущий пост , но основная ситуация такая, как показано ниже.

Я написал скрипт для подсчета очков, а также частоту для списка генетических профилей.

Генетический профиль здесь состоит из комбинации SNP. Каждый SNP имеет два аллеля. Следовательно, входной файл для 3 SNP выглядит примерно так, как показано ниже, где показаны все возможные комбинации всех аллелей для всех трех SNP. Эта таблица была сгенерирована с использованием продукта itertool в другом скрипте:

    AA   CC   TT
    AT   CC   TT
    TT   CC   TT
    AA   CG   TT
    AT   CG   TT
    TT   CG   TT
    AA   GG   TT
    AT   GG   TT
    TT   GG   TT
    AA   CC   TA
    AT   CC   TA
    TT   CC   TA
    AA   CG   TA
    AT   CG   TA
    TT   CG   TA
    AA   GG   TA
    AT   GG   TA
    TT   GG   TA
    AA   CC   AA
    AT   CC   AA
    TT   CC   AA
    AA   CG   AA
    AT   CG   AA
    TT   CG   AA
    AA   GG   AA
    AT   GG   AA
    TT   GG   AA

Затем у меня есть другой файл с таблицей, содержащей веса и частоты для трех SNP, как показано ниже:

SNP1             A       T       1.25    0.223143551314     0.97273 
SNP2             C       G       1.07    0.0676586484738    0.3     
SNP3             T       A       1.08    0.0769610411361    0.1136  

Столбцами являются идентификаторы SNP, аллель риска, эталонный аллель, OR, log (OR) и частота населения. Вес дан для аллеля риска.

Основной сценарий берет эти два файла и вычисляет оценку на основе суммы коэффициентов логарифмов для каждого аллеля риска в каждом SNP для каждого генетического профиля, а также частоты, основанной на умножении частот аллелей, предполагая, что Hardy Равновесие Вайнберга.

import sys

snp={}
riskall={}
weights={}
freqs={}    # effect allele, *MAY NOT BE MINOR ALLELE

pop = int(int(sys.argv[4]) + 4) # for additional columns due to additional populations. the example table given only has one population (column 6)

# read in OR table
pos = 0
with open(sys.argv[1], 'r') as f:
    for line in f:
        snp[pos]=(line.split()[0])
        riskall[line.split()[0]]=line.split()[1]
        weights[line.split()[0]]=line.split()[4]
        freqs[line.split()[0]]=line.split()[pop]

        pos+=1



### compute scores for each combination
with open(sys.argv[2], 'r') as f:
    for line in f:
        score=0
        freq=1
        for j in range(len(line.split())):
            rsid=snp[j]
            riskallele=riskall[rsid]
            frequency=freqs[rsid]
            wei=weights[rsid]
            allele1=line.split()[j][0]
            allele2=line.split()[j][1]
            if allele2 != riskallele:      # homozygous for ref
                score+=0
                freq*=(1-float(frequency))*(1-float(frequency))
            elif allele1 != riskallele and allele2 == riskallele:  # heterozygous, be sure that A2 is risk allele!
                score+=float(wei)
                freq*=2*(1-float(frequency))*(float(frequency))
            elif allele1 == riskallele: # and allele2 == riskall[snp[j]]:      # homozygous for risk, be sure to limit risk to second allele!
                score+=2*float(wei)
                freq*=float(frequency)*float(frequency)

            if freq < float(sys.argv[3]):   # threshold to stop loop in interest of efficiency 
                break

        print(','.join(line.split()) + "\t" + str(score) + "\t" + str(freq))

Я установил переменную, в которой я могу указать порог для разрыва цикла, когда частота становится крайне низкой. Какие улучшения могут быть сделаны для ускорения работы скрипта?

Я пытался использовать Pandas, который все еще намного медленнее, так как я не уверен, возможна ли в этом случае векторизация. У меня проблемы с установкой Dask на мой сервер Unix. Я также позаботился о том, чтобы использовать только словари Python, а не списки, и это дало небольшое улучшение.

Ожидаемый результат из вышеперечисленного будет таким:

GG,AA,GG        0       0.000286302968304
GG,AA,GA        0.0769610411361 7.33845153414e-05
GG,AA,AA        0.153922082272  4.70243735491e-06
GG,AG,GG        0.0676586484738 0.00024540254426
GG,AG,GA        0.14461968961   6.29010131498e-05
GG,AG,AA        0.221580730746  4.03066058992e-06
GG,GG,GG        0.135317296948  5.25862594844e-05
GG,GG,GA        0.212278338084  1.34787885321e-05
GG,GG,AA        0.28923937922   8.63712983555e-07
GA,AA,GG        0.223143551314  0.0204250448374
GA,AA,GA        0.30010459245   0.00523530030129
GA,AA,AA        0.377065633586  0.000335475019306
GA,AG,GG        0.290802199788  0.0175071812892
GA,AG,GA        0.367763240924  0.00448740025824
GA,AG,AA        0.44472428206   0.000287550016548
GA,GG,GG        0.358460848262  0.00375153884769
GA,GG,GA        0.435421889398  0.000961585769624
GA,GG,AA        0.512382930534  6.16178606889e-05
AA,AA,GG        0.446287102628  0.364284082594
AA,AA,GA        0.523248143764  0.0933724543834
AA,AA,AA        0.6002091849    0.00598325294334
AA,AG,GG        0.513945751102  0.312243499367
AA,AG,GA        0.590906792238  0.0800335323286
AA,AG,AA        0.667867833374  0.00512850252286
AA,GG,GG        0.581604399576  0.0669093212928
AA,GG,GA        0.658565440712  0.0171500426418
AA,GG,AA        0.735526481848  0.00109896482633

РЕДАКТИРОВАТЬ: Добавлена ​​ссылка на предыдущий пост, а также ожидаемый результат.

Answer 1

Отказ от ответственности: я не проверял это, это скорее псевдокод.

Я даю некоторые общие идеи о том, что медленно / быстро в программировании и особенно в python:

Вы должны попытаться вывести из циклов все, что не меняется в этом цикле. Кроме того, в Python, вы должны попытаться заменить циклы с пониманием https://www.pythonforbeginners.com/basics/list-comprehensions-in-python

[ expression for item in list if conditional ]

вы должны попытаться использовать функции карты / фильтра, если это возможно, и вы также можете подготовить свои данные, чтобы программа была более эффективной

    rsid=snp[j]
    riskallele=riskall[rsid]

по сути является двойным отображением, и, возможно, его можно будет сделать лучше, если вы сможете создать свою структуру snp следующим образом (вы можете использовать -1 индекс для последнего столбца и избавиться от pop):

snp = [{"riskall": line[1],"freq": float(line[4]),"weight": float(line[-1])} 
         for line in map(split,f)]

и ваш вычислительный цикл может выглядеть примерно так:

### compute scores for each combination
stop = sys.argv[3]
with open(sys.argv[2], 'r') as f:
    for fline in f:
        score=0.0 # work with floats from the start
        freq=1.0
        line = fline.split() # do it only once

        for j,field in line:
            s=snp[j]
            riskallele=s["riskall"]
            frequency=s["freq"]
            wei=s["weight"]
            (allele1,allele2) = line[j]

            if allele2 != riskallele:      # homozygous for ref
                score+=0
                freq*=(1-frequency)*(1-frequency)
            elif allele1 != riskallele and allele2 == riskallele:  # heterozygous, be sure that A2 is risk allele!
                score+=wei
                freq*=2*(1-frequency)*frequency
            elif allele1 == riskallele: # and allele2 == riskallele:      # homozygous for risk, be sure to limit risk to second allele!
                score+=2*wei
                freq*=frequency*frequency

            if freq < stop):   # threshold to stop loop in interest of efficiency 
                break

        print(','.join(line.split()) + "\t" + str(score) + "\t" + str(freq))

Конечная цель, которую я хотел бы достичь, - преобразовать ее в какую-то карту / уменьшить форму:

аллель может иметь [A, C, G, T] [A, C, G, T] 16 комбинаций, и мы проверяем на это [A, C, G, T] это только 64 комбинации, поэтому я могу создать карта в форме [AC, C] -> Score, freq_function и я могу избавиться от всего блока if.

Иногда лучшим подходом является разделение кода на небольшие функции, реорганизация и последующее слияние.