feat(IPVC-2276): skip accessions from gff that are not in the tx_info file, log missing acs to file

sbin/exonset-to-seqinfo

@@ -1,20 +1,20 @@
 #!/usr/bin/env python
 
-"""write SeqInfo files from fasta"""
+"""Writes SeqInfo files from fasta to stdout. Writes transcript accessions that are missing from
+SeqRepo to a file called exonset_missing_accessions.txt in the same directory as the input files."""
 
 import argparse
-import configparser as ConfigParser
 import gzip
 import itertools
-import logging
 import logging.config
-import pkg_resources
+import os
 import re
 import sys
 
-from bioutils.digests import seq_md5
+import configparser as ConfigParser
+import pkg_resources
 from biocommons.seqrepo import SeqRepo
-# from multifastadb import MultiFastaDB
+from bioutils.digests import seq_md5
 
 from uta.formats.exonset import ExonSetReader
 from uta.formats.seqinfo import SeqInfo, SeqInfoWriter
@@ -49,6 +49,8 @@ if __name__ == "__main__":
     ac_re = re.compile("[NX][CGMPR]_")
 
     opts = parse_args(sys.argv[1:])
+    input_dir = os.path.dirname(opts.FILES[0])
+    missing_accessions_file_path = os.path.join(input_dir, "exonset_missing_accessions.txt")
 
     cf = ConfigParser.ConfigParser()
     for conf_fn in opts.conf:
@@ -72,25 +74,27 @@ if __name__ == "__main__":
     # this is just a fancy way to make a set of all tx_ac and alt_ac accessions
     acs = sorted(
         set(itertools.chain.from_iterable((es.tx_ac, es.alt_ac) for es in esr)))
-
-    acs_not_found = set()
-    for ac in acs:
-        try:
-            seq = str(sr[ac])
-        except KeyError:
-            logging.warning("Sequence not found: " + ac)
-            acs_not_found.update([ac])
-            continue
-
-        si = SeqInfo(
-            ac=ac,
-            descr=None,
-            len=len(seq),
-            md5=seq_md5(seq),
-            origin=opts.origin,
-            seq=None,
-        )
-        siw.write(si)
-
-    if acs_not_found:
-        raise RuntimeError("Sequences for {} accessions not found".format(len(acs_not_found)))
+    missing_acs_count = 0
+    with open(missing_accessions_file_path, "w") as f:
+        f.write("tx_ac\n")
+        for ac in acs:
+            try:
+                seq = str(sr[ac])
+            except KeyError:
+                logging.warning("Sequence not found: " + ac)
+                missing_acs_count += 1
+                f.write(ac + "\n")
+                continue
+
+            si = SeqInfo(
+                ac=ac,
+                descr=None,
+                len=len(seq),
+                md5=seq_md5(seq),
+                origin=opts.origin,
+                seq=None,
+            )
+            siw.write(si)
+
+    # Log the number of missing accessions written to the file
+    logging.info(f"{missing_acs_count} accessions were missing from Seqrepo. See {missing_accessions_file_path}.")

sbin/ncbi_parse_genomic_gff.py

@@ -115,7 +115,7 @@ def _get_exon_number_from_id(alignment_id: str) -> int:
     return int(alignment_id.split("-")[-1])
 
 
-def parse_gff_file(file_paths: List[str]) -> dict[str, List[GFFRecord]]:
+def parse_gff_files(file_paths: List[str]) -> dict[str, List[GFFRecord]]:
     tx_data = defaultdict(list)
     for file_path in file_paths:
         with open_file(file_path) as f:
@@ -152,7 +152,7 @@ def get_zero_based_exon_ranges(transcript_exons: List[GFFRecord]) -> str:
     gff_files = args.gff_files
     esw = ExonSetWriter(sys.stdout)
 
-    transcript_alignments = parse_gff_file(gff_files)
+    transcript_alignments = parse_gff_files(gff_files)
     logger.info(
         f"read {len(transcript_alignments)} transcript alignments from file(s): {', '.join(gff_files)}"
     )

tests/test_ncbi_parse_genomic_gff.py

@@ -7,7 +7,7 @@
 from sbin.ncbi_parse_genomic_gff import (
     get_zero_based_exon_ranges,
     GFFRecord,
-    parse_gff_file,
+    parse_gff_files,
     parse_gff_record,
 )
 
@@ -160,7 +160,7 @@ def test_parse_gff_record_raises_unparseable_id(self):
     def test_parse_gff_file(self):
         # Test parsing the entire uncompressed GFF file
         expected_result = {"rna-NR_046018.2": self.gff_records}
-        parsed_result = parse_gff_file([self.temp_gff.name])
+        parsed_result = parse_gff_files([self.temp_gff.name])
         self.assertEqual(parsed_result, expected_result)
 
     def test_parse_gff_file_accepts_gzipped_files(self):
@@ -171,7 +171,7 @@ def test_parse_gff_file_accepts_gzipped_files(self):
 
         # Test parsing the gzipped GFF file
         expected_result = {"rna-NR_046018.2": self.gff_records}
-        parsed_result = parse_gff_file([self.temp_gff.name + ".gz"])
+        parsed_result = parse_gff_files([self.temp_gff.name + ".gz"])
         self.assertEqual(parsed_result, expected_result)
 
     def test_get_zero_based_exon_ranges(self):