1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
|
import lmdb
import json
import click
from dataclasses import dataclass
from contextlib import contextmanager
import numpy as np
@dataclass
class GenotypeMatrix:
"""Store the actual Genotype Matrix"""
matrix: np.ndarray
transpose: np.ndarray
metadata: dict
@contextmanager
def save_excursion(stream):
"""Context manager to restore stream position after use."""
position = stream.tell()
try:
yield
finally:
stream.seek(position)
def count_lines(stream) -> int:
"""Count the number of lines in the stream from the current
position."""
count = 0
with save_excursion(stream):
while stream.readline().strip():
count += 1
return count
def create_database(db_path: str) -> lmdb.Environment:
"""Create or open an LMDB environment."""
return lmdb.open(db_path, map_size=100 * 1024 * 1024, create=True)
def genotype_db_put(db: lmdb.Environment, genotype: GenotypeMatrix) -> bool:
metadata = json.dumps(genotype.metadata).encode("utf-8")
with db.begin(write=True) as txn:
txn.put(b"matrix", genotype.matrix.tobytes())
txn.put(b"transpose", genotype.transpose.tobytes())
txn.put(b"metadata", metadata)
return True
def read_geno_file(genotype_file: str) -> GenotypeMatrix:
"""Read a geno file and return a GenotypeMatrix object."""
with open(genotype_file, 'r') as stream:
# Read file metadata
file_metadata = {}
while True:
line = stream.readline().strip()
if not line:
break
if line.startswith('#'):
continue
if line.startswith('@'):
key, value = line[1:].split(':', 1)
file_metadata[key] = value
else:
stream.seek(stream.tell() - len(line) - 1)
break
# Read header
header = stream.readline().strip().split()
metadata_columns = ["Chr", "Locus", "cM",
"Mb"] if "Mb" in header else ["Chr", "Locus", "cM"]
individuals = header[len(metadata_columns):]
# Read data
nrows = count_lines(stream)
ncols = len(individuals)
matrix = np.zeros((nrows, ncols), dtype=np.uint8)
maternal = file_metadata.get("mat")
paternal = file_metadata.get("pat")
heterozygous = file_metadata.get("het")
unknown = file_metadata.get("unk")
metadata = {
"nrows": nrows,
"ncols": ncols,
"individuals": individuals,
"metadata_keys": metadata_columns + ["individuals"]
}
for key in metadata_columns[2:]:
metadata[key] = []
locus, chromosomes = [], []
for i in range(nrows):
line = stream.readline().strip().split()
meta, data = line[:len(metadata_columns)
], line[len(metadata_columns):]
for j, element in enumerate(data):
if element.isdigit():
matrix[i, j] = int(element)
elif element == maternal:
matrix[i, j] = 0
elif element == paternal:
matrix[i, j] = 1
elif element == heterozygous:
matrix[i, j] = 2
elif element == unknown:
matrix[i, j] = 3
data = dict(zip(metadata_columns, meta))
locus.append(data.get("Locus"))
chromosomes.append(data.get("Chr"))
for col in metadata_columns[2:]:
metadata[col].append(float(data.get(col)))
metadata["locus"] = locus
metadata["chromosomes"] = chromosomes
genotype_matrix = GenotypeMatrix(matrix, matrix.T, metadata)
return genotype_matrix
def genotype_db_get(db: lmdb.Environment) -> GenotypeMatrix:
with db.begin() as txn:
metadata = json.loads(txn.get(b"metadata").decode("utf-8"))
nrows, ncols = metadata.get("nrows"), metadata.get("ncols")
matrix, transpose = txn.get(b"matrix"), txn.get(b"transpose")
# Double check this
matrix = np.frombuffer(matrix, dtype=np.uint8).reshape(nrows, ncols)
transpose = np.frombuffer(transpose, dtype=np.uint8)\
.reshape(ncols, nrows)
return GenotypeMatrix(matrix, transpose, metadata)
@click.command(help="Import the genotype file")
@click.argument("geno_file")
@click.argument("genotype_database")
def import_genotype(geno_file: str, genotype_database: str):
with create_database(genotype_database) as db:
genotype_db_put(db, read_geno_file(geno_file))
@click.command(help="Print the current matrix")
@click.argument("database_directory")
def print_current_matrix(database_directory: str):
"""Print the current matrix in the database."""
with create_database(database_directory) as db:
current = genotype_db_get(db)
print(f"Matrix: {current.matrix}")
print(f"Transpose: {current.transpose}")
print(f"Metadata: {current.metadata}")
@click.group()
def cli():
pass
cli.add_command(print_current_matrix)
cli.add_command(import_genotype)
if __name__ == "__main__":
cli()
|
