Mercurial > repos > dfornika > micall_lite

comparison amino2consensus.py @ 15:61667ff2c8b5 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/dfornika/galaxytools/tree/master/tools/micall-lite commit 822e7e1c2de31a72c2a13bcc15b9df06b699561f-dirty
author dfornika
date Wed, 25 Sep 2019 17:49:41 -0400
parents
children 43a987c03ec5
comparison
equal deleted inserted replaced
14:79d9866da30f 15:61667ff2c8b5
1 #!/usr/bin/env python
2
3 from __future__ import print_function
4
5 import csv
6 import argparse
7
8 AMINO_ACIDS = ['A','C','D','E','F','G','H','I','K','L','M','N','P','Q','R','S','T','V','W','Y','*']
9
10 def determine_amino(amino_counts, threshold):
11 amino = ""
12 total_count = sum(amino_counts.values())
13 amino_with_max_counts = sorted(amino_counts, key=amino_counts.get, reverse=True)[0]
14 if total_count == 0:
15 amino = "#"
16 elif (amino_counts[amino_with_max_counts] / total_count) > threshold:
17 amino = amino_with_max_counts
18 else:
19 amino = "@"
20 return amino
21
22 def determine_first_region(amino_file):
23 with open(amino_file, newline='') as f:
24 reader = csv.DictReader(f)
25 row = next(reader)
26 region = row['region']
27 return region
28
29 def main(args):
30 current_region = determine_first_region(args.amino)
31 seq = []
32 with open(args.amino, newline='') as f:
33 reader = csv.DictReader(f)
34 for row in reader:
35 if row['region'] == current_region:
36 amino_counts = {}
37 for amino_acid in AMINO_ACIDS:
38 amino_counts[amino_acid] = int(row[amino_acid])
39 amino = determine_amino(amino_counts, args.threshold)
40 seq.append(amino)
41 else:
42 print(">" + current_region)
43 print(''.join(seq))
44 current_region = row['region']
45 seq = []
46 amino_counts = {}
47 for amino_acid in AMINO_ACIDS:
48 amino_counts[amino_acid] = int(row[amino_acid])
49 amino = determine_amino(amino_counts, args.threshold)
50 seq.append(amino)
51 print(">" + current_region)
52 print(''.join(seq))
53
54
55 if __name__ == '__main__':
56 parser = argparse.ArgumentParser()
57 parser.add_argument("amino", help="MiCall amino.csv output file")
58 parser.add_argument("--threshold", default=0.15, help="Threshold for calling")
59 args = parser.parse_args()
60 main(args)