Я написал сценарий для вычисления балла, а также частоты для списка генетических профилей.
Генетический профиль здесь состоит из комбинации 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.Теперь я изменил все структуры данных на словари, насколько это возможно.Что еще я могу сделать?