Повышение эффективности зацикливания, математический скрипт - PullRequest
0 голосов
/ 16 октября 2018

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

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

    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        RiskAll  RefAll         OR            log(OR)    RiskAllFreq) # example header, not in file
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 в генетическом профиле, а также частоты, основанной на умножении частот аллелей, при условии равновесия Харди Вайнберга.

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))

До сих пор я установил переменную, в которой я могу указать порог для прерывания цикла, когда частота становится очень низкой, например, приблизительно 1e-10.Я хочу увеличить это, чтобы включить по крайней мере 20 SNP.Какие улучшения могут быть сделаны для ускорения работы скрипта?

РЕДАКТИРОВАТЬ: Добавлен пример файла таблицы.

РЕДАКТИРОВАТЬ 2: Я сейчас попытался запустить скрипт с порогом частоты 1e-4.Прошло уже около шести дней, и он все еще работает, и это слишком медленно, поэтому я ищу дополнительные предложения!

РЕДАКТИРОВАТЬ 3: Уточненные входные файлы, заголовки на самом деле не во входном файле, онипросто указание.

РЕДАКТИРОВАТЬ 4: Попытка с использованием панд, но это было медленно, и не уверен, что для этой цели возможна векторизация.У Dask есть проблемы с установкой на мой сервер Unix.Теперь я изменил все структуры данных на словари, насколько это возможно.Что еще я могу сделать?

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...