|
0
|
1 #!/usr/bin/env python
|
|
|
2
|
|
|
3 import argparse
|
|
|
4 import gzip
|
|
|
5 import os
|
|
4
|
6 from collections import OrderedDict
|
|
|
7
|
|
0
|
8 import yaml
|
|
|
9 from Bio.SeqIO.QualityIO import FastqGeneralIterator
|
|
|
10
|
|
|
11 OUTPUT_DBKEY_DIR = 'output_dbkey'
|
|
|
12 OUTPUT_METRICS_DIR = 'output_metrics'
|
|
|
13
|
|
|
14
|
|
4
|
15 def get_sample_name(file_path):
|
|
0
|
16 base_file_name = os.path.basename(file_path)
|
|
|
17 if base_file_name.find(".") > 0:
|
|
|
18 # Eliminate the extension.
|
|
|
19 return os.path.splitext(base_file_name)[0]
|
|
|
20 return base_file_name
|
|
|
21
|
|
|
22
|
|
|
23 def get_dbkey(dnaprints_dict, key, s):
|
|
|
24 # dnaprints_dict looks something like this:
|
|
|
25 # {'brucella': {'NC_002945v4': ['11001110', '11011110', '11001100']}
|
|
|
26 # {'bovis': {'NC_006895': ['11111110', '00010010', '01111011']}}
|
|
|
27 d = dnaprints_dict.get(key, {})
|
|
|
28 for data_table_value, v_list in d.items():
|
|
|
29 if s in v_list:
|
|
|
30 return data_table_value
|
|
|
31 return ""
|
|
|
32
|
|
|
33
|
|
|
34 def get_dnaprints_dict(dnaprint_fields):
|
|
|
35 # A dndprint_fields entry looks something liek this.
|
|
|
36 # [['AF2122', '/galaxy/tool-data/vsnp/AF2122/dnaprints/NC_002945v4.yml']]
|
|
|
37 dnaprints_dict = {}
|
|
|
38 for item in dnaprint_fields:
|
|
|
39 # Here item is a 2-element list of data
|
|
|
40 # table components, # value and path.
|
|
|
41 value = item[0]
|
|
|
42 path = item[1].strip()
|
|
|
43 with open(path, "rt") as fh:
|
|
|
44 # The format of all dnaprints yaml
|
|
|
45 # files is something like this:
|
|
|
46 # brucella:
|
|
|
47 # - 0111111111111111
|
|
|
48 print_dict = yaml.load(fh, Loader=yaml.Loader)
|
|
|
49 for print_dict_k, print_dict_v in print_dict.items():
|
|
|
50 dnaprints_v_dict = dnaprints_dict.get(print_dict_k, {})
|
|
|
51 if len(dnaprints_v_dict) > 0:
|
|
|
52 # dnaprints_dict already contains k (e.g., 'brucella',
|
|
|
53 # and dnaprints_v_dict will be a dictionary # that
|
|
|
54 # looks something like this:
|
|
|
55 # {'NC_002945v4': ['11001110', '11011110', '11001100']}
|
|
|
56 value_list = dnaprints_v_dict.get(value, [])
|
|
|
57 value_list = value_list + print_dict_v
|
|
|
58 dnaprints_v_dict[value] = value_list
|
|
|
59 else:
|
|
|
60 # dnaprints_v_dict is an empty dictionary.
|
|
|
61 dnaprints_v_dict[value] = print_dict_v
|
|
|
62 dnaprints_dict[print_dict_k] = dnaprints_v_dict
|
|
|
63 # dnaprints_dict looks something like this:
|
|
|
64 # {'brucella': {'NC_002945v4': ['11001110', '11011110', '11001100']}
|
|
|
65 # {'bovis': {'NC_006895': ['11111110', '00010010', '01111011']}}
|
|
|
66 return dnaprints_dict
|
|
|
67
|
|
|
68
|
|
|
69 def get_group_and_dbkey(dnaprints_dict, brucella_string, brucella_sum, bovis_string, bovis_sum, para_string, para_sum):
|
|
|
70 if brucella_sum > 3:
|
|
|
71 group = "Brucella"
|
|
|
72 dbkey = get_dbkey(dnaprints_dict, "brucella", brucella_string)
|
|
|
73 elif bovis_sum > 3:
|
|
|
74 group = "TB"
|
|
|
75 dbkey = get_dbkey(dnaprints_dict, "bovis", bovis_string)
|
|
|
76 elif para_sum >= 1:
|
|
|
77 group = "paraTB"
|
|
|
78 dbkey = get_dbkey(dnaprints_dict, "para", para_string)
|
|
|
79 else:
|
|
|
80 group = ""
|
|
|
81 dbkey = ""
|
|
|
82 return group, dbkey
|
|
|
83
|
|
|
84
|
|
|
85 def get_oligo_dict():
|
|
|
86 oligo_dict = {}
|
|
|
87 oligo_dict["01_ab1"] = "AATTGTCGGATAGCCTGGCGATAACGACGC"
|
|
|
88 oligo_dict["02_ab3"] = "CACACGCGGGCCGGAACTGCCGCAAATGAC"
|
|
|
89 oligo_dict["03_ab5"] = "GCTGAAGCGGCAGACCGGCAGAACGAATAT"
|
|
|
90 oligo_dict["04_mel"] = "TGTCGCGCGTCAAGCGGCGTGAAATCTCTG"
|
|
|
91 oligo_dict["05_suis1"] = "TGCGTTGCCGTGAAGCTTAATTCGGCTGAT"
|
|
|
92 oligo_dict["06_suis2"] = "GGCAATCATGCGCAGGGCTTTGCATTCGTC"
|
|
|
93 oligo_dict["07_suis3"] = "CAAGGCAGATGCACATAATCCGGCGACCCG"
|
|
|
94 oligo_dict["08_ceti1"] = "GTGAATATAGGGTGAATTGATCTTCAGCCG"
|
|
|
95 oligo_dict["09_ceti2"] = "TTACAAGCAGGCCTATGAGCGCGGCGTGAA"
|
|
|
96 oligo_dict["10_canis4"] = "CTGCTACATAAAGCACCCGGCGACCGAGTT"
|
|
|
97 oligo_dict["11_canis"] = "ATCGTTTTGCGGCATATCGCTGACCACAGC"
|
|
|
98 oligo_dict["12_ovis"] = "CACTCAATCTTCTCTACGGGCGTGGTATCC"
|
|
|
99 oligo_dict["13_ether2"] = "CGAAATCGTGGTGAAGGACGGGACCGAACC"
|
|
|
100 oligo_dict["14_63B1"] = "CCTGTTTAAAAGAATCGTCGGAACCGCTCT"
|
|
|
101 oligo_dict["15_16M0"] = "TCCCGCCGCCATGCCGCCGAAAGTCGCCGT"
|
|
|
102 oligo_dict["16_mel1b"] = "TCTGTCCAAACCCCGTGACCGAACAATAGA"
|
|
|
103 oligo_dict["17_tb157"] = "CTCTTCGTATACCGTTCCGTCGTCACCATGGTCCT"
|
|
|
104 oligo_dict["18_tb7"] = "TCACGCAGCCAACGATATTCGTGTACCGCGACGGT"
|
|
|
105 oligo_dict["19_tbbov"] = "CTGGGCGACCCGGCCGACCTGCACACCGCGCATCA"
|
|
|
106 oligo_dict["20_tb5"] = "CCGTGGTGGCGTATCGGGCCCCTGGATCGCGCCCT"
|
|
|
107 oligo_dict["21_tb2"] = "ATGTCTGCGTAAAGAAGTTCCATGTCCGGGAAGTA"
|
|
|
108 oligo_dict["22_tb3"] = "GAAGACCTTGATGCCGATCTGGGTGTCGATCTTGA"
|
|
|
109 oligo_dict["23_tb4"] = "CGGTGTTGAAGGGTCCCCCGTTCCAGAAGCCGGTG"
|
|
|
110 oligo_dict["24_tb6"] = "ACGGTGATTCGGGTGGTCGACACCGATGGTTCAGA"
|
|
|
111 oligo_dict["25_para"] = "CCTTTCTTGAAGGGTGTTCG"
|
|
|
112 oligo_dict["26_para_sheep"] = "CGTGGTGGCGACGGCGGCGGGCCTGTCTAT"
|
|
|
113 oligo_dict["27_para_cattle"] = "TCTCCTCGGTCGGTGATTCGGGGGCGCGGT"
|
|
|
114 return oligo_dict
|
|
|
115
|
|
|
116
|
|
|
117 def get_seq_counts(value, fastq_list, gzipped):
|
|
|
118 count = 0
|
|
|
119 for fastq_file in fastq_list:
|
|
4
|
120 if gzipped:
|
|
0
|
121 with gzip.open(fastq_file, 'rt') as fh:
|
|
|
122 for title, seq, qual in FastqGeneralIterator(fh):
|
|
|
123 count += seq.count(value)
|
|
|
124 else:
|
|
|
125 with open(fastq_file, 'r') as fh:
|
|
|
126 for title, seq, qual in FastqGeneralIterator(fh):
|
|
|
127 count += seq.count(value)
|
|
|
128 return(value, count)
|
|
|
129
|
|
|
130
|
|
|
131 def get_species_counts(fastq_list, gzipped):
|
|
|
132 count_summary = {}
|
|
|
133 oligo_dict = get_oligo_dict()
|
|
|
134 for v1 in oligo_dict.values():
|
|
|
135 returned_value, count = get_seq_counts(v1, fastq_list, gzipped)
|
|
|
136 for key, v2 in oligo_dict.items():
|
|
|
137 if returned_value == v2:
|
|
|
138 count_summary.update({key: count})
|
|
|
139 count_list = []
|
|
|
140 for v in count_summary.values():
|
|
|
141 count_list.append(v)
|
|
|
142 brucella_sum = sum(count_list[:16])
|
|
|
143 bovis_sum = sum(count_list[16:24])
|
|
|
144 para_sum = sum(count_list[24:])
|
|
|
145 return count_summary, count_list, brucella_sum, bovis_sum, para_sum
|
|
|
146
|
|
|
147
|
|
|
148 def get_species_strings(count_summary):
|
|
|
149 binary_dictionary = {}
|
|
|
150 for k, v in count_summary.items():
|
|
|
151 if v > 1:
|
|
|
152 binary_dictionary.update({k: 1})
|
|
|
153 else:
|
|
|
154 binary_dictionary.update({k: 0})
|
|
|
155 binary_dictionary = OrderedDict(sorted(binary_dictionary.items()))
|
|
|
156 binary_list = []
|
|
|
157 for v in binary_dictionary.values():
|
|
|
158 binary_list.append(v)
|
|
|
159 brucella_binary = binary_list[:16]
|
|
|
160 brucella_string = ''.join(str(e) for e in brucella_binary)
|
|
|
161 bovis_binary = binary_list[16:24]
|
|
|
162 bovis_string = ''.join(str(e) for e in bovis_binary)
|
|
|
163 para_binary = binary_list[24:]
|
|
|
164 para_string = ''.join(str(e) for e in para_binary)
|
|
|
165 return brucella_string, bovis_string, para_string
|
|
|
166
|
|
|
167
|
|
4
|
168 def output_dbkey(file_name, dbkey, output_file):
|
|
0
|
169 # Output the dbkey.
|
|
|
170 with open(output_file, "w") as fh:
|
|
|
171 fh.write("%s" % dbkey)
|
|
|
172
|
|
|
173
|
|
4
|
174 def output_files(fastq_file, count_list, group, dbkey, dbkey_file, metrics_file):
|
|
|
175 base_file_name = get_sample_name(fastq_file)
|
|
0
|
176 output_dbkey(base_file_name, dbkey, dbkey_file)
|
|
|
177 output_metrics(base_file_name, count_list, group, dbkey, metrics_file)
|
|
|
178
|
|
|
179
|
|
4
|
180 def output_metrics(file_name, count_list, group, dbkey, output_file):
|
|
0
|
181 # Output the metrics.
|
|
|
182 with open(output_file, "w") as fh:
|
|
|
183 fh.write("Sample: %s\n" % file_name)
|
|
|
184 fh.write("Brucella counts: ")
|
|
|
185 for i in count_list[:16]:
|
|
|
186 fh.write("%d," % i)
|
|
|
187 fh.write("\nTB counts: ")
|
|
|
188 for i in count_list[16:24]:
|
|
|
189 fh.write("%d," % i)
|
|
|
190 fh.write("\nPara counts: ")
|
|
|
191 for i in count_list[24:]:
|
|
|
192 fh.write("%d," % i)
|
|
|
193 fh.write("\nGroup: %s" % group)
|
|
|
194 fh.write("\ndbkey: %s\n" % dbkey)
|
|
|
195
|
|
|
196
|
|
|
197 if __name__ == '__main__':
|
|
|
198 parser = argparse.ArgumentParser()
|
|
|
199
|
|
|
200 parser.add_argument('--dnaprint_fields', action='append', dest='dnaprint_fields', nargs=2, help="List of dnaprints data table value, name and path fields")
|
|
4
|
201 parser.add_argument('--read1', action='store', dest='read1', help='Required: single read')
|
|
0
|
202 parser.add_argument('--read2', action='store', dest='read2', required=False, default=None, help='Optional: paired read')
|
|
4
|
203 parser.add_argument('--gzipped', action='store_true', dest='gzipped', help='Input files are gzipped')
|
|
|
204 parser.add_argument('--output_dbkey', action='store', dest='output_dbkey', help='Output reference file')
|
|
|
205 parser.add_argument('--output_metrics', action='store', dest='output_metrics', help='Output metrics file')
|
|
0
|
206
|
|
|
207 args = parser.parse_args()
|
|
|
208
|
|
4
|
209 fastq_list = [args.read1]
|
|
|
210 if args.read2 is not None:
|
|
|
211 fastq_list.append(args.read2)
|
|
0
|
212
|
|
|
213 # The value of dnaprint_fields is a list of lists, where each list is
|
|
|
214 # the [value, name, path] components of the vsnp_dnaprints data table.
|
|
|
215 # The data_manager_vsnp_dnaprints tool assigns the dbkey column from the
|
|
|
216 # all_fasta data table to the value column in the vsnp_dnaprints data
|
|
|
217 # table to ensure a proper mapping for discovering the dbkey.
|
|
|
218 dnaprints_dict = get_dnaprints_dict(args.dnaprint_fields)
|
|
|
219
|
|
4
|
220 # Here fastq_list consists of either a single read
|
|
|
221 # or a set of paired reads, producing single outputs.
|
|
|
222 count_summary, count_list, brucella_sum, bovis_sum, para_sum = get_species_counts(fastq_list, args.gzipped)
|
|
|
223 brucella_string, bovis_string, para_string = get_species_strings(count_summary)
|
|
|
224 group, dbkey = get_group_and_dbkey(dnaprints_dict, brucella_string, brucella_sum, bovis_string, bovis_sum, para_string, para_sum)
|
|
|
225 output_files(args.read1, count_list, group, dbkey, dbkey_file=args.output_dbkey, metrics_file=args.output_metrics)
|
