Allow gene predictions to overlap on opposite strands

Author: lisavader

src/geneml/main.py

@@ -239,6 +239,7 @@ def parse_args(argv=None):
     advanced.add_argument('--gene-candidates', type=int, default=100, help="Maximum number of gene candidates to consider (default: %(default)s).")
     advanced.add_argument('--gene-range-time-limit', type=int, default=None, help="Time limit for each gene prediction, in seconds (default: %(default)s)")
     advanced.add_argument('--hardmask-repeats-min-size', type=int, default=None, help="Minimum size of softmasked repeat stretches to hardmask (default: %(default)s).")
+    advanced.add_argument('--allow_opposite_strand_overlaps', type=bool, default=True, help="Predict overlapping genes on opposite strands (default: %(default)s).")
 
     args = parser.parse_args(argv)
     check_args(parser, args)

src/geneml/params.py

@@ -27,6 +27,7 @@ class Params(namedtuple('Params', (
     'contigs_filter', 'output_segs', 'output_genes', 'output_proteins',
     'num_cores', 'debug', 'verbose', 'basepath', 'inpath', 'outpath',
     'hardmask_repeats_min_size', 'single_recurse_max_num_ops', 'recurse_region_max_num_ops',
+    'allow_opposite_strand_overlaps',
 ))):
     def to_json(self, **kwargs):
         return json.dumps(self._asdict(), cls=EnhancedJSONEncoder, **kwargs)
@@ -69,6 +70,8 @@ def build_params_namedtuple(args: Namespace) -> Params:
         'min_gene_score': args.min_gene_score,
         'gene_candidates': args.gene_candidates,
         'hardmask_repeats_min_size': args.hardmask_repeats_min_size,
+        'allow_opposite_strand_overlaps': args.allow_opposite_strand_overlaps,
+
         'single_recurse_max_num_ops': 100000,
         'recurse_region_max_num_ops': 200000,
     }

src/geneml/produce_genes.py

@@ -90,7 +90,8 @@ def create_transcripts(scored_gene_calls: list[tuple[float, list[GeneEvent]]],
     return transcripts
 
 
-def filter_transcripts(transcripts: list[Transcript], min_score: float):
+def filter_transcripts(transcripts: list[Transcript], min_score: float,
+                       allow_opposite_strand_overlaps: bool) -> list[Transcript]:
     """ This takes potential gene calls on both forward strand and reverse complement, filters them to remove
     overlapping calls, starting from highest scores first, and returns the best ones"""
 
@@ -103,7 +104,7 @@ def filter_transcripts(transcripts: list[Transcript], min_score: float):
     valid_transcripts.sort(key=lambda x: x.start)
     non_overlapping = [valid_transcripts[0]]
     for t in valid_transcripts[1:]:
-        if non_overlapping[-1].overlaps_with(t, ignore_strand=True):
+        if non_overlapping[-1].overlaps_with(t, ignore_strand=not allow_opposite_strand_overlaps):
             if t.score > non_overlapping[-1].score:
                 non_overlapping[-1] = t
         else:
@@ -141,7 +142,7 @@ def build_transcripts(preds: Optional[np.ndarray], rc_preds: Optional[np.ndarray
 
     logger.info('Selecting best gene calls for %s', contig_id)
     filtered_transcripts = filter_transcripts(
-        transcripts, params.min_gene_score)
+        transcripts, params.min_gene_score, params.allow_opposite_strand_overlaps)
 
     return filtered_transcripts