Последние несколько недель я пытался загрузить набор геномных данных в Neo4j, используя библиотеку Scikit Allel. Мне удалось загрузить все варианты экзомов в файл VCF и всех субъектов с соответствующими данными фенотипа, теперь я просто пытаюсь создать отношения между вариантами и субъектами. Я не очень разбираюсь в python, и я не думаю, что этот вопрос требует хорошего понимания ни геномики, ни библиотеки Scikit-Allel, так что не пугайтесь этого.
Код ниже работает, но невероятно медленно. Я думаю, что создание фрейма данных или наборов списков для каждого субъекта может быть лучшим способом сделать это, а не создавать транзакцию графа для каждого элемента в цикле j for, но любые советы о том, как лучше всего это сделать, будут оценены. Я также рассматривал возможность включения Numba для этого, но не знал, с чего начать.
Этот процесс создает ~ 1 нового субъекта в час на каждую хромосому из ~ 1750 субъектов и 23 хромосом, поэтому для того, чтобы текущая настройка заработала, потребуется очень много времени.
import pandas as pd
import csv
import math
import allel
import zarr
from py2neo import Graph, Node, Relationship, NodeMatcher
zarr_path = '.../chroms.zarr'
callset = zarr.open_group(zarr_path, mode='r')
samples = callset[chrom]['samples']
graph = Graph(user="neo4j", password="password")
chrom_list = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,'X']
for chrom in chrom_list:
variants = allel.VariantChunkedTable(callset[chrom]['variants'], names=['AC','AF_AFR', 'AF_AMR', 'AF_ASN', 'AF_EUR', 'AF_MAX', 'CGT', 'CLR', 'CSQ', 'DP', 'DP4', 'ESP_MAF', 'FILTER_LowQual', 'FILTER_MinHWE', 'FILTER_MinVQSLOD', 'FILTER_PASS', 'HWE', 'ICF', 'ID', 'IS', 'PC2', 'PCHI2', 'POS', 'PR', 'QCHI2', 'QUAL', 'REF', 'ALT', 'INDEL', 'SHAPEIT', 'SNP_ID', 'TYPE', 'UGT', 'VQSLOD', 'dbSNPmismatch', 'is_snp', 'numalt'], index='POS')
pos = variants['POS'][:]
SNPid = variants['ID'][:]
ref = variants['REF'][:]
alt = variants['ALT'][:]
dp = variants['DP'][:]
ac = variants['AC'][:]
vartype = variants['TYPE'][:]
qual = variants['QUAL'][:]
vq = variants['VQSLOD'][:]
numalt = variants['numalt'][:]
csq = variants['CSQ'][:]
vcfv = 'VCFv4.1'
refv = 'file:///lustre/scratch105/projects/g1k/ref/main_project/human_g1k_v37.fasta'
dpz = callset[chrom]['calldata/DP']
psz = callset[chrom]['calldata/PS']
plz = callset[chrom]['calldata/PL']
gpz = callset[chrom]['calldata/GP']
calldata = callset[chrom]['calldata']
gt = allel.GenotypeDaskArray(calldata['GT'])
hap = gt.to_haplotypes()
hap1 = hap[:, ::2]
hap2 = hap[:, 1::2]
i = 0
for i in range(len(samples)):
subject = samples[i]
subject_node = matcher.match("Subject", subject_id= subject)
if subject_node.first() is None:
continue
seq_tech = 'Illumina HiSeq 2000 (ILLUMINA)'
dp = dpz[:, i]
ps = psz[:, i]
pl = plz[:, i]
gp = gpz[:, i]
list_h1 = hap1[:, i].compute()
list_h2 = hap2[:, i].compute()
chrom_label = "Chromosome_" + str(chrom)
j = 0
for j in range(len(pos)):
h1 = int(list_h1[j])
h2 = int(list_h2[j])
read_depth = int(dp[j])
ps1 = int(ps[j])
PL0 = int(pl[j][0])
PL1 = int(pl[j][1])
PL2 = int(pl[j][2])
genotype = str(h1) + '|' + str(h2)
GP0 = float(gp[j][0])
GP1 = float(gp[j][1])
GP2 = float(gp[j][2])
k = int(pos[j])
l = str(ref[j])
m = str(alt[j][h1-1])
o = str(alt[j][h2-1])
if h1 == 0 and h2 == 0:
a1 = matcher.match(chrom_label, "Allele", pos= k, bp = l)
r2 = Relationship(subject_node.first(), "Homozygous", a1.first(), HTA=h1, HTB=h2, GT=genotype, seq_tech=seq_tech, dp=read_depth, phase_set=ps1, PL0=PL0, PL1=PL1, PL2=PL2, GP0=GP0, GP1=GP1, GP2=GP2)
graph.create(r2)
elif h1 == 0 and h2 > 0:
a1 = matcher.match(chrom_label, "Allele", pos= k, bp = l)
r2 = Relationship(subject_node.first(), "Heterozygous", a1.first(), HTA=h1, GT=genotype, seq_tech=seq_tech, dp=read_depth, phase_set=ps1, PL0=PL0, PL1=PL1, PL2=PL2, GP0=GP0, GP1=GP1, GP2=GP2)
graph.create(r2)
a2 = matcher.match(chrom_label, "Allele", pos= k, bp = o)
r3 = Relationship(subject_node.first(), "Heterozygous", a2.first(), HTB=h2, GT=genotype, seq_tech=seq_tech, dp=read_depth, phase_set=ps1, PL0=PL0, PL1=PL1, PL2=PL2, GP0=GP0, GP1=GP1, GP2=GP2)
graph.create(r3)
elif h1 > 0 and h2 == 0:
a1 = matcher.match(chrom_label, "Allele", pos= k, bp = m)
r2 = Relationship(subject_node.first(), "Heterozygous", a1.first(), HTA=h1, GT=genotype, seq_tech=seq_tech, dp=read_depth, phase_set=ps1, PL0=PL0, PL1=PL1, PL2=PL2, GP0=GP0, GP1=GP1, GP2=GP2)
graph.create(r2)
a2 = matcher.match(chrom_label, "Allele", pos= k, bp = l)
r3 = Relationship(subject_node.first(), "Heterozygous", a2.first(), HTB=h2, GT=genotype, seq_tech=seq_tech, dp=read_depth, phase_set=ps1, PL0=PL0, PL1=PL1, PL2=PL2, GP0=GP0, GP1=GP1, GP2=GP2)
graph.create(r3)
elif h1 == h2 and h1 > 0:
a1 = matcher.match(chrom_label, "Allele", pos= k, bp = m)
r2 = Relationship(subject_node.first(), "Homozygous", a1.first(), HTA = h1, HTB = h2, GT=genotype, seq_tech=seq_tech, dp=read_depth, phase_set=ps1, PL0=PL0, PL1=PL1, PL2=PL2, GP0=GP0, GP1=GP1, GP2=GP2)
graph.create(r2)
else:
a1 = matcher.match(chrom_label, "Allele", pos= k, bp = m)
r2 = Relationship(subject_node.first(), "Heterozygous", a1.first(), HTA=h1, GT=genotype, seq_tech=seq_tech, dp=read_depth, phase_set=ps1, PL0=PL0, PL1=PL1, PL2=PL2, GP0=GP0, GP1=GP1, GP2=GP2)
graph.create(r2)
a2 = matcher.match(chrom_label, "Allele", pos= k, bp = o)
r3 = Relationship(subject_node.first(), "Heterozygous", a2.first(), HTB=h2, GT=genotype, seq_tech=seq_tech, dp=read_depth, phase_set=ps1, PL0=PL0, PL1=PL1, PL2=PL2, GP0=GP0, GP1=GP1, GP2=GP2)
graph.create(r3)
print("Subject " + subject + " completed.")
print(chrom_label + "completed.")
Любая помощь очень ценится!
