kkonganti@11: #!/usr/bin/env python3
kkonganti@11:
kkonganti@11: # Kranti Konganti
kkonganti@11:
kkonganti@11: import argparse
kkonganti@11: import glob
kkonganti@11: import inspect
kkonganti@11: import json
kkonganti@11: import logging
kkonganti@11: import os
kkonganti@11: import pickle
kkonganti@11: import pprint
kkonganti@11: import re
kkonganti@11: from collections import defaultdict
kkonganti@11:
kkonganti@11: import yaml
kkonganti@11:
kkonganti@11:
kkonganti@11: # Multiple inheritence for pretty printing of help text.
kkonganti@11: class MultiArgFormatClasses(argparse.RawTextHelpFormatter, argparse.ArgumentDefaultsHelpFormatter):
kkonganti@11: pass
kkonganti@11:
kkonganti@11:
kkonganti@11: # Main
kkonganti@11: def main() -> None:
kkonganti@11: """
kkonganti@11: The succesful execution of this script requires access to bettercallsal formatted
kkonganti@11: db flat files. On raven2, they are at /hpc/db/bettercallsall/PDGXXXXXXXXXX.XXXXX
kkonganti@11:
kkonganti@11: It takes the ACC2SERO.pickle file and *.reference_target.cluster_list.tsv file
kkonganti@11: for that particular NCBI Pathogens release from the db directory mentioned with
kkonganti@11: -db option and a root parent directory of the `salmon quant` results mentioned
kkonganti@11: with -sal option and generates a final results table with number of reads
kkonganti@11: mapped and a .json file to be used with MultiQC to generate a stacked bar plot.
kkonganti@11:
kkonganti@11: Using -url option optionally adds an extra column of NCBI Pathogens Isolates
kkonganti@11: Browser, which directly links out to NCBI Pathogens Isolates SNP viewer tool.
kkonganti@11: """
kkonganti@11: # Set logging.
kkonganti@11: logging.basicConfig(
kkonganti@11: format="\n" + "=" * 55 + "\n%(asctime)s - %(levelname)s\n" + "=" * 55 + "\n%(message)s\n\n",
kkonganti@11: level=logging.DEBUG,
kkonganti@11: )
kkonganti@11:
kkonganti@11: # Debug print.
kkonganti@11: ppp = pprint.PrettyPrinter(width=55)
kkonganti@11: prog_name = inspect.stack()[0].filename
kkonganti@11:
kkonganti@11: parser = argparse.ArgumentParser(
kkonganti@11: prog=prog_name, description=main.__doc__, formatter_class=MultiArgFormatClasses
kkonganti@11: )
kkonganti@11:
kkonganti@11: required = parser.add_argument_group("required arguments")
kkonganti@11:
kkonganti@11: required.add_argument(
kkonganti@11: "-sal",
kkonganti@11: dest="salmon_res_dir",
kkonganti@11: default=False,
kkonganti@11: required=True,
kkonganti@11: help="Absolute UNIX path to the parent directory that contains the\n"
kkonganti@11: + "`salmon quant` results directory. For example, if path to\n"
kkonganti@11: + "`quant.sf` is in /hpc/john_doe/test/salmon_res/quant.sf, then\n"
kkonganti@11: + "use this command-line option as:\n"
kkonganti@11: + "-sal /hpc/john_doe/test",
kkonganti@11: )
kkonganti@11: required.add_argument(
kkonganti@11: "-snp",
kkonganti@11: dest="rtc",
kkonganti@11: default=False,
kkonganti@11: required=True,
kkonganti@11: help="Absolute UNIX Path to the PDG SNP reference target cluster\n"
kkonganti@11: + "metadata file. On raven2, these are located at\n"
kkonganti@11: + "/hpc/db/bettercallsal/PDGXXXXXXXXXX.XXXXX\n"
kkonganti@11: + "Required if -sal is on.",
kkonganti@11: )
kkonganti@11: required.add_argument(
kkonganti@11: "-pickle",
kkonganti@11: dest="acc2sero",
kkonganti@11: default=False,
kkonganti@11: required=True,
kkonganti@11: help="Absolute UNIX Path to the *ACC2SERO.pickle\n"
kkonganti@11: + "metadata file. On raven2, these are located at\n"
kkonganti@11: + "/hpc/db/bettercallsal/PDGXXXXXXXXXX.XXXXX\n"
kkonganti@11: + "Required if -sal is on.",
kkonganti@11: )
kkonganti@11: parser.add_argument(
kkonganti@11: "-op",
kkonganti@11: dest="out_prefix",
kkonganti@11: default="bettercallsal.tblsum",
kkonganti@11: required=False,
kkonganti@11: help="Set the output file(s) prefix for output(s) generated\n" + "by this program.",
kkonganti@11: )
kkonganti@11: parser.add_argument(
kkonganti@11: "-url",
kkonganti@11: dest="show_snp_clust_info",
kkonganti@11: default=False,
kkonganti@11: required=False,
kkonganti@11: action="store_true",
kkonganti@11: help="Show SNP cluster participation information of the final genome hit.\n"
kkonganti@11: + "This may be useful to see a relative placement of your sample in\n"
kkonganti@11: + "NCBI Isolates SNP Tree Viewer based on genome similarity but however\n"
kkonganti@11: + "due to rapid nature of the updates at NCBI Pathogen Detection Project,\n"
kkonganti@11: + "the placement may be in an outdated cluster.",
kkonganti@11: )
kkonganti@11:
kkonganti@11: args = parser.parse_args()
kkonganti@11: salmon_res_dir = args.salmon_res_dir
kkonganti@11: out_prefix = args.out_prefix
kkonganti@11: show_snp_clust_col = args.show_snp_clust_info
kkonganti@11: rtc = args.rtc
kkonganti@11: pickled_sero = args.acc2sero
kkonganti@11: no_hit = "No genome hit"
kkonganti@11: no_presence = "Salmonella presence not detected"
kkonganti@11: bcs_sal_yn_prefix = "bettercallsal_salyn"
kkonganti@11: sal_y = "Detected"
kkonganti@11: sal_n = "Not detected"
kkonganti@11: null_value = "NULL"
kkonganti@11: assm_pat = re.compile(r"GC[AF]\_[0-9]+\.*[0-9]*")
kkonganti@11: ncbi_pathogens_base_url = "https://www.ncbi.nlm.nih.gov/pathogens/"
kkonganti@11: ncbi_pathogens_genome_base = "https://www.ncbi.nlm.nih.gov/datasets/genome/"
kkonganti@11:
kkonganti@11: sample2salmon, snp_clusters, multiqc_salmon_counts, seen_sero, sal_yn = (
kkonganti@11: defaultdict(defaultdict),
kkonganti@11: defaultdict(defaultdict),
kkonganti@11: defaultdict(defaultdict),
kkonganti@11: defaultdict(int),
kkonganti@11: defaultdict(int),
kkonganti@11: )
kkonganti@11:
kkonganti@11: cell_colors_yml = {
kkonganti@11: bcs_sal_yn_prefix: {sal_y: "#c8e6c9 !important;", sal_n: "#ffcdd2 !important;"}
kkonganti@11: }
kkonganti@11:
kkonganti@11: salmon_comb_res = os.path.join(os.getcwd(), out_prefix + ".txt")
kkonganti@11: bcs_sal_yn = re.sub(out_prefix, bcs_sal_yn_prefix + ".tblsum", salmon_comb_res)
kkonganti@11: cell_colors_yml_file = re.sub(
kkonganti@11: out_prefix + ".txt", bcs_sal_yn_prefix + ".cellcolors.yml", salmon_comb_res
kkonganti@11: )
kkonganti@11: salmon_comb_res_mqc = os.path.join(os.getcwd(), str(out_prefix).split(".")[0] + "_mqc.json")
kkonganti@11: salmon_res_files = glob.glob(os.path.join(salmon_res_dir, "*", "quant.sf"), recursive=True)
kkonganti@11: salmon_res_file_failed = glob.glob(os.path.join(salmon_res_dir, "BCS_NO_CALLS.txt"))
kkonganti@11:
kkonganti@11: if rtc and (not os.path.exists(rtc) or not os.path.getsize(rtc) > 0):
kkonganti@11: logging.error(
kkonganti@11: "The reference target cluster metadata file,\n"
kkonganti@11: + f"{os.path.basename(rtc)} does not exist or is empty!"
kkonganti@11: )
kkonganti@11: exit(1)
kkonganti@11:
kkonganti@11: if rtc and (not salmon_res_dir or not pickled_sero):
kkonganti@11: logging.error("When -rtc is on, -sal and -ps are also required.")
kkonganti@11: exit(1)
kkonganti@11:
kkonganti@11: if pickled_sero and (not os.path.exists(pickled_sero) or not os.path.getsize(pickled_sero)):
kkonganti@11: logging.error(
kkonganti@11: "The pickle file,\n" + f"{os.path.basename(pickled_sero)} does not exist or is empty!"
kkonganti@11: )
kkonganti@11: exit(1)
kkonganti@11:
kkonganti@11: if salmon_res_dir:
kkonganti@11: if not os.path.isdir(salmon_res_dir):
kkonganti@11: logging.error("UNIX path\n" + f"{salmon_res_dir}\n" + "does not exist!")
kkonganti@11: exit(1)
kkonganti@11: if len(salmon_res_files) <= 0:
kkonganti@11: # logging.error(
kkonganti@11: # "Parent directory,\n"
kkonganti@11: # + f"{salmon_res_dir}"
kkonganti@11: # + "\ndoes not seem to have any directories that contain\n"
kkonganti@11: # + "the `quant.sf` file(s)."
kkonganti@11: # )
kkonganti@11: # exit(1)
kkonganti@11: with open(salmon_comb_res, "w") as salmon_comb_res_fh:
kkonganti@11: salmon_comb_res_fh.write(f"Sample\n{no_hit}s in any samples\n")
kkonganti@11: salmon_comb_res_fh.close()
kkonganti@11:
kkonganti@11: with open(bcs_sal_yn, "w") as bcs_sal_yn_fh:
kkonganti@11: bcs_sal_yn_fh.write(f"Sample\n{no_presence} in any samples\n")
kkonganti@11: bcs_sal_yn_fh.close()
kkonganti@11:
kkonganti@11: exit(0)
kkonganti@11:
kkonganti@11: if rtc and os.path.exists(rtc) and os.path.getsize(rtc) > 0:
kkonganti@11:
kkonganti@11: # pdg_release = re.match(r"(^PDG\d+\.\d+)\..+", os.path.basename(rtc))[1] + "/"
kkonganti@11: acc2sero = pickle.load(file=open(pickled_sero, "rb"))
kkonganti@11:
kkonganti@11: with open(rtc, "r") as rtc_fh:
kkonganti@11:
kkonganti@11: for line in rtc_fh:
kkonganti@11: cols = line.strip().split("\t")
kkonganti@11:
kkonganti@11: if len(cols) < 4:
kkonganti@11: logging.error(
kkonganti@11: f"The file {os.path.basename(rtc)} seems to\n"
kkonganti@11: + "be malformed. It contains less than required 4 columns."
kkonganti@11: )
kkonganti@11: exit(1)
kkonganti@11: elif cols[3] != null_value:
kkonganti@11: snp_clusters[cols[0]].setdefault("assembly_accs", []).append(cols[3])
kkonganti@11: snp_clusters[cols[3]].setdefault("snp_clust_id", []).append(cols[0])
kkonganti@11: snp_clusters[cols[3]].setdefault("pathdb_acc_id", []).append(cols[1])
kkonganti@11: if len(snp_clusters[cols[3]]["snp_clust_id"]) > 1:
kkonganti@11: logging.error(
kkonganti@11: f"There is a duplicate reference accession [{cols[3]}]"
kkonganti@11: + f"in the metadata file{os.path.basename(rtc)}!"
kkonganti@11: )
kkonganti@11: exit(1)
kkonganti@11:
kkonganti@11: rtc_fh.close()
kkonganti@11:
kkonganti@11: for salmon_res_file in salmon_res_files:
kkonganti@11: sample_name = re.match(
kkonganti@11: r"(^.+?)((\_salmon\_res)|(\.salmon))$",
kkonganti@11: os.path.basename(os.path.dirname(salmon_res_file)),
kkonganti@11: )[1]
kkonganti@11: salmon_meta_json = os.path.join(
kkonganti@11: os.path.dirname(salmon_res_file), "aux_info", "meta_info.json"
kkonganti@11: )
kkonganti@11:
kkonganti@11: if not os.path.exists(salmon_meta_json) or not os.path.getsize(salmon_meta_json) > 0:
kkonganti@11: logging.error(
kkonganti@11: "The file\n"
kkonganti@11: + f"{salmon_meta_json}\ndoes not exist or is empty!\n"
kkonganti@11: + "Did `salmon quant` fail?"
kkonganti@11: )
kkonganti@11: exit(1)
kkonganti@11:
kkonganti@11: if not os.path.exists(salmon_res_file) or not os.path.getsize(salmon_res_file):
kkonganti@11: logging.error(
kkonganti@11: "The file\n"
kkonganti@11: + f"{salmon_res_file}\ndoes not exist or is empty!\n"
kkonganti@11: + "Did `salmon quant` fail?"
kkonganti@11: )
kkonganti@11: exit(1)
kkonganti@11:
kkonganti@11: with open(salmon_res_file, "r") as salmon_res_fh:
kkonganti@11: for line in salmon_res_fh.readlines():
kkonganti@11: if re.match(r"^Name.+", line):
kkonganti@11: continue
kkonganti@11: cols = line.strip().split("\t")
kkonganti@11: ref_acc = "_".join(cols[0].split("_")[:2])
kkonganti@11:
kkonganti@11: if ref_acc not in snp_clusters.keys():
kkonganti@11: snp_clusters[ref_acc]["snp_clust_id"] = ref_acc
kkonganti@11: snp_clusters[ref_acc]["pathdb_acc_id"] = ref_acc
kkonganti@11:
kkonganti@11: (
kkonganti@11: sample2salmon[sample_name]
kkonganti@11: .setdefault(acc2sero[cols[0]], [])
kkonganti@11: .append(int(round(float(cols[4]), 2)))
kkonganti@11: )
kkonganti@11: (
kkonganti@11: sample2salmon[sample_name]
kkonganti@11: .setdefault("snp_clust_ids", {})
kkonganti@11: .setdefault("".join(snp_clusters[ref_acc]["snp_clust_id"]), [])
kkonganti@11: .append("".join(snp_clusters[ref_acc]["pathdb_acc_id"]))
kkonganti@11: )
kkonganti@11: seen_sero[acc2sero[cols[0]]] = 1
kkonganti@11:
kkonganti@11: salmon_meta_json_read = json.load(open(salmon_meta_json, "r"))
kkonganti@11: (
kkonganti@11: sample2salmon[sample_name]
kkonganti@11: .setdefault("tot_reads", [])
kkonganti@11: .append(salmon_meta_json_read["num_processed"])
kkonganti@11: )
kkonganti@11:
kkonganti@11: with open(salmon_comb_res, "w") as salmon_comb_res_fh:
kkonganti@11:
kkonganti@11: # snp_clust_col_header = (
kkonganti@11: # "\tSNP Cluster(s) by Genome Hit\n" if show_snp_clust_col else "\n"
kkonganti@11: # )
kkonganti@11: snp_clust_col_header = (
kkonganti@11: "\tNCBI Pathogens Isolate Browser\n" if show_snp_clust_col else "\n"
kkonganti@11: )
kkonganti@11: serotypes = sorted(seen_sero.keys())
kkonganti@11: formatted_serotypes = [
kkonganti@11: re.sub(r"\,antigen_formula=", " | ", s)
kkonganti@11: for s in [re.sub(r"serotype=", "", s) for s in serotypes]
kkonganti@11: ]
kkonganti@11: salmon_comb_res_fh.write(
kkonganti@11: "Sample\t" + "\t".join(formatted_serotypes) + snp_clust_col_header
kkonganti@11: )
kkonganti@11: # sample_snp_relation = (
kkonganti@11: # ncbi_pathogens_base_url
kkonganti@11: # + pdg_release
kkonganti@11: # + "".join(snp_clusters[ref_acc]["snp_clust_id"])
kkonganti@11: # + "?accessions="
kkonganti@11: # )
kkonganti@11: if len(salmon_res_file_failed) == 1:
kkonganti@11: with (open("".join(salmon_res_file_failed), "r")) as no_calls_fh:
kkonganti@11: for line in no_calls_fh.readlines():
kkonganti@11: if line in ["\n", "\n\r", "\r"]:
kkonganti@11: continue
kkonganti@11: salmon_comb_res_fh.write(line.strip())
kkonganti@11: sal_yn[line.strip()] += 0
kkonganti@11: for serotype in serotypes:
kkonganti@11: salmon_comb_res_fh.write("\t-")
kkonganti@11: salmon_comb_res_fh.write(
kkonganti@11: "\t-\n"
kkonganti@11: ) if show_snp_clust_col else salmon_comb_res_fh.write("\n")
kkonganti@11: no_calls_fh.close()
kkonganti@11:
kkonganti@11: for sample, counts in sorted(sample2salmon.items()):
kkonganti@11: salmon_comb_res_fh.write(sample)
kkonganti@11: snp_cluster_res_col = list()
kkonganti@11:
kkonganti@11: for snp_clust_id in sample2salmon[sample]["snp_clust_ids"].keys():
kkonganti@11: # print(snp_clust_id)
kkonganti@11: # print(",".join(sample2salmon[sample]["snp_clust_ids"][snp_clust_id]))
kkonganti@11: # ppp.pprint(sample2salmon[sample]["snp_clust_ids"])
kkonganti@11: # ppp.pprint(sample2salmon[sample]["snp_clust_ids"][snp_clust_id])
kkonganti@11: # final_url_text = ",".join(
kkonganti@11: # sample2salmon[sample]["snp_clust_ids"][snp_clust_id]
kkonganti@11: # )
kkonganti@11: # final_url_text_to_show = snp_clust_id
kkonganti@11: # snp_cluster_res_col.append(
kkonganti@11: # "".join(
kkonganti@11: # [
kkonganti@11: # f'{snp_clust_id}',
kkonganti@11: # ]
kkonganti@11: # )
kkonganti@11: # )
kkonganti@11: # ppp.pprint(sample2salmon[sample])
kkonganti@11: for pathdbacc in sample2salmon[sample]["snp_clust_ids"][snp_clust_id]:
kkonganti@11: # final_url_text_to_show = " ".join(
kkonganti@11: # sample2salmon[sample]["snp_clust_ids"][snp_clust_id]
kkonganti@11: # )
kkonganti@11: sample_snp_relation = (
kkonganti@11: ncbi_pathogens_genome_base
kkonganti@11: if assm_pat.match(pathdbacc)
kkonganti@11: else ncbi_pathogens_base_url + "isolates/#"
kkonganti@11: )
kkonganti@11:
kkonganti@11: snp_cluster_res_col.append(
kkonganti@11: "".join(
kkonganti@11: [
kkonganti@11: f'{pathdbacc}',
kkonganti@11: ]
kkonganti@11: )
kkonganti@11: )
kkonganti@11:
kkonganti@11: per_serotype_counts = 0
kkonganti@11: for serotype in serotypes:
kkonganti@11:
kkonganti@11: if serotype in sample2salmon[sample].keys():
kkonganti@11: # ppp.pprint(counts)
kkonganti@11: sample_perc_mapped = round(
kkonganti@11: sum(counts[serotype]) / sum(counts["tot_reads"]) * 100, 2
kkonganti@11: )
kkonganti@11: salmon_comb_res_fh.write(
kkonganti@11: f"\t{sum(counts[serotype])} ({sample_perc_mapped}%)"
kkonganti@11: )
kkonganti@11: multiqc_salmon_counts[sample].setdefault(
kkonganti@11: re.match(r"^serotype=(.+?)\,antigen_formula.*", serotype)[1],
kkonganti@11: sum(counts[serotype]),
kkonganti@11: )
kkonganti@11: per_serotype_counts += sum(counts[serotype])
kkonganti@11: sal_yn[sample] += 1
kkonganti@11: else:
kkonganti@11: salmon_comb_res_fh.write(f"\t-")
kkonganti@11: sal_yn[sample] += 0
kkonganti@11:
kkonganti@11: multiqc_salmon_counts[sample].setdefault(
kkonganti@11: no_hit, sum(counts["tot_reads"]) - per_serotype_counts
kkonganti@11: )
kkonganti@11: snp_clust_col_val = (
kkonganti@11: f'\t{" ".join(snp_cluster_res_col)}\n' if show_snp_clust_col else "\n"
kkonganti@11: )
kkonganti@11: # ppp.pprint(multiqc_salmon_counts)
kkonganti@11: salmon_comb_res_fh.write(snp_clust_col_val)
kkonganti@11:
kkonganti@11: with open(bcs_sal_yn, "w") as bcs_sal_yn_fh:
kkonganti@11: bcs_sal_yn_fh.write("Sample\tSalmonella Presence\tNo. of Serotypes\n")
kkonganti@11: for sample in sal_yn.keys():
kkonganti@11: if sal_yn[sample] > 0:
kkonganti@11: bcs_sal_yn_fh.write(f"{sample}\tDetected\t{sal_yn[sample]}\n")
kkonganti@11: else:
kkonganti@11: bcs_sal_yn_fh.write(f"{sample}\tNot detected\t{sal_yn[sample]}\n")
kkonganti@11:
kkonganti@11: with open(cell_colors_yml_file, "w") as cell_colors_fh:
kkonganti@11: yaml.dump(cell_colors_yml, cell_colors_fh, default_flow_style=False)
kkonganti@11:
kkonganti@11: salmon_plot_json(salmon_comb_res_mqc, multiqc_salmon_counts, no_hit)
kkonganti@11:
kkonganti@11: salmon_comb_res_fh.close()
kkonganti@11: bcs_sal_yn_fh.close()
kkonganti@11: cell_colors_fh.close()
kkonganti@11:
kkonganti@11:
kkonganti@11: def salmon_plot_json(file: None, sample_salmon_counts: None, no_hit: None) -> None:
kkonganti@11: """
kkonganti@11: This method will take a dictionary of salmon counts per sample
kkonganti@11: and will dump a JSON that will be used by MultiQC.
kkonganti@11: """
kkonganti@11:
kkonganti@11: if file is None or sample_salmon_counts is None:
kkonganti@11: logging.error(
kkonganti@11: "Neither an output file to dump the JSON for MultiQC or the"
kkonganti@11: + "dictionary holding the salmon counts was not passed."
kkonganti@11: )
kkonganti@11:
kkonganti@11: # Credit: http://phrogz.net/tmp/24colors.html
kkonganti@11: # Will cycle through 20 distinct colors.
kkonganti@11: distinct_color_palette = [
kkonganti@11: "#FF0000",
kkonganti@11: "#FFFF00",
kkonganti@11: "#00EAFF",
kkonganti@11: "#AA00FF",
kkonganti@11: "#FF7F00",
kkonganti@11: "#BFFF00",
kkonganti@11: "#0095FF",
kkonganti@11: "#FF00AA",
kkonganti@11: "#FFD400",
kkonganti@11: "#6AFF00",
kkonganti@11: "#0040FF",
kkonganti@11: "#EDB9B9",
kkonganti@11: "#B9D7ED",
kkonganti@11: "#E7E9B9",
kkonganti@11: "#DCB9ED",
kkonganti@11: "#B9EDE0",
kkonganti@11: "#8F2323",
kkonganti@11: "#23628F",
kkonganti@11: "#8F6A23",
kkonganti@11: "#6B238F",
kkonganti@11: "#4F8F23",
kkonganti@11: ]
kkonganti@11:
kkonganti@11: # Credit: https://mokole.com/palette.html
kkonganti@11: # Will use this palette if we run out ouf
kkonganti@11: # 20 serotypes. More than 50 serotypes
kkonganti@11: # per run is probably rare but if not,
kkonganti@11: # will cycle through about 45.
kkonganti@11: distinct_color_palette2 = [
kkonganti@11: "#2F4F4F", # darkslategray
kkonganti@11: "#556B2F", # darkolivegreen
kkonganti@11: "#A0522D", # sienna
kkonganti@11: "#2E8B57", # seagreen
kkonganti@11: "#006400", # darkgreen
kkonganti@11: "#8B0000", # darkred
kkonganti@11: "#808000", # olive
kkonganti@11: "#BC8F8F", # rosybrown
kkonganti@11: "#663399", # rebeccapurple
kkonganti@11: "#B8860B", # darkgoldenrod
kkonganti@11: "#4682B4", # steelblue
kkonganti@11: "#000080", # navy
kkonganti@11: "#D2691E", # chocolate
kkonganti@11: "#9ACD32", # yellowgreen
kkonganti@11: "#20B2AA", # lightseagreen
kkonganti@11: "#CD5C5C", # indianred
kkonganti@11: "#8FBC8F", # darkseagreen
kkonganti@11: "#800080", # purple
kkonganti@11: "#B03060", # maroon3
kkonganti@11: "#FF8C00", # darkorange
kkonganti@11: "#FFD700", # gold
kkonganti@11: "#FFFF00", # yellow
kkonganti@11: "#DEB887", # burlywood
kkonganti@11: "#00FF00", # lime
kkonganti@11: "#BA55D3", # mediumorchid
kkonganti@11: "#00FA9A", # mediumspringgreen
kkonganti@11: "#4169E1", # royalblue
kkonganti@11: "#E9967A", # darksalmon
kkonganti@11: "#DC143C", # crimson
kkonganti@11: "#00FFFF", # aqua
kkonganti@11: "#F4A460", # sandybrown
kkonganti@11: "#9370DB", # mediumpurple
kkonganti@11: "#0000FF", # blue
kkonganti@11: "#ADFF2F", # greenyellow
kkonganti@11: "#FF6347", # tomato
kkonganti@11: "#D8BFD8", # thistle
kkonganti@11: "#FF00FF", # fuchsia
kkonganti@11: "#DB7093", # palevioletred
kkonganti@11: "#F0E68C", # khaki
kkonganti@11: "#6495ED", # cornflower
kkonganti@11: "#DDA0DD", # plum
kkonganti@11: "#EE82EE", # violet
kkonganti@11: "#7FFFD4", # aquamarine
kkonganti@11: "#FAFAD2", # lightgoldenrod
kkonganti@11: "#FF69B4", # hotpink
kkonganti@11: "#FFB6C1", # lightpink
kkonganti@11: ]
kkonganti@11:
kkonganti@11: no_hit_color = "#434348"
kkonganti@11: col_count = 0
kkonganti@11: serotypes = set()
kkonganti@11: salmon_counts = defaultdict(defaultdict)
kkonganti@11: salmon_counts["id"] = "BETTERCALLSAL_SALMON_COUNTS"
kkonganti@11: salmon_counts["section_name"] = "Salmon read counts"
kkonganti@11: salmon_counts["description"] = (
kkonganti@11: "This section shows the read counts from running salmon "
kkonganti@11: + "in --meta mode using SE, merged PE or concatenated PE reads against "
kkonganti@11: + "an on-the-fly salmon index generated from the genome hits "
kkonganti@11: + "of kma."
kkonganti@11: )
kkonganti@11: salmon_counts["plot_type"] = "bargraph"
kkonganti@11: salmon_counts["pconfig"]["id"] = "bettercallsal_salmon_counts_plot"
kkonganti@11: salmon_counts["pconfig"]["title"] = "Salmon: Read counts"
kkonganti@11: salmon_counts["pconfig"]["ylab"] = "Number of reads"
kkonganti@11: salmon_counts["pconfig"]["xDecimals"] = "false"
kkonganti@11: salmon_counts["pconfig"]["cpswitch_counts_label"] = "Number of reads (Counts)"
kkonganti@11: salmon_counts["pconfig"]["cpswitch_percent_label"] = "Number of reads (Percentages)"
kkonganti@11:
kkonganti@11: for sample in sorted(sample_salmon_counts.keys()):
kkonganti@11: serotypes.update(list(sample_salmon_counts[sample].keys()))
kkonganti@11: salmon_counts["data"][sample] = sample_salmon_counts[sample]
kkonganti@11:
kkonganti@11: if len(serotypes) > len(distinct_color_palette):
kkonganti@11: distinct_color_palette = distinct_color_palette2
kkonganti@11:
kkonganti@11: for serotype in sorted(serotypes):
kkonganti@11: if serotype == no_hit:
kkonganti@11: continue
kkonganti@11: if col_count == len(distinct_color_palette) - 1:
kkonganti@11: col_count = 0
kkonganti@11:
kkonganti@11: col_count += 1
kkonganti@11: salmon_counts["categories"][serotype] = {"color": distinct_color_palette[col_count]}
kkonganti@11:
kkonganti@11: salmon_counts["categories"][no_hit] = {"color": no_hit_color}
kkonganti@11: json.dump(salmon_counts, open(file, "w"))
kkonganti@11:
kkonganti@11:
kkonganti@11: if __name__ == "__main__":
kkonganti@11: main()