add comments to clade script
@@ -41,8 +41,10 @@ if __name__ == '__main__':
output = open(args.output, 'w')
print('name\tclade\tparent clades', file=output)
+ # break the sequences into chunks, align each to the reference, and assign clades one-by-one
done = False
while not done:
+ # generate a chunk with chunk-size sequences
chunk = []
while len(chunk)<args.chunk_size and (not done):
try:
@@ -55,6 +57,7 @@ if __name__ == '__main__':
with open(in_fname, 'wt') as fh:
SeqIO.write(chunk + [ref], fh, 'fasta')
+ # align using augur command structure
cmd = generate_alignment_cmd('mafft', args.nthreads, None, in_fname, out_fname, log_fname)
success = run_shell_command(cmd)
@@ -65,6 +68,7 @@ if __name__ == '__main__':
for seq in alignment:
if seq.id==ref.id:
continue
+ # read sequence and all its annotated features
seq_container = tmpNode()
seq_str = str(seq.seq)
seq_container.sequences['nuc'] = {i:c for i,c in enumerate(seq_str)}
@@ -72,11 +76,16 @@ if __name__ == '__main__':
if feat.type != 'source':
seq_container.sequences[fname] = {i:c for i,c in enumerate(safe_translate(feat.extract(seq_str)))}
+ # for each clade, check whether it matches any of the clade definitions in the tsv
matches = []
for clade_name, clade_alleles in clade_designations.items():
if is_node_in_clade(clade_alleles, seq_container, ref):
matches.append(clade_name)
+ # print the last match as clade assignment and all others as ancestral clades
+ # note that this assumes that clades in the tsv are ordered by order of appearence.
+ # furthermore, this will only work if parent clades don't have definitions that exclude
+ # child clades, i.e. positions can only be additive for this to work.
if matches:
print(f"{seq.description}\t{matches[-1]}\t{', '.join(matches[:-1])}", file=output)
else:
