|
kkonganti@11
|
1 #!/usr/bin/env python3
|
|
kkonganti@11
|
2
|
|
kkonganti@11
|
3 # Kranti Konganti
|
|
kkonganti@11
|
4
|
|
kkonganti@11
|
5 import argparse
|
|
kkonganti@11
|
6 import glob
|
|
kkonganti@11
|
7 import inspect
|
|
kkonganti@11
|
8 import json
|
|
kkonganti@11
|
9 import logging
|
|
kkonganti@11
|
10 import os
|
|
kkonganti@11
|
11 import pickle
|
|
kkonganti@11
|
12 import pprint
|
|
kkonganti@11
|
13 import re
|
|
kkonganti@11
|
14 from collections import defaultdict
|
|
kkonganti@11
|
15
|
|
kkonganti@11
|
16 import yaml
|
|
kkonganti@11
|
17
|
|
kkonganti@11
|
18
|
|
kkonganti@11
|
19 # Multiple inheritence for pretty printing of help text.
|
|
kkonganti@11
|
20 class MultiArgFormatClasses(argparse.RawTextHelpFormatter, argparse.ArgumentDefaultsHelpFormatter):
|
|
kkonganti@11
|
21 pass
|
|
kkonganti@11
|
22
|
|
kkonganti@11
|
23
|
|
kkonganti@11
|
24 # Main
|
|
kkonganti@11
|
25 def main() -> None:
|
|
kkonganti@11
|
26 """
|
|
kkonganti@11
|
27 The succesful execution of this script requires access to bettercallsal formatted
|
|
kkonganti@11
|
28 db flat files. On raven2, they are at /hpc/db/bettercallsall/PDGXXXXXXXXXX.XXXXX
|
|
kkonganti@11
|
29
|
|
kkonganti@11
|
30 It takes the ACC2SERO.pickle file and *.reference_target.cluster_list.tsv file
|
|
kkonganti@11
|
31 for that particular NCBI Pathogens release from the db directory mentioned with
|
|
kkonganti@11
|
32 -db option and a root parent directory of the `salmon quant` results mentioned
|
|
kkonganti@11
|
33 with -sal option and generates a final results table with number of reads
|
|
kkonganti@11
|
34 mapped and a .json file to be used with MultiQC to generate a stacked bar plot.
|
|
kkonganti@11
|
35
|
|
kkonganti@11
|
36 Using -url option optionally adds an extra column of NCBI Pathogens Isolates
|
|
kkonganti@11
|
37 Browser, which directly links out to NCBI Pathogens Isolates SNP viewer tool.
|
|
kkonganti@11
|
38 """
|
|
kkonganti@11
|
39 # Set logging.
|
|
kkonganti@11
|
40 logging.basicConfig(
|
|
kkonganti@11
|
41 format="\n" + "=" * 55 + "\n%(asctime)s - %(levelname)s\n" + "=" * 55 + "\n%(message)s\n\n",
|
|
kkonganti@11
|
42 level=logging.DEBUG,
|
|
kkonganti@11
|
43 )
|
|
kkonganti@11
|
44
|
|
kkonganti@11
|
45 # Debug print.
|
|
kkonganti@11
|
46 ppp = pprint.PrettyPrinter(width=55)
|
|
kkonganti@11
|
47 prog_name = inspect.stack()[0].filename
|
|
kkonganti@11
|
48
|
|
kkonganti@11
|
49 parser = argparse.ArgumentParser(
|
|
kkonganti@11
|
50 prog=prog_name, description=main.__doc__, formatter_class=MultiArgFormatClasses
|
|
kkonganti@11
|
51 )
|
|
kkonganti@11
|
52
|
|
kkonganti@11
|
53 required = parser.add_argument_group("required arguments")
|
|
kkonganti@11
|
54
|
|
kkonganti@11
|
55 required.add_argument(
|
|
kkonganti@11
|
56 "-sal",
|
|
kkonganti@11
|
57 dest="salmon_res_dir",
|
|
kkonganti@11
|
58 default=False,
|
|
kkonganti@11
|
59 required=True,
|
|
kkonganti@11
|
60 help="Absolute UNIX path to the parent directory that contains the\n"
|
|
kkonganti@11
|
61 + "`salmon quant` results directory. For example, if path to\n"
|
|
kkonganti@11
|
62 + "`quant.sf` is in /hpc/john_doe/test/salmon_res/quant.sf, then\n"
|
|
kkonganti@11
|
63 + "use this command-line option as:\n"
|
|
kkonganti@11
|
64 + "-sal /hpc/john_doe/test",
|
|
kkonganti@11
|
65 )
|
|
kkonganti@11
|
66 required.add_argument(
|
|
kkonganti@11
|
67 "-snp",
|
|
kkonganti@11
|
68 dest="rtc",
|
|
kkonganti@11
|
69 default=False,
|
|
kkonganti@11
|
70 required=True,
|
|
kkonganti@11
|
71 help="Absolute UNIX Path to the PDG SNP reference target cluster\n"
|
|
kkonganti@11
|
72 + "metadata file. On raven2, these are located at\n"
|
|
kkonganti@11
|
73 + "/hpc/db/bettercallsal/PDGXXXXXXXXXX.XXXXX\n"
|
|
kkonganti@11
|
74 + "Required if -sal is on.",
|
|
kkonganti@11
|
75 )
|
|
kkonganti@11
|
76 required.add_argument(
|
|
kkonganti@11
|
77 "-pickle",
|
|
kkonganti@11
|
78 dest="acc2sero",
|
|
kkonganti@11
|
79 default=False,
|
|
kkonganti@11
|
80 required=True,
|
|
kkonganti@11
|
81 help="Absolute UNIX Path to the *ACC2SERO.pickle\n"
|
|
kkonganti@11
|
82 + "metadata file. On raven2, these are located at\n"
|
|
kkonganti@11
|
83 + "/hpc/db/bettercallsal/PDGXXXXXXXXXX.XXXXX\n"
|
|
kkonganti@11
|
84 + "Required if -sal is on.",
|
|
kkonganti@11
|
85 )
|
|
kkonganti@11
|
86 parser.add_argument(
|
|
kkonganti@11
|
87 "-op",
|
|
kkonganti@11
|
88 dest="out_prefix",
|
|
kkonganti@11
|
89 default="bettercallsal.tblsum",
|
|
kkonganti@11
|
90 required=False,
|
|
kkonganti@11
|
91 help="Set the output file(s) prefix for output(s) generated\n" + "by this program.",
|
|
kkonganti@11
|
92 )
|
|
kkonganti@11
|
93 parser.add_argument(
|
|
kkonganti@11
|
94 "-url",
|
|
kkonganti@11
|
95 dest="show_snp_clust_info",
|
|
kkonganti@11
|
96 default=False,
|
|
kkonganti@11
|
97 required=False,
|
|
kkonganti@11
|
98 action="store_true",
|
|
kkonganti@11
|
99 help="Show SNP cluster participation information of the final genome hit.\n"
|
|
kkonganti@11
|
100 + "This may be useful to see a relative placement of your sample in\n"
|
|
kkonganti@11
|
101 + "NCBI Isolates SNP Tree Viewer based on genome similarity but however\n"
|
|
kkonganti@11
|
102 + "due to rapid nature of the updates at NCBI Pathogen Detection Project,\n"
|
|
kkonganti@11
|
103 + "the placement may be in an outdated cluster.",
|
|
kkonganti@11
|
104 )
|
|
kkonganti@11
|
105
|
|
kkonganti@11
|
106 args = parser.parse_args()
|
|
kkonganti@11
|
107 salmon_res_dir = args.salmon_res_dir
|
|
kkonganti@11
|
108 out_prefix = args.out_prefix
|
|
kkonganti@11
|
109 show_snp_clust_col = args.show_snp_clust_info
|
|
kkonganti@11
|
110 rtc = args.rtc
|
|
kkonganti@11
|
111 pickled_sero = args.acc2sero
|
|
kkonganti@11
|
112 no_hit = "No genome hit"
|
|
kkonganti@11
|
113 no_presence = "Salmonella presence not detected"
|
|
kkonganti@11
|
114 bcs_sal_yn_prefix = "bettercallsal_salyn"
|
|
kkonganti@11
|
115 sal_y = "Detected"
|
|
kkonganti@11
|
116 sal_n = "Not detected"
|
|
kkonganti@11
|
117 null_value = "NULL"
|
|
kkonganti@11
|
118 assm_pat = re.compile(r"GC[AF]\_[0-9]+\.*[0-9]*")
|
|
kkonganti@11
|
119 ncbi_pathogens_base_url = "https://www.ncbi.nlm.nih.gov/pathogens/"
|
|
kkonganti@11
|
120 ncbi_pathogens_genome_base = "https://www.ncbi.nlm.nih.gov/datasets/genome/"
|
|
kkonganti@11
|
121
|
|
kkonganti@11
|
122 sample2salmon, snp_clusters, multiqc_salmon_counts, seen_sero, sal_yn = (
|
|
kkonganti@11
|
123 defaultdict(defaultdict),
|
|
kkonganti@11
|
124 defaultdict(defaultdict),
|
|
kkonganti@11
|
125 defaultdict(defaultdict),
|
|
kkonganti@11
|
126 defaultdict(int),
|
|
kkonganti@11
|
127 defaultdict(int),
|
|
kkonganti@11
|
128 )
|
|
kkonganti@11
|
129
|
|
kkonganti@11
|
130 cell_colors_yml = {
|
|
kkonganti@11
|
131 bcs_sal_yn_prefix: {sal_y: "#c8e6c9 !important;", sal_n: "#ffcdd2 !important;"}
|
|
kkonganti@11
|
132 }
|
|
kkonganti@11
|
133
|
|
kkonganti@11
|
134 salmon_comb_res = os.path.join(os.getcwd(), out_prefix + ".txt")
|
|
kkonganti@11
|
135 bcs_sal_yn = re.sub(out_prefix, bcs_sal_yn_prefix + ".tblsum", salmon_comb_res)
|
|
kkonganti@11
|
136 cell_colors_yml_file = re.sub(
|
|
kkonganti@11
|
137 out_prefix + ".txt", bcs_sal_yn_prefix + ".cellcolors.yml", salmon_comb_res
|
|
kkonganti@11
|
138 )
|
|
kkonganti@11
|
139 salmon_comb_res_mqc = os.path.join(os.getcwd(), str(out_prefix).split(".")[0] + "_mqc.json")
|
|
kkonganti@11
|
140 salmon_res_files = glob.glob(os.path.join(salmon_res_dir, "*", "quant.sf"), recursive=True)
|
|
kkonganti@11
|
141 salmon_res_file_failed = glob.glob(os.path.join(salmon_res_dir, "BCS_NO_CALLS.txt"))
|
|
kkonganti@11
|
142
|
|
kkonganti@11
|
143 if rtc and (not os.path.exists(rtc) or not os.path.getsize(rtc) > 0):
|
|
kkonganti@11
|
144 logging.error(
|
|
kkonganti@11
|
145 "The reference target cluster metadata file,\n"
|
|
kkonganti@11
|
146 + f"{os.path.basename(rtc)} does not exist or is empty!"
|
|
kkonganti@11
|
147 )
|
|
kkonganti@11
|
148 exit(1)
|
|
kkonganti@11
|
149
|
|
kkonganti@11
|
150 if rtc and (not salmon_res_dir or not pickled_sero):
|
|
kkonganti@11
|
151 logging.error("When -rtc is on, -sal and -ps are also required.")
|
|
kkonganti@11
|
152 exit(1)
|
|
kkonganti@11
|
153
|
|
kkonganti@11
|
154 if pickled_sero and (not os.path.exists(pickled_sero) or not os.path.getsize(pickled_sero)):
|
|
kkonganti@11
|
155 logging.error(
|
|
kkonganti@11
|
156 "The pickle file,\n" + f"{os.path.basename(pickled_sero)} does not exist or is empty!"
|
|
kkonganti@11
|
157 )
|
|
kkonganti@11
|
158 exit(1)
|
|
kkonganti@11
|
159
|
|
kkonganti@11
|
160 if salmon_res_dir:
|
|
kkonganti@11
|
161 if not os.path.isdir(salmon_res_dir):
|
|
kkonganti@11
|
162 logging.error("UNIX path\n" + f"{salmon_res_dir}\n" + "does not exist!")
|
|
kkonganti@11
|
163 exit(1)
|
|
kkonganti@11
|
164 if len(salmon_res_files) <= 0:
|
|
kkonganti@11
|
165 # logging.error(
|
|
kkonganti@11
|
166 # "Parent directory,\n"
|
|
kkonganti@11
|
167 # + f"{salmon_res_dir}"
|
|
kkonganti@11
|
168 # + "\ndoes not seem to have any directories that contain\n"
|
|
kkonganti@11
|
169 # + "the `quant.sf` file(s)."
|
|
kkonganti@11
|
170 # )
|
|
kkonganti@11
|
171 # exit(1)
|
|
kkonganti@11
|
172 with open(salmon_comb_res, "w") as salmon_comb_res_fh:
|
|
kkonganti@11
|
173 salmon_comb_res_fh.write(f"Sample\n{no_hit}s in any samples\n")
|
|
kkonganti@11
|
174 salmon_comb_res_fh.close()
|
|
kkonganti@11
|
175
|
|
kkonganti@11
|
176 with open(bcs_sal_yn, "w") as bcs_sal_yn_fh:
|
|
kkonganti@11
|
177 bcs_sal_yn_fh.write(f"Sample\n{no_presence} in any samples\n")
|
|
kkonganti@11
|
178 bcs_sal_yn_fh.close()
|
|
kkonganti@11
|
179
|
|
kkonganti@11
|
180 exit(0)
|
|
kkonganti@11
|
181
|
|
kkonganti@11
|
182 if rtc and os.path.exists(rtc) and os.path.getsize(rtc) > 0:
|
|
kkonganti@11
|
183
|
|
kkonganti@11
|
184 # pdg_release = re.match(r"(^PDG\d+\.\d+)\..+", os.path.basename(rtc))[1] + "/"
|
|
kkonganti@11
|
185 acc2sero = pickle.load(file=open(pickled_sero, "rb"))
|
|
kkonganti@11
|
186
|
|
kkonganti@11
|
187 with open(rtc, "r") as rtc_fh:
|
|
kkonganti@11
|
188
|
|
kkonganti@11
|
189 for line in rtc_fh:
|
|
kkonganti@11
|
190 cols = line.strip().split("\t")
|
|
kkonganti@11
|
191
|
|
kkonganti@11
|
192 if len(cols) < 4:
|
|
kkonganti@11
|
193 logging.error(
|
|
kkonganti@11
|
194 f"The file {os.path.basename(rtc)} seems to\n"
|
|
kkonganti@11
|
195 + "be malformed. It contains less than required 4 columns."
|
|
kkonganti@11
|
196 )
|
|
kkonganti@11
|
197 exit(1)
|
|
kkonganti@11
|
198 elif cols[3] != null_value:
|
|
kkonganti@11
|
199 snp_clusters[cols[0]].setdefault("assembly_accs", []).append(cols[3])
|
|
kkonganti@11
|
200 snp_clusters[cols[3]].setdefault("snp_clust_id", []).append(cols[0])
|
|
kkonganti@11
|
201 snp_clusters[cols[3]].setdefault("pathdb_acc_id", []).append(cols[1])
|
|
kkonganti@11
|
202 if len(snp_clusters[cols[3]]["snp_clust_id"]) > 1:
|
|
kkonganti@11
|
203 logging.error(
|
|
kkonganti@11
|
204 f"There is a duplicate reference accession [{cols[3]}]"
|
|
kkonganti@11
|
205 + f"in the metadata file{os.path.basename(rtc)}!"
|
|
kkonganti@11
|
206 )
|
|
kkonganti@11
|
207 exit(1)
|
|
kkonganti@11
|
208
|
|
kkonganti@11
|
209 rtc_fh.close()
|
|
kkonganti@11
|
210
|
|
kkonganti@11
|
211 for salmon_res_file in salmon_res_files:
|
|
kkonganti@11
|
212 sample_name = re.match(
|
|
kkonganti@11
|
213 r"(^.+?)((\_salmon\_res)|(\.salmon))$",
|
|
kkonganti@11
|
214 os.path.basename(os.path.dirname(salmon_res_file)),
|
|
kkonganti@11
|
215 )[1]
|
|
kkonganti@11
|
216 salmon_meta_json = os.path.join(
|
|
kkonganti@11
|
217 os.path.dirname(salmon_res_file), "aux_info", "meta_info.json"
|
|
kkonganti@11
|
218 )
|
|
kkonganti@11
|
219
|
|
kkonganti@11
|
220 if not os.path.exists(salmon_meta_json) or not os.path.getsize(salmon_meta_json) > 0:
|
|
kkonganti@11
|
221 logging.error(
|
|
kkonganti@11
|
222 "The file\n"
|
|
kkonganti@11
|
223 + f"{salmon_meta_json}\ndoes not exist or is empty!\n"
|
|
kkonganti@11
|
224 + "Did `salmon quant` fail?"
|
|
kkonganti@11
|
225 )
|
|
kkonganti@11
|
226 exit(1)
|
|
kkonganti@11
|
227
|
|
kkonganti@11
|
228 if not os.path.exists(salmon_res_file) or not os.path.getsize(salmon_res_file):
|
|
kkonganti@11
|
229 logging.error(
|
|
kkonganti@11
|
230 "The file\n"
|
|
kkonganti@11
|
231 + f"{salmon_res_file}\ndoes not exist or is empty!\n"
|
|
kkonganti@11
|
232 + "Did `salmon quant` fail?"
|
|
kkonganti@11
|
233 )
|
|
kkonganti@11
|
234 exit(1)
|
|
kkonganti@11
|
235
|
|
kkonganti@11
|
236 with open(salmon_res_file, "r") as salmon_res_fh:
|
|
kkonganti@11
|
237 for line in salmon_res_fh.readlines():
|
|
kkonganti@11
|
238 if re.match(r"^Name.+", line):
|
|
kkonganti@11
|
239 continue
|
|
kkonganti@11
|
240 cols = line.strip().split("\t")
|
|
kkonganti@11
|
241 ref_acc = "_".join(cols[0].split("_")[:2])
|
|
kkonganti@11
|
242
|
|
kkonganti@11
|
243 if ref_acc not in snp_clusters.keys():
|
|
kkonganti@11
|
244 snp_clusters[ref_acc]["snp_clust_id"] = ref_acc
|
|
kkonganti@11
|
245 snp_clusters[ref_acc]["pathdb_acc_id"] = ref_acc
|
|
kkonganti@11
|
246
|
|
kkonganti@11
|
247 (
|
|
kkonganti@11
|
248 sample2salmon[sample_name]
|
|
kkonganti@11
|
249 .setdefault(acc2sero[cols[0]], [])
|
|
kkonganti@11
|
250 .append(int(round(float(cols[4]), 2)))
|
|
kkonganti@11
|
251 )
|
|
kkonganti@11
|
252 (
|
|
kkonganti@11
|
253 sample2salmon[sample_name]
|
|
kkonganti@11
|
254 .setdefault("snp_clust_ids", {})
|
|
kkonganti@11
|
255 .setdefault("".join(snp_clusters[ref_acc]["snp_clust_id"]), [])
|
|
kkonganti@11
|
256 .append("".join(snp_clusters[ref_acc]["pathdb_acc_id"]))
|
|
kkonganti@11
|
257 )
|
|
kkonganti@11
|
258 seen_sero[acc2sero[cols[0]]] = 1
|
|
kkonganti@11
|
259
|
|
kkonganti@11
|
260 salmon_meta_json_read = json.load(open(salmon_meta_json, "r"))
|
|
kkonganti@11
|
261 (
|
|
kkonganti@11
|
262 sample2salmon[sample_name]
|
|
kkonganti@11
|
263 .setdefault("tot_reads", [])
|
|
kkonganti@11
|
264 .append(salmon_meta_json_read["num_processed"])
|
|
kkonganti@11
|
265 )
|
|
kkonganti@11
|
266
|
|
kkonganti@11
|
267 with open(salmon_comb_res, "w") as salmon_comb_res_fh:
|
|
kkonganti@11
|
268
|
|
kkonganti@11
|
269 # snp_clust_col_header = (
|
|
kkonganti@11
|
270 # "\tSNP Cluster(s) by Genome Hit\n" if show_snp_clust_col else "\n"
|
|
kkonganti@11
|
271 # )
|
|
kkonganti@11
|
272 snp_clust_col_header = (
|
|
kkonganti@11
|
273 "\tNCBI Pathogens Isolate Browser\n" if show_snp_clust_col else "\n"
|
|
kkonganti@11
|
274 )
|
|
kkonganti@11
|
275 serotypes = sorted(seen_sero.keys())
|
|
kkonganti@11
|
276 formatted_serotypes = [
|
|
kkonganti@11
|
277 re.sub(r"\,antigen_formula=", " | ", s)
|
|
kkonganti@11
|
278 for s in [re.sub(r"serotype=", "", s) for s in serotypes]
|
|
kkonganti@11
|
279 ]
|
|
kkonganti@11
|
280 salmon_comb_res_fh.write(
|
|
kkonganti@11
|
281 "Sample\t" + "\t".join(formatted_serotypes) + snp_clust_col_header
|
|
kkonganti@11
|
282 )
|
|
kkonganti@11
|
283 # sample_snp_relation = (
|
|
kkonganti@11
|
284 # ncbi_pathogens_base_url
|
|
kkonganti@11
|
285 # + pdg_release
|
|
kkonganti@11
|
286 # + "".join(snp_clusters[ref_acc]["snp_clust_id"])
|
|
kkonganti@11
|
287 # + "?accessions="
|
|
kkonganti@11
|
288 # )
|
|
kkonganti@11
|
289 if len(salmon_res_file_failed) == 1:
|
|
kkonganti@11
|
290 with (open("".join(salmon_res_file_failed), "r")) as no_calls_fh:
|
|
kkonganti@11
|
291 for line in no_calls_fh.readlines():
|
|
kkonganti@11
|
292 if line in ["\n", "\n\r", "\r"]:
|
|
kkonganti@11
|
293 continue
|
|
kkonganti@11
|
294 salmon_comb_res_fh.write(line.strip())
|
|
kkonganti@11
|
295 sal_yn[line.strip()] += 0
|
|
kkonganti@11
|
296 for serotype in serotypes:
|
|
kkonganti@11
|
297 salmon_comb_res_fh.write("\t-")
|
|
kkonganti@11
|
298 salmon_comb_res_fh.write(
|
|
kkonganti@11
|
299 "\t-\n"
|
|
kkonganti@11
|
300 ) if show_snp_clust_col else salmon_comb_res_fh.write("\n")
|
|
kkonganti@11
|
301 no_calls_fh.close()
|
|
kkonganti@11
|
302
|
|
kkonganti@11
|
303 for sample, counts in sorted(sample2salmon.items()):
|
|
kkonganti@11
|
304 salmon_comb_res_fh.write(sample)
|
|
kkonganti@11
|
305 snp_cluster_res_col = list()
|
|
kkonganti@11
|
306
|
|
kkonganti@11
|
307 for snp_clust_id in sample2salmon[sample]["snp_clust_ids"].keys():
|
|
kkonganti@11
|
308 # print(snp_clust_id)
|
|
kkonganti@11
|
309 # print(",".join(sample2salmon[sample]["snp_clust_ids"][snp_clust_id]))
|
|
kkonganti@11
|
310 # ppp.pprint(sample2salmon[sample]["snp_clust_ids"])
|
|
kkonganti@11
|
311 # ppp.pprint(sample2salmon[sample]["snp_clust_ids"][snp_clust_id])
|
|
kkonganti@11
|
312 # final_url_text = ",".join(
|
|
kkonganti@11
|
313 # sample2salmon[sample]["snp_clust_ids"][snp_clust_id]
|
|
kkonganti@11
|
314 # )
|
|
kkonganti@11
|
315 # final_url_text_to_show = snp_clust_id
|
|
kkonganti@11
|
316 # snp_cluster_res_col.append(
|
|
kkonganti@11
|
317 # "".join(
|
|
kkonganti@11
|
318 # [
|
|
kkonganti@11
|
319 # f'<a href="',
|
|
kkonganti@11
|
320 # sample_snp_relation,
|
|
kkonganti@11
|
321 # ",".join(sample2salmon[sample]["snp_clust_ids"][snp_clust_id]),
|
|
kkonganti@11
|
322 # f'" target="_blank">{snp_clust_id}</a>',
|
|
kkonganti@11
|
323 # ]
|
|
kkonganti@11
|
324 # )
|
|
kkonganti@11
|
325 # )
|
|
kkonganti@11
|
326 # ppp.pprint(sample2salmon[sample])
|
|
kkonganti@11
|
327 for pathdbacc in sample2salmon[sample]["snp_clust_ids"][snp_clust_id]:
|
|
kkonganti@11
|
328 # final_url_text_to_show = " ".join(
|
|
kkonganti@11
|
329 # sample2salmon[sample]["snp_clust_ids"][snp_clust_id]
|
|
kkonganti@11
|
330 # )
|
|
kkonganti@11
|
331 sample_snp_relation = (
|
|
kkonganti@11
|
332 ncbi_pathogens_genome_base
|
|
kkonganti@11
|
333 if assm_pat.match(pathdbacc)
|
|
kkonganti@11
|
334 else ncbi_pathogens_base_url + "isolates/#"
|
|
kkonganti@11
|
335 )
|
|
kkonganti@11
|
336
|
|
kkonganti@11
|
337 snp_cluster_res_col.append(
|
|
kkonganti@11
|
338 "".join(
|
|
kkonganti@11
|
339 [
|
|
kkonganti@11
|
340 f'<a href="',
|
|
kkonganti@11
|
341 sample_snp_relation,
|
|
kkonganti@11
|
342 pathdbacc,
|
|
kkonganti@11
|
343 f'" target="_blank">{pathdbacc}</a>',
|
|
kkonganti@11
|
344 ]
|
|
kkonganti@11
|
345 )
|
|
kkonganti@11
|
346 )
|
|
kkonganti@11
|
347
|
|
kkonganti@11
|
348 per_serotype_counts = 0
|
|
kkonganti@11
|
349 for serotype in serotypes:
|
|
kkonganti@11
|
350
|
|
kkonganti@11
|
351 if serotype in sample2salmon[sample].keys():
|
|
kkonganti@11
|
352 # ppp.pprint(counts)
|
|
kkonganti@11
|
353 sample_perc_mapped = round(
|
|
kkonganti@11
|
354 sum(counts[serotype]) / sum(counts["tot_reads"]) * 100, 2
|
|
kkonganti@11
|
355 )
|
|
kkonganti@11
|
356 salmon_comb_res_fh.write(
|
|
kkonganti@11
|
357 f"\t{sum(counts[serotype])} ({sample_perc_mapped}%)"
|
|
kkonganti@11
|
358 )
|
|
kkonganti@11
|
359 multiqc_salmon_counts[sample].setdefault(
|
|
kkonganti@11
|
360 re.match(r"^serotype=(.+?)\,antigen_formula.*", serotype)[1],
|
|
kkonganti@11
|
361 sum(counts[serotype]),
|
|
kkonganti@11
|
362 )
|
|
kkonganti@11
|
363 per_serotype_counts += sum(counts[serotype])
|
|
kkonganti@11
|
364 sal_yn[sample] += 1
|
|
kkonganti@11
|
365 else:
|
|
kkonganti@11
|
366 salmon_comb_res_fh.write(f"\t-")
|
|
kkonganti@11
|
367 sal_yn[sample] += 0
|
|
kkonganti@11
|
368
|
|
kkonganti@11
|
369 multiqc_salmon_counts[sample].setdefault(
|
|
kkonganti@11
|
370 no_hit, sum(counts["tot_reads"]) - per_serotype_counts
|
|
kkonganti@11
|
371 )
|
|
kkonganti@11
|
372 snp_clust_col_val = (
|
|
kkonganti@11
|
373 f'\t{" ".join(snp_cluster_res_col)}\n' if show_snp_clust_col else "\n"
|
|
kkonganti@11
|
374 )
|
|
kkonganti@11
|
375 # ppp.pprint(multiqc_salmon_counts)
|
|
kkonganti@11
|
376 salmon_comb_res_fh.write(snp_clust_col_val)
|
|
kkonganti@11
|
377
|
|
kkonganti@11
|
378 with open(bcs_sal_yn, "w") as bcs_sal_yn_fh:
|
|
kkonganti@11
|
379 bcs_sal_yn_fh.write("Sample\tSalmonella Presence\tNo. of Serotypes\n")
|
|
kkonganti@11
|
380 for sample in sal_yn.keys():
|
|
kkonganti@11
|
381 if sal_yn[sample] > 0:
|
|
kkonganti@11
|
382 bcs_sal_yn_fh.write(f"{sample}\tDetected\t{sal_yn[sample]}\n")
|
|
kkonganti@11
|
383 else:
|
|
kkonganti@11
|
384 bcs_sal_yn_fh.write(f"{sample}\tNot detected\t{sal_yn[sample]}\n")
|
|
kkonganti@11
|
385
|
|
kkonganti@11
|
386 with open(cell_colors_yml_file, "w") as cell_colors_fh:
|
|
kkonganti@11
|
387 yaml.dump(cell_colors_yml, cell_colors_fh, default_flow_style=False)
|
|
kkonganti@11
|
388
|
|
kkonganti@11
|
389 salmon_plot_json(salmon_comb_res_mqc, multiqc_salmon_counts, no_hit)
|
|
kkonganti@11
|
390
|
|
kkonganti@11
|
391 salmon_comb_res_fh.close()
|
|
kkonganti@11
|
392 bcs_sal_yn_fh.close()
|
|
kkonganti@11
|
393 cell_colors_fh.close()
|
|
kkonganti@11
|
394
|
|
kkonganti@11
|
395
|
|
kkonganti@11
|
396 def salmon_plot_json(file: None, sample_salmon_counts: None, no_hit: None) -> None:
|
|
kkonganti@11
|
397 """
|
|
kkonganti@11
|
398 This method will take a dictionary of salmon counts per sample
|
|
kkonganti@11
|
399 and will dump a JSON that will be used by MultiQC.
|
|
kkonganti@11
|
400 """
|
|
kkonganti@11
|
401
|
|
kkonganti@11
|
402 if file is None or sample_salmon_counts is None:
|
|
kkonganti@11
|
403 logging.error(
|
|
kkonganti@11
|
404 "Neither an output file to dump the JSON for MultiQC or the"
|
|
kkonganti@11
|
405 + "dictionary holding the salmon counts was not passed."
|
|
kkonganti@11
|
406 )
|
|
kkonganti@11
|
407
|
|
kkonganti@11
|
408 # Credit: http://phrogz.net/tmp/24colors.html
|
|
kkonganti@11
|
409 # Will cycle through 20 distinct colors.
|
|
kkonganti@11
|
410 distinct_color_palette = [
|
|
kkonganti@11
|
411 "#FF0000",
|
|
kkonganti@11
|
412 "#FFFF00",
|
|
kkonganti@11
|
413 "#00EAFF",
|
|
kkonganti@11
|
414 "#AA00FF",
|
|
kkonganti@11
|
415 "#FF7F00",
|
|
kkonganti@11
|
416 "#BFFF00",
|
|
kkonganti@11
|
417 "#0095FF",
|
|
kkonganti@11
|
418 "#FF00AA",
|
|
kkonganti@11
|
419 "#FFD400",
|
|
kkonganti@11
|
420 "#6AFF00",
|
|
kkonganti@11
|
421 "#0040FF",
|
|
kkonganti@11
|
422 "#EDB9B9",
|
|
kkonganti@11
|
423 "#B9D7ED",
|
|
kkonganti@11
|
424 "#E7E9B9",
|
|
kkonganti@11
|
425 "#DCB9ED",
|
|
kkonganti@11
|
426 "#B9EDE0",
|
|
kkonganti@11
|
427 "#8F2323",
|
|
kkonganti@11
|
428 "#23628F",
|
|
kkonganti@11
|
429 "#8F6A23",
|
|
kkonganti@11
|
430 "#6B238F",
|
|
kkonganti@11
|
431 "#4F8F23",
|
|
kkonganti@11
|
432 ]
|
|
kkonganti@11
|
433
|
|
kkonganti@11
|
434 # Credit: https://mokole.com/palette.html
|
|
kkonganti@11
|
435 # Will use this palette if we run out ouf
|
|
kkonganti@11
|
436 # 20 serotypes. More than 50 serotypes
|
|
kkonganti@11
|
437 # per run is probably rare but if not,
|
|
kkonganti@11
|
438 # will cycle through about 45.
|
|
kkonganti@11
|
439 distinct_color_palette2 = [
|
|
kkonganti@11
|
440 "#2F4F4F", # darkslategray
|
|
kkonganti@11
|
441 "#556B2F", # darkolivegreen
|
|
kkonganti@11
|
442 "#A0522D", # sienna
|
|
kkonganti@11
|
443 "#2E8B57", # seagreen
|
|
kkonganti@11
|
444 "#006400", # darkgreen
|
|
kkonganti@11
|
445 "#8B0000", # darkred
|
|
kkonganti@11
|
446 "#808000", # olive
|
|
kkonganti@11
|
447 "#BC8F8F", # rosybrown
|
|
kkonganti@11
|
448 "#663399", # rebeccapurple
|
|
kkonganti@11
|
449 "#B8860B", # darkgoldenrod
|
|
kkonganti@11
|
450 "#4682B4", # steelblue
|
|
kkonganti@11
|
451 "#000080", # navy
|
|
kkonganti@11
|
452 "#D2691E", # chocolate
|
|
kkonganti@11
|
453 "#9ACD32", # yellowgreen
|
|
kkonganti@11
|
454 "#20B2AA", # lightseagreen
|
|
kkonganti@11
|
455 "#CD5C5C", # indianred
|
|
kkonganti@11
|
456 "#8FBC8F", # darkseagreen
|
|
kkonganti@11
|
457 "#800080", # purple
|
|
kkonganti@11
|
458 "#B03060", # maroon3
|
|
kkonganti@11
|
459 "#FF8C00", # darkorange
|
|
kkonganti@11
|
460 "#FFD700", # gold
|
|
kkonganti@11
|
461 "#FFFF00", # yellow
|
|
kkonganti@11
|
462 "#DEB887", # burlywood
|
|
kkonganti@11
|
463 "#00FF00", # lime
|
|
kkonganti@11
|
464 "#BA55D3", # mediumorchid
|
|
kkonganti@11
|
465 "#00FA9A", # mediumspringgreen
|
|
kkonganti@11
|
466 "#4169E1", # royalblue
|
|
kkonganti@11
|
467 "#E9967A", # darksalmon
|
|
kkonganti@11
|
468 "#DC143C", # crimson
|
|
kkonganti@11
|
469 "#00FFFF", # aqua
|
|
kkonganti@11
|
470 "#F4A460", # sandybrown
|
|
kkonganti@11
|
471 "#9370DB", # mediumpurple
|
|
kkonganti@11
|
472 "#0000FF", # blue
|
|
kkonganti@11
|
473 "#ADFF2F", # greenyellow
|
|
kkonganti@11
|
474 "#FF6347", # tomato
|
|
kkonganti@11
|
475 "#D8BFD8", # thistle
|
|
kkonganti@11
|
476 "#FF00FF", # fuchsia
|
|
kkonganti@11
|
477 "#DB7093", # palevioletred
|
|
kkonganti@11
|
478 "#F0E68C", # khaki
|
|
kkonganti@11
|
479 "#6495ED", # cornflower
|
|
kkonganti@11
|
480 "#DDA0DD", # plum
|
|
kkonganti@11
|
481 "#EE82EE", # violet
|
|
kkonganti@11
|
482 "#7FFFD4", # aquamarine
|
|
kkonganti@11
|
483 "#FAFAD2", # lightgoldenrod
|
|
kkonganti@11
|
484 "#FF69B4", # hotpink
|
|
kkonganti@11
|
485 "#FFB6C1", # lightpink
|
|
kkonganti@11
|
486 ]
|
|
kkonganti@11
|
487
|
|
kkonganti@11
|
488 no_hit_color = "#434348"
|
|
kkonganti@11
|
489 col_count = 0
|
|
kkonganti@11
|
490 serotypes = set()
|
|
kkonganti@11
|
491 salmon_counts = defaultdict(defaultdict)
|
|
kkonganti@11
|
492 salmon_counts["id"] = "BETTERCALLSAL_SALMON_COUNTS"
|
|
kkonganti@11
|
493 salmon_counts["section_name"] = "Salmon read counts"
|
|
kkonganti@11
|
494 salmon_counts["description"] = (
|
|
kkonganti@11
|
495 "This section shows the read counts from running <code>salmon</code> "
|
|
kkonganti@11
|
496 + "in <code>--meta</code> mode using SE, merged PE or concatenated PE reads against "
|
|
kkonganti@11
|
497 + "an on-the-fly <code>salmon</code> index generated from the genome hits "
|
|
kkonganti@11
|
498 + "of <code>kma</code>."
|
|
kkonganti@11
|
499 )
|
|
kkonganti@11
|
500 salmon_counts["plot_type"] = "bargraph"
|
|
kkonganti@11
|
501 salmon_counts["pconfig"]["id"] = "bettercallsal_salmon_counts_plot"
|
|
kkonganti@11
|
502 salmon_counts["pconfig"]["title"] = "Salmon: Read counts"
|
|
kkonganti@11
|
503 salmon_counts["pconfig"]["ylab"] = "Number of reads"
|
|
kkonganti@11
|
504 salmon_counts["pconfig"]["xDecimals"] = "false"
|
|
kkonganti@11
|
505 salmon_counts["pconfig"]["cpswitch_counts_label"] = "Number of reads (Counts)"
|
|
kkonganti@11
|
506 salmon_counts["pconfig"]["cpswitch_percent_label"] = "Number of reads (Percentages)"
|
|
kkonganti@11
|
507
|
|
kkonganti@11
|
508 for sample in sorted(sample_salmon_counts.keys()):
|
|
kkonganti@11
|
509 serotypes.update(list(sample_salmon_counts[sample].keys()))
|
|
kkonganti@11
|
510 salmon_counts["data"][sample] = sample_salmon_counts[sample]
|
|
kkonganti@11
|
511
|
|
kkonganti@11
|
512 if len(serotypes) > len(distinct_color_palette):
|
|
kkonganti@11
|
513 distinct_color_palette = distinct_color_palette2
|
|
kkonganti@11
|
514
|
|
kkonganti@11
|
515 for serotype in sorted(serotypes):
|
|
kkonganti@11
|
516 if serotype == no_hit:
|
|
kkonganti@11
|
517 continue
|
|
kkonganti@11
|
518 if col_count == len(distinct_color_palette) - 1:
|
|
kkonganti@11
|
519 col_count = 0
|
|
kkonganti@11
|
520
|
|
kkonganti@11
|
521 col_count += 1
|
|
kkonganti@11
|
522 salmon_counts["categories"][serotype] = {"color": distinct_color_palette[col_count]}
|
|
kkonganti@11
|
523
|
|
kkonganti@11
|
524 salmon_counts["categories"][no_hit] = {"color": no_hit_color}
|
|
kkonganti@11
|
525 json.dump(salmon_counts, open(file, "w"))
|
|
kkonganti@11
|
526
|
|
kkonganti@11
|
527
|
|
kkonganti@11
|
528 if __name__ == "__main__":
|
|
kkonganti@11
|
529 main()
|
