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
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
|
#!/usr/bin/python3
import optparse
import sys
import os
import gzip
import re
#
# Global configuration variables.
#
path = ""
aln_path = ""
out_path = ""
batch_id = -1
def parse_command_line():
global path
global batch_id
p = optparse.OptionParser()
#
# Add options.
#
p.add_option("-p", action="store", dest="path",
help="path to Stacks directory.")
p.add_option("-b", action="store", dest="batch_id",
help="Stacks batch ID.")
#
# Parse the command line
#
(opts, args) = p.parse_args()
if opts.path != None:
path = opts.path
if opts.batch_id != None:
batch_id = int(opts.batch_id)
if len(path) == 0 or os.path.exists(path) == False:
print("You must specify a valid path to Stacks input files.", file=sys.stderr)
p.print_help()
sys.exit()
if batch_id < 0:
pritn >> sys.stderr, "You must specify the batch ID that was supplied to Stacks."
p.print_help()
sys.exit()
if path.endswith("/") == False:
path += "/"
def find_stacks_files(path, files):
try:
entries = os.listdir(path)
for entry in entries:
pos = entry.find(".matches.tsv.gz")
if (pos == -1):
pos = entry.find(".matches.tsv")
if (pos != -1):
files.append(entry[0:pos])
print("Found", len(files), "Stacks samples.", file=sys.stderr)
except:
print("Unable to read files from Stacks directory, '" + path + "'", file=sys.stderr)
def parse_cigar(cigar, components):
#
# Parse a cigar string, e.g. 48M1D47M or 43M3D52M3D.
#
start = 0
end = 0
dist = ""
for c in cigar:
if c.isalpha():
dist = int(cigar[start:end])
op = c.upper();
end += 1
start = end
components.append((op, dist))
else:
end += 1
def adjust_snps(cigar, sample_snps):
#
# Adjust SNP positions according to how this sample was aligned to the catalog
# (which is described by the CIGAR string).
#
if len(sample_snps) == 0:
return
comp = []
parse_cigar(cigar, comp)
index = 0
offset = 0
bp = 0
for (op, dist) in comp:
if op == 'M' or op == 'I':
bp += dist
while index < len(sample_snps) and sample_snps[index] < bp:
sample_snps[index] += offset
index += 1
if op == 'D':
offset += dist
def count_sample_snps(path, file, sample_snps):
matches = {}
cigars = {}
#
# Open the matches file and load the matches to the catalog.
#
p = path + file + ".matches.tsv.gz"
if os.path.exists(p):
gzipped = True;
fh = gzip.open(p, 'rb')
else:
gzipped = False;
fh = open(path + file + ".matches.tsv", "r")
for line in fh:
line = line.strip("\n")
if len(line) == 0 or line[0] == "#":
continue
parts = line.split("\t")
cat_locus = int(parts[2])
sample_locus = int(parts[4])
if len(parts) == 9:
cigar = parts[8] if len(parts[8]) > 0 else ""
else:
cigar = ""
if sample_locus not in matches:
matches[sample_locus] = cat_locus
else:
if cat_locus != matches[sample_locus]:
print("Error: sample locus", sample_locus, "matches more than one catalog locus.", file=sys.stderr)
if len(cigar) > 0:
if sample_locus not in cigars:
cigars[sample_locus] = cigar
else:
if cigar != cigars[sample_locus]:
print("Error: sample locus", sample_locus, "has multiple cigar alignments.", file=sys.stderr)
fh.close()
#
# Open the SNPs file and record all the SNP positions found in this sample.
#
if gzipped:
fh = gzip.open(path + file + ".snps.tsv.gz", "rb")
else:
fh = open(path + file + ".snps.tsv", "r")
snps = {}
for line in fh:
if line[0] == "#":
continue
if len(line) == 0:
continue
parts = line.split("\t")
sample_locus = int(parts[2])
col = int(parts[3])
model = parts[4]
if model != "E":
continue
if sample_locus not in matches:
continue
if sample_locus not in snps:
snps[sample_locus] = []
snps[sample_locus].append(col)
fh.close()
#
# Adjust SNP positions according to the gapped alignments recorded in the CIGAR string.
#
for sample_locus in snps:
if sample_locus in cigars:
adjust_snps(cigars[sample_locus], snps[sample_locus])
snp_cnt = 0
#
# Transfer this sample's SNPs to the catalog level dictionary.
#
for sample_locus in snps:
cat_locus = matches[sample_locus]
for col in snps[sample_locus]:
if cat_locus in sample_snps:
if col in sample_snps[cat_locus]:
sample_snps[cat_locus][col] += 1
# print >> sys.stderr, "CatLocus:", cat_locus, "; col:", col
else:
sample_snps[cat_locus][col] = 1
snp_cnt += 1
else:
sample_snps[cat_locus] = {}
sample_snps[cat_locus][col] = 1
snp_cnt += 1
return snp_cnt
def count_catalog_snps(path, batch_id, catalog_snps):
#
# Open the tags file and rewrite it with the alignment coordinates.
#
p = path + "batch_" + str(batch_id) + ".catalog.snps.tsv.gz"
if os.path.exists(p):
gzipped = True;
fh = gzip.open(path + "batch_" + str(batch_id) + ".catalog.snps.tsv.gz", "rb")
else:
gzipped = False;
fh = open(path + "batch_" + str(batch_id) + ".catalog.snps.tsv", "r")
snp_cnt = 0
for line in fh:
if line[0] == "#":
continue
if len(line) == 0:
continue
parts = line.split("\t")
cat_locus = int(parts[2])
col = int(parts[3])
model = parts[4]
if model != "E":
continue
snp_cnt += 1
if cat_locus not in catalog_snps:
catalog_snps[cat_locus] = []
catalog_snps[cat_locus].append(col)
fh.close()
return snp_cnt
# # # #
# # ------------------------------------------------------------------------------------------- # #
# # # #
parse_command_line()
files = []
catalog_snps = {}
sample_snps = {}
#
# Find all the sample files by looking for the matches files in the Stacks directory.
#
find_stacks_files(path, files)
#
# Count all the SNPs identified in the samples.
#
i = 1
for file in files:
print("Processing file", str(i), "of", len(files), "['" + file + "']", file=sys.stderr)
cnt = count_sample_snps(path, file, sample_snps)
print(" Found", cnt, "heterozygous SNPs in sample.", file=sys.stderr)
i += 1
total_snps = 0
for locus in sample_snps:
for col in sample_snps[locus]:
total_snps += 1
print("Found", total_snps, "variable sites across the population.", file=sys.stderr)
#
# Count all the SNPs found in the catalog.
#
print("Processing the catalog", file=sys.stderr)
cnt = count_catalog_snps(path, batch_id, catalog_snps)
print(" Found", cnt, "heterozygous SNPs in the catalog.", file=sys.stderr)
#
# Count all the SNPs in the catalog but not in any sample: these are the fixed differences cstacks identified.
#
total_snps = 0
fixed_snps = 0
for locus in catalog_snps:
if locus not in sample_snps:
continue
c = {}
for col in catalog_snps[locus]:
c[col] = 1
total_snps += 1
if col not in sample_snps[locus]:
# print >> sys.stderr, "Locus:", locus, "Catalog SNPs:", catalog_snps[locus], "Sample SNPs:", sample_snps[locus]
fixed_snps += 1
for col in sample_snps[locus]:
if col not in c:
print("Locus:", locus, "col:", col, "Catalog SNPs:", catalog_snps[locus], "Sample SNPs:", sample_snps[locus])
print("Found", total_snps, "SNPs across all samples and in the catalog.", file=sys.stderr)
print("Found", fixed_snps, "fixed SNPs only in the catalog.", file=sys.stderr)
|
