Mercurial > repos > kkonganti > hfp_nowayout

annotate 0.5.0/bin/gen_otf_genome.py @ 0:97cd2f532efe

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload
author kkonganti
date Mon, 31 Mar 2025 14:50:40 -0400
parents
children
rev   line source
kkonganti@0 1 #!/usr/bin/env python3
kkonganti@0 2
kkonganti@0 3 # Kranti Konganti
kkonganti@0 4
kkonganti@0 5 import argparse
kkonganti@0 6 import glob
kkonganti@0 7 import gzip
kkonganti@0 8 import inspect
kkonganti@0 9 import logging
kkonganti@0 10 import os
kkonganti@0 11 import pprint
kkonganti@0 12 import re
kkonganti@0 13
kkonganti@0 14 # Set logging.
kkonganti@0 15 logging.basicConfig(
kkonganti@0 16 format="\n"
kkonganti@0 17 + "=" * 55
kkonganti@0 18 + "\n%(asctime)s - %(levelname)s\n"
kkonganti@0 19 + "=" * 55
kkonganti@0 20 + "\n%(message)s\n\n",
kkonganti@0 21 level=logging.DEBUG,
kkonganti@0 22 )
kkonganti@0 23
kkonganti@0 24 # Debug print.
kkonganti@0 25 ppp = pprint.PrettyPrinter(width=50, indent=4)
kkonganti@0 26
kkonganti@0 27
kkonganti@0 28 # Multiple inheritence for pretty printing of help text.
kkonganti@0 29 class MultiArgFormatClasses(
kkonganti@0 30 argparse.RawTextHelpFormatter, argparse.ArgumentDefaultsHelpFormatter
kkonganti@0 31 ):
kkonganti@0 32 pass
kkonganti@0 33
kkonganti@0 34
kkonganti@0 35 def main() -> None:
kkonganti@0 36 """
kkonganti@0 37 This script works only in the context of a Nextflow workflow.
kkonganti@0 38 It takes:
kkonganti@0 39 1. A text file containing accessions or FASTA IDs, one per line and
kkonganti@0 40 then,
kkonganti@0 41 2. Optionally, searches for a genome FASTA file in gzipped format in specified
kkonganti@0 42 search path, where the prefix of the filename is the accession or
kkonganti@0 43 FASTA ID from 1. and then, creates a new concatenated gzipped genome FASTA
kkonganti@0 44 file with all the genomes in the text file from 1.
kkonganti@0 45 3. Creates a new FASTQ file with reads aligned to the accessions in the text
kkonganti@0 46 file from 1.
kkonganti@0 47 """
kkonganti@0 48
kkonganti@0 49 prog_name = os.path.basename(inspect.stack()[0].filename)
kkonganti@0 50
kkonganti@0 51 parser = argparse.ArgumentParser(
kkonganti@0 52 prog=prog_name, description=main.__doc__, formatter_class=MultiArgFormatClasses
kkonganti@0 53 )
kkonganti@0 54
kkonganti@0 55 required = parser.add_argument_group("required arguments")
kkonganti@0 56
kkonganti@0 57 required.add_argument(
kkonganti@0 58 "-txt",
kkonganti@0 59 dest="accs_txt",
kkonganti@0 60 default=False,
kkonganti@0 61 required=True,
kkonganti@0 62 help="Absolute UNIX path to .txt file containing accessions\n"
kkonganti@0 63 + "FASTA IDs, one per line.",
kkonganti@0 64 )
kkonganti@0 65 required.add_argument(
kkonganti@0 66 "-op",
kkonganti@0 67 dest="out_prefix",
kkonganti@0 68 default="CATTED_GENOMES",
kkonganti@0 69 required=True,
kkonganti@0 70 help="Set the output file prefix for .fna.gz and .txt\n" + "files.",
kkonganti@0 71 )
kkonganti@0 72 parser.add_argument(
kkonganti@0 73 "-gd",
kkonganti@0 74 dest="genomes_dir",
kkonganti@0 75 default=False,
kkonganti@0 76 required=False,
kkonganti@0 77 help="Absolute UNIX path to a directory containing\n"
kkonganti@0 78 + "gzipped genome FASTA files or a file.\n",
kkonganti@0 79 )
kkonganti@0 80 parser.add_argument(
kkonganti@0 81 "-gds",
kkonganti@0 82 dest="genomes_dir_suffix",
kkonganti@0 83 default="_scaffolded_genomic.fna.gz",
kkonganti@0 84 required=False,
kkonganti@0 85 help="Genome FASTA file suffix to search for\nin the directory mentioned using\n-gd.",
kkonganti@0 86 )
kkonganti@0 87 parser.add_argument(
kkonganti@0 88 "-query",
kkonganti@0 89 dest="id_is_query",
kkonganti@0 90 default=False,
kkonganti@0 91 action="store_true",
kkonganti@0 92 required=False,
kkonganti@0 93 help="In the produced FASTQ file, should the FASTA ID should be of KMA query ID\n"
kkonganti@0 94 + "or template ID.",
kkonganti@0 95 )
kkonganti@0 96 parser.add_argument(
kkonganti@0 97 "-txts",
kkonganti@0 98 dest="accs_suffix",
kkonganti@0 99 default="_template_hits.txt",
kkonganti@0 100 required=False,
kkonganti@0 101 help="The suffix of the file supplied with -txt option. It is assumed that the\n"
kkonganti@0 102 + "sample name is present in the file supplied with -txt option and the suffix\n"
kkonganti@0 103 + "will be stripped and stored in a file that logs samples which have no hits.",
kkonganti@0 104 )
kkonganti@0 105 parser.add_argument(
kkonganti@0 106 "-frag_delim",
kkonganti@0 107 dest="frag_delim",
kkonganti@0 108 default="\t",
kkonganti@0 109 required=False,
kkonganti@0 110 help="The delimitor by which the fields are separated in *_frag.gz file.",
kkonganti@0 111 )
kkonganti@0 112
kkonganti@0 113 args = parser.parse_args()
kkonganti@0 114 accs_txt = args.accs_txt
kkonganti@0 115 genomes_dir = args.genomes_dir
kkonganti@0 116 genomes_dir_suffix = args.genomes_dir_suffix
kkonganti@0 117 id_is_query = args.id_is_query
kkonganti@0 118 out_prefix = args.out_prefix
kkonganti@0 119 accs_suffix = args.accs_suffix
kkonganti@0 120 frag_delim = args.frag_delim
kkonganti@0 121 accs_seen = dict()
kkonganti@0 122 cat_genomes_gz = os.path.join(os.getcwd(), out_prefix + "_" + genomes_dir_suffix)
kkonganti@0 123 cat_genomes_gz = re.sub("__", "_", str(cat_genomes_gz))
kkonganti@0 124 frags_gz = os.path.join(os.getcwd(), out_prefix + ".frag.gz")
kkonganti@0 125 cat_reads_gz = os.path.join(os.getcwd(), out_prefix + "_aln_reads.fna.gz")
kkonganti@0 126 cat_reads_gz = re.sub("__", "_", cat_reads_gz)
kkonganti@0 127
kkonganti@0 128 if (
kkonganti@0 129 accs_txt
kkonganti@0 130 and os.path.exists(cat_genomes_gz)
kkonganti@0 131 and os.path.getsize(cat_genomes_gz) > 0
kkonganti@0 132 ):
kkonganti@0 133 logging.error(
kkonganti@0 134 "A concatenated genome FASTA file,\n"
kkonganti@0 135 + f"{os.path.basename(cat_genomes_gz)} already exists in:\n"
kkonganti@0 136 + f"{os.getcwd()}\n"
kkonganti@0 137 + "Please remove or move it as we will not "
kkonganti@0 138 + "overwrite it."
kkonganti@0 139 )
kkonganti@0 140 exit(1)
kkonganti@0 141
kkonganti@0 142 if accs_txt and (not os.path.exists(accs_txt) or not os.path.getsize(accs_txt) > 0):
kkonganti@0 143 logging.error("File,\n" + f"{accs_txt}\ndoes not exist " + "or is empty!")
kkonganti@0 144 failed_sample_name = re.sub(accs_suffix, "", os.path.basename(accs_txt))
kkonganti@0 145 with open(
kkonganti@0 146 os.path.join(os.getcwd(), "_".join([out_prefix, "FAILED.txt"])), "w"
kkonganti@0 147 ) as failed_sample_fh:
kkonganti@0 148 failed_sample_fh.write(f"{failed_sample_name}\n")
kkonganti@0 149 failed_sample_fh.close()
kkonganti@0 150 exit(0)
kkonganti@0 151
kkonganti@0 152 # ppp.pprint(mash_hits)
kkonganti@0 153 empty_lines = 0
kkonganti@0 154 empty_lines_msg = ""
kkonganti@0 155
kkonganti@0 156 with open(accs_txt, "r") as accs_txt_fh:
kkonganti@0 157 for line in accs_txt_fh:
kkonganti@0 158 if line in ["\n", "\n\r"]:
kkonganti@0 159 empty_lines += 1
kkonganti@0 160 continue
kkonganti@0 161 else:
kkonganti@0 162 line = line.strip()
kkonganti@0 163
kkonganti@0 164 if line in accs_seen.keys():
kkonganti@0 165 continue
kkonganti@0 166 else:
kkonganti@0 167 accs_seen[line] = 1
kkonganti@0 168 accs_txt_fh.close()
kkonganti@0 169
kkonganti@0 170 if genomes_dir:
kkonganti@0 171 if not os.path.isdir(genomes_dir):
kkonganti@0 172 logging.error("UNIX path\n" + f"{genomes_dir}\n" + "does not exist!")
kkonganti@0 173 exit(1)
kkonganti@0 174 if len(glob.glob(os.path.join(genomes_dir, "*" + genomes_dir_suffix))) <= 0:
kkonganti@0 175 logging.error(
kkonganti@0 176 "Genomes directory"
kkonganti@0 177 + f"{genomes_dir}"
kkonganti@0 178 + "\ndoes not seem to have any\n"
kkonganti@0 179 + f"files ending with suffix: {genomes_dir_suffix}"
kkonganti@0 180 )
kkonganti@0 181 exit(1)
kkonganti@0 182
kkonganti@0 183 with open(cat_genomes_gz, "wb") as genomes_out_gz:
kkonganti@0 184 for line in accs_seen.keys():
kkonganti@0 185 genome_file = os.path.join(genomes_dir, line + genomes_dir_suffix)
kkonganti@0 186
kkonganti@0 187 if not os.path.exists(genome_file) or os.path.getsize(genome_file) <= 0:
kkonganti@0 188 logging.error(
kkonganti@0 189 f"Genome file {os.path.basename(genome_file)} does not\n"
kkonganti@0 190 + "exits or is empty!"
kkonganti@0 191 )
kkonganti@0 192 exit(1)
kkonganti@0 193 else:
kkonganti@0 194 with open(genome_file, "rb") as genome_file_h:
kkonganti@0 195 genomes_out_gz.writelines(genome_file_h.readlines())
kkonganti@0 196 genome_file_h.close()
kkonganti@0 197 genomes_out_gz.close()
kkonganti@0 198
kkonganti@0 199 if (
kkonganti@0 200 len(accs_seen.keys()) > 0
kkonganti@0 201 and os.path.exists(frags_gz)
kkonganti@0 202 and os.path.getsize(frags_gz) > 0
kkonganti@0 203 ):
kkonganti@0 204 with gzip.open(
kkonganti@0 205 cat_reads_gz, "wt", encoding="utf-8", compresslevel=6
kkonganti@0 206 ) as cat_reads_gz_fh:
kkonganti@0 207 with gzip.open(frags_gz, "rb", compresslevel=6) as fragz_gz_fh:
kkonganti@0 208 fasta_id = 7 if id_is_query else 6
kkonganti@0 209 for frag_line in fragz_gz_fh:
kkonganti@0 210 frag_lines = frag_line.decode("utf-8").strip().split(frag_delim)
kkonganti@0 211 # Per KMA specification, 6=template, 7=query, 1=read
kkonganti@0 212 cat_reads_gz_fh.write(f">{frag_lines[fasta_id]}\n{frag_lines[0]}\n")
kkonganti@0 213 fragz_gz_fh.close()
kkonganti@0 214 cat_reads_gz_fh.close()
kkonganti@0 215
kkonganti@0 216 if empty_lines > 0:
kkonganti@0 217 empty_lines_msg = f"Skipped {empty_lines} empty line(s).\n"
kkonganti@0 218
kkonganti@0 219 logging.info(
kkonganti@0 220 empty_lines_msg
kkonganti@0 221 + f"File {os.path.basename(cat_genomes_gz)}\n"
kkonganti@0 222 + f"written in:\n{os.getcwd()}\nDone! Bye!"
kkonganti@0 223 )
kkonganti@0 224
kkonganti@0 225
kkonganti@0 226 if __name__ == "__main__":
kkonganti@0 227 main()