0

我正在尝试访问基于 Numpy 的用于基因组数据的工具SciKit Allele中的基因组数据。

我对 python 不是很好，但我试图遍历每个变体并提取数组中的相关列，然后使用 Neo4j Rest Client 在 Neo4j 数据库中创建节点。

下面的代码生成一个包含所有变体和所有数据类型的数组：

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', 'svlen'], index='POS')

我（我想我）已经以数组形式声明了变量，如下所示：

pos = variants['POS'][:]
alt = variants['ALT'][:]
dp = variants['DP'][:]
ac = variants['AC'][:]
type = variants['TYPE'][:]
svlen = variants['svlen'][:]
qual = variants['QUAL'][:]
vq = variants['VQSLOD'][:]

这些变量创建数组，例如：

In: pos
Out: array([    28590,     50481,     52152, ..., 249225077, 249229702,
       249231222], dtype=int32)

我现在正在尝试访问每一行的变量，但似乎无法弄清楚如何做到这一点。我当前的尝试看起来像这样（前 10 行）：

for variant in variants[0:10]:
    a1 = db.nodes.create(pos=pos[variant], bp=alt[variant][0], DP=dp[variant], AC=ac[variant][0], type=type[variant][0], svlen=svlen[variant][0], qual=qual[variant], vqslod=vq[variant])
    a1.relationships.create("belongs_to", c1)

不幸的是，这会出现以下错误：

IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices

谁能帮我弄清楚如何获取每个属性的特定变量？

4

1 回答 1

1

由于没有可重现的示例，我无法遵循您的代码，因此我必须根据 scikit-allel 文档创建一个：

https://scikit-allel.readthedocs.io/en/stable/model/chunked.html#variantchunkedtable

import h5py
import allel
import os

#cleanup
h5_file = 'callset.h5'
os.remove(h5_file) if os.path.exists(h5_file) else None

chrom = [b'chr1', b'chr1', b'chr2', b'chr2', b'chr3', b'chr3']
pos = [2, 7, 3, 9, 6, 11]
dp = [35, 12, 78, 22, 99, 96]
qd = [4.5, 6.7, 1.2, 4.4, 2.8, 3.2]
ac = [(1, 2), (3, 4), (5, 6), (7, 8), (9, 10), (11, 12)]

with h5py.File(h5_file, mode='w') as h5f:
    h5g = h5f.create_group('/3L/variants')
    h5g.create_dataset('CHROM', data=chrom, chunks=True)
    h5g.create_dataset('POS', data=pos, chunks=True)
    h5g.create_dataset('DP', data=dp, chunks=True)
    h5g.create_dataset('QD', data=qd, chunks=True)
    h5g.create_dataset('AC', data=ac, chunks=True)


callset = h5py.File(h5_file, mode='r')
variants = allel.VariantChunkedTable(callset['/3L/variants'],
                                names=['CHROM', 'POS', 'AC', 'QD', 'DP'])

因此，variants变量（本例中只有 6 行）如下所示：

>>> variants
<VariantChunkedTable shape=(6,) dtype=[('CHROM', 'S4'), ('POS', '<i8'), ('AC', '<i8', (2,)), ('QD', '<f8'), ('DP', '<i8')] 
nbytes=264 cbytes=264 cratio=1.0 values=h5py._hl.group.Group>
    CHROM   POS AC      QD  DP
0   b'chr1' 2   [1 2]   4.5 35
1   b'chr1' 7   [3 4]   6.7 12
2   b'chr2' 3   [5 6]   1.2 78
3   b'chr2' 9   [7 8]   4.4 22
4   b'chr3' 6   [ 9 10] 2.8 99
5   b'chr3' 11  [11 12] 3.2 96

您已正确定义pos, alt,dp等数组（即pos = variants['POS'][:], 等）

然后，在您的循环中，我假设您的目标是遍历variants变量的前 10 行，从每行获取一些值，例如pos[variant], ac[variant][0]，dp[variant]并在 GraphDatabase 中使用这些属性创建一个新节点。

您当前编写循环的方式是variants，每次迭代都会得到一整行，并尝试将其用作索引来访问pos, alt, ... 数组中的元素，这会引发错误。

正确的做法是遍历数字索引；在我的示例中，要遍历variants您应该运行的变量的所有 6 行：

for i in range(len(variants)):
    print(f"> Row {i}")
    print(pos[i])
    print(dp[i])
    print(ac[i][0])

然后可以将 pos[i]、dp[i] 等值成对输入db.nodes.create。name=value当然，对于前 10 行，您只需要使用for i in range(10).

于 2019-01-11T18:39:46.423 回答