Mercurial > repos > kkonganti > cfsan_bettercallsal

annotate 0.5.0/bin/gen_otf_genome.py @ 2:4678c2cd1c9e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload"
author kkonganti
date Tue, 06 Jun 2023 09:46:25 -0400
parents 365849f031fd
children
rev   line source
kkonganti@1 1 #!/usr/bin/env python3
kkonganti@1 2
kkonganti@1 3 # Kranti Konganti
kkonganti@1 4
kkonganti@1 5 import os
kkonganti@1 6 import argparse
kkonganti@1 7 import inspect
kkonganti@1 8 import logging
kkonganti@1 9 import re
kkonganti@1 10 import pprint
kkonganti@1 11 import glob
kkonganti@1 12
kkonganti@1 13 # Set logging.
kkonganti@1 14 logging.basicConfig(
kkonganti@1 15 format="\n" + "=" * 55 + "\n%(asctime)s - %(levelname)s\n" + "=" * 55 + "\n%(message)s\n\n",
kkonganti@1 16 level=logging.DEBUG,
kkonganti@1 17 )
kkonganti@1 18
kkonganti@1 19 # Debug print.
kkonganti@1 20 ppp = pprint.PrettyPrinter(width=50, indent=4)
kkonganti@1 21
kkonganti@1 22 # Multiple inheritence for pretty printing of help text.
kkonganti@1 23 class MultiArgFormatClasses(argparse.RawTextHelpFormatter, argparse.ArgumentDefaultsHelpFormatter):
kkonganti@1 24 pass
kkonganti@1 25
kkonganti@1 26
kkonganti@1 27 def main() -> None:
kkonganti@1 28 """
kkonganti@1 29 This script works only in the context of `bettercallsal` Nextflow workflow.
kkonganti@1 30 It takes:
kkonganti@1 31 1. A text file containing accessions or FASTA IDs, one per line and
kkonganti@1 32 then,
kkonganti@1 33 2. Searches for a genome FASTA file in gzipped format in specified
kkonganti@1 34 search path, where the prefix of the filename is the accession or
kkonganti@1 35 FASTA ID from 1. and then,
kkonganti@1 36 creates a new concatenated gzipped genome FASTA file with all the genomes
kkonganti@1 37 in the text file from 1.
kkonganti@1 38 """
kkonganti@1 39
kkonganti@1 40 prog_name = os.path.basename(inspect.stack()[0].filename)
kkonganti@1 41
kkonganti@1 42 parser = argparse.ArgumentParser(
kkonganti@1 43 prog=prog_name, description=main.__doc__, formatter_class=MultiArgFormatClasses
kkonganti@1 44 )
kkonganti@1 45
kkonganti@1 46 required = parser.add_argument_group("required arguments")
kkonganti@1 47
kkonganti@1 48 required.add_argument(
kkonganti@1 49 "-txt",
kkonganti@1 50 dest="accs_txt",
kkonganti@1 51 default=False,
kkonganti@1 52 required=True,
kkonganti@1 53 help="Absolute UNIX path to .txt file containing accessions\n" + "FASTA IDs, one per line.",
kkonganti@1 54 )
kkonganti@1 55 required.add_argument(
kkonganti@1 56 "-gd",
kkonganti@1 57 dest="genomes_dir",
kkonganti@1 58 default=False,
kkonganti@1 59 required=True,
kkonganti@1 60 help="Absolute UNIX path to a directory containing\n"
kkonganti@1 61 + "gzipped genome FASTA files.\n"
kkonganti@1 62 + "Required if -m is on.",
kkonganti@1 63 )
kkonganti@1 64 parser.add_argument(
kkonganti@1 65 "-gds",
kkonganti@1 66 dest="genomes_dir_suffix",
kkonganti@1 67 default="_scaffolded_genomic.fna.gz",
kkonganti@1 68 required=False,
kkonganti@1 69 help="Genome FASTA file suffix to search for\nin the directory mentioned using\n-gd.",
kkonganti@1 70 )
kkonganti@1 71 parser.add_argument(
kkonganti@1 72 "-op",
kkonganti@1 73 dest="out_prefix",
kkonganti@1 74 default="CATTED_GENOMES",
kkonganti@1 75 help="Set the output file prefix for .fna.gz and .txt\n" + "files.",
kkonganti@1 76 )
kkonganti@1 77 parser.add_argument(
kkonganti@1 78 "-txts",
kkonganti@1 79 dest="accs_suffix",
kkonganti@1 80 default="_template_hits.txt",
kkonganti@1 81 required=False,
kkonganti@1 82 help="The suffix of the file supplied with -txt option. It is assumed that the\n"
kkonganti@1 83 + "sample name is present in the file supplied with -txt option and the suffix\n"
kkonganti@1 84 + "will be stripped and stored in a file that logs samples which have no hits.",
kkonganti@1 85 )
kkonganti@1 86
kkonganti@1 87 args = parser.parse_args()
kkonganti@1 88 accs_txt = args.accs_txt
kkonganti@1 89 genomes_dir = args.genomes_dir
kkonganti@1 90 genomes_dir_suffix = args.genomes_dir_suffix
kkonganti@1 91 out_prefix = args.out_prefix
kkonganti@1 92 accs_suffix = args.accs_suffix
kkonganti@1 93 accs_seen = dict()
kkonganti@1 94 cat_genomes_gz = os.path.join(os.getcwd(), out_prefix + "_" + genomes_dir_suffix)
kkonganti@1 95 cat_genomes_gz = re.sub("__", "_", str(cat_genomes_gz))
kkonganti@1 96
kkonganti@1 97 if accs_txt and os.path.exists(cat_genomes_gz) and os.path.getsize(cat_genomes_gz) > 0:
kkonganti@1 98 logging.error(
kkonganti@1 99 "A concatenated genome FASTA file,\n"
kkonganti@1 100 + f"{os.path.basename(cat_genomes_gz)} already exists in:\n"
kkonganti@1 101 + f"{os.getcwd()}\n"
kkonganti@1 102 + "Please remove or move it as we will not "
kkonganti@1 103 + "overwrite it."
kkonganti@1 104 )
kkonganti@1 105 exit(1)
kkonganti@1 106
kkonganti@1 107 if accs_txt and (not os.path.exists(accs_txt) or not os.path.getsize(accs_txt) > 0):
kkonganti@1 108 logging.error("File,\n" + f"{accs_txt}\ndoes not exist " + "or is empty!")
kkonganti@1 109 failed_sample_name = re.sub(accs_suffix, "", os.path.basename(accs_txt))
kkonganti@1 110 with open(
kkonganti@1 111 os.path.join(os.getcwd(), "_".join([out_prefix, "FAILED.txt"])), "w"
kkonganti@1 112 ) as failed_sample_fh:
kkonganti@1 113 failed_sample_fh.write(f"{failed_sample_name}\n")
kkonganti@1 114 failed_sample_fh.close()
kkonganti@1 115 exit(0)
kkonganti@1 116
kkonganti@1 117 if genomes_dir:
kkonganti@1 118 if not os.path.isdir(genomes_dir):
kkonganti@1 119 logging.error("UNIX path\n" + f"{genomes_dir}\n" + "does not exist!")
kkonganti@1 120 exit(1)
kkonganti@1 121 if len(glob.glob(os.path.join(genomes_dir, "*" + genomes_dir_suffix))) <= 0:
kkonganti@1 122 logging.error(
kkonganti@1 123 "Genomes directory"
kkonganti@1 124 + f"{genomes_dir}"
kkonganti@1 125 + "\ndoes not seem to have any\n"
kkonganti@1 126 + f"files ending with suffix: {genomes_dir_suffix}"
kkonganti@1 127 )
kkonganti@1 128 exit(1)
kkonganti@1 129
kkonganti@1 130 # ppp.pprint(mash_hits)
kkonganti@1 131 empty_lines = 0
kkonganti@1 132 empty_lines_msg = ""
kkonganti@1 133 with open(cat_genomes_gz, "wb") as genomes_out_gz:
kkonganti@1 134 with open(accs_txt, "r") as accs_txt_fh:
kkonganti@1 135 for line in accs_txt_fh:
kkonganti@1 136 if line in ["\n", "\n\r"]:
kkonganti@1 137 empty_lines += 1
kkonganti@1 138 continue
kkonganti@1 139 else:
kkonganti@1 140 line = line.strip()
kkonganti@1 141
kkonganti@1 142 if line in accs_seen.keys():
kkonganti@1 143 continue
kkonganti@1 144 else:
kkonganti@1 145 accs_seen[line] = 1
kkonganti@1 146
kkonganti@1 147 genome_file = os.path.join(genomes_dir, line + genomes_dir_suffix)
kkonganti@1 148
kkonganti@1 149 if not os.path.exists(genome_file) or os.path.getsize(genome_file) <= 0:
kkonganti@1 150 logging.error(
kkonganti@1 151 f"Genome file {os.path.basename(genome_file)} does not\n"
kkonganti@1 152 + "exits or is empty!"
kkonganti@1 153 )
kkonganti@1 154 exit(1)
kkonganti@1 155 else:
kkonganti@1 156 with open(genome_file, "rb") as genome_file_h:
kkonganti@1 157 genomes_out_gz.writelines(genome_file_h.readlines())
kkonganti@1 158 genome_file_h.close()
kkonganti@1 159 accs_txt_fh.close()
kkonganti@1 160 genomes_out_gz.close()
kkonganti@1 161
kkonganti@1 162 if empty_lines > 0:
kkonganti@1 163 empty_lines_msg = f"Skipped {empty_lines} empty line(s).\n"
kkonganti@1 164
kkonganti@1 165 logging.info(
kkonganti@1 166 empty_lines_msg
kkonganti@1 167 + f"File {os.path.basename(cat_genomes_gz)}\n"
kkonganti@1 168 + f"written in:\n{os.getcwd()}\nDone! Bye!"
kkonganti@1 169 )
kkonganti@1 170 exit(0)
kkonganti@1 171
kkonganti@1 172
kkonganti@1 173 if __name__ == "__main__":
kkonganti@1 174 main()