Transcriptome polish

.flake8

@@ -1,19 +1,11 @@
 [flake8]
-# E722: do not use bare except; maybe fix later
+# E302: expected 2 blank lines
+# E741: ambiguous variable name
+# F523: .format unused named arguments
+# F524: .format missing named arguments
+# W503: line break before binary operator (PEP 8 preferred style)
+# W504: line break after binary operator
 ignore = E302,E741,F523,F524,W503,W504
 max-complexity = 10
 max-line-length = 256
-
-per-file-ignores =
-    # imported but unused
-    src/flair/count_sam_transcripts.py: E502
-    # fixing W504 causes E502
-    src/flair/filter_collapsed_isoforms.py,src/flair/collapse_isoforms_precise.py,src/flair/call_diffsplice_events.py: W504
-    # openblas setting must be in middle of import section
-    deFLAIR.py,es_as_inc_excl_to_counts.py,flair.py,runDE.py,runDS.py,runDU.py,samJuncs.py: E402
-    src/flair/ssPrep.py: F401 
-    # hard to fix in-loop variable definitions (maybe later)
-    src/flair/identify_gene_isoform.py: F821
-    # FIXME: make these less complex
-    src/flair/flair_transcriptome.py: C901
-    src/flair/filter_transcriptome_align.py: C901
+exclude = .claude.bak,src/flair/pycbio

.github/workflows/lint.yml

@@ -0,0 +1,12 @@
+name: Lint
+on: [push, pull_request]
+jobs:
+  flake8:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
+        with:
+          python-version: '3.x'
+      - run: pip install flake8
+      - run: make flake8

CHANGELOG.md

@@ -1,6 +1,9 @@
 # Major user-visible changes
 
 ## [v3.x.xx]
+* General
+  * fixed problems with some diffsplice_fishers_exact, and other
+    auxiliary console script installs
 * Incompatibles
   * Removed flair correct and collapse modules, the functionality is replaced
     by flair transcriptome.

Makefile

@@ -9,6 +9,11 @@ include ${root}/defs.mk
 #   make test-installed
 #
 
+# get files to include in flake8
+PYPROGS = $(shell file -F $$'\t' test/bin/* | awk '/Python script/{print $$1}')
+
+FLAKE8_CHECK = src/flair/*.py test/*.py ${PYPROGS}
+
 default:
 
 doc:
@@ -34,7 +39,6 @@ test-base-installed:
 #   see .flake8 for configuration
 #   due to to gradual cleanup of code, flake8.mk is the list of files to check
 ##
-include flake8.mk
 lint: flake8
 pycbio-lint: pycbio-flake8
 

docs/source/modules.rst

@@ -70,7 +70,8 @@ Optional arguments
   --stringent           [HIGHLY RECOMMENDED] specify if all supporting reads need to be full-length 
                         (spanning 25 bp of the first and last exons)
   --check_splice        [HIGHLY RECOMMENDED] enforce coverage of 4 out of 6 bp around each splice site 
-                        and no insertions greater than 3 bp at the splice site DON'T USE WITH DATA WITH HIGH ERROR RATES (old direct-RNA)
+                        and no insertions greater than 3 bp at the splice site. DON'T USE WITH DATA WITH ERROR RATES of 5% or more, such
+                        as older direct-RNA (anything before using the RNA-specific flow-cell with the Dorado basecaller).
   -w --end_window
                         window size for comparing TSS/TES (100)
   --noaligntoannot      related to old annotation_reliant, now specify if you don't want an initial alignment 

pyproject.toml

@@ -18,13 +18,12 @@ dependencies = [
     "pysam (>=0.23.0,<0.24.0)",
     "pipettor (>=1.1.0,<2.0.0)",
     "ncls (>=0.0.70,<0.1.0)",
-    "intervaltree (>=3.0,<4.0)",
+    "ruranges (>=0.1.4,<0.2.0)",
     "setuptools (>=80.2.0,<81.0.0)",
     "mappy (>=2.28,<3.0)",
     "numpy (>=2.2,<3.0)",
     "scipy (>=1.15.1,<2.0.0)",
     "matplotlib (>=3.10.0,<4.0.0)",
-    "intervaltree (>=3.2.1,<4.0.0)",
 ]
 
 [project.optional-dependencies]
@@ -46,7 +45,6 @@ repository = "https://github.com/BrooksLabUCSC/flair/"
 flair = "flair.flair_cli:main"
 junctions_from_sam = "flair.junctions_from_sam:main"
 mark_intron_retention = "flair.mark_intron_retention:main"
-mark_productivity = "flair.mark_productivity:main"
 normalize_counts_matrix = "flair.normalize_counts_matrix:main"
 plot_isoform_usage = "flair.plot_isoform_usage:main"
 predictProductivity = "flair.predictProductivity:main"

src/flair/annotate_aaseq_with_uniprot.py

@@ -5,7 +5,8 @@
 def parse_args():
     parser = argparse.ArgumentParser(description='''for annotating FLAIR aaseq predictions with UniProt names''')
     parser.add_argument('-i', '--input_aaseq', required=True, help='protein sequence fasta file - sequence names should be aaseqID_geneID')
-    parser.add_argument('-r', '--reference_seq', required=True, help='reference protein sequence fasta file - sequence names should be db|refID|geneName_organism, ex: >sp|Q8NH21|OR4F5_HUMAN. This script will extract the refID. this can be fasta or fasta.gz')
+    parser.add_argument('-r', '--reference_seq', required=True,
+                        help='reference protein sequence fasta file - sequence names should be db|refID|geneName_organism, ex: >sp|Q8NH21|OR4F5_HUMAN. This script will extract the refID. this can be fasta or fasta.gz')
     parser.add_argument('-o', '--output', required=True, help='output name - should be a fasta file')
     args = parser.parse_args()
     return args
@@ -44,10 +45,12 @@ def annotate_input(input_aaseq, ref_seq_to_name, output_name):
         out.write('>' + info + '\n' + seq + '\n')
     out.close()
 
+
 def main():
     args = parse_args()
     ref_seq_to_name = process_reference(args.reference_seq)
     annotate_input(args.input_aaseq, ref_seq_to_name, args.output)
 
+
 if __name__ == '__main__':
-    main()
+    main()

src/flair/annotate_group_vcf_vars.py

@@ -32,7 +32,7 @@ def main():
     out = open(args.output, 'w')
     for chrom, region in vcfvars:
         annotate_vars_in_region(vcfvars[(chrom, region)], chrom, region, out)
-        
+
 
 if __name__ == "__main__":
     main()

src/flair/annotation_data.py

@@ -0,0 +1,139 @@
+"""External gene annotation data loaded from GTF.
+
+Provides indexed lookups for junction-to-gene mapping, gene strand,
+transcript exon structure, and single-exon gene tracking.  Used by
+flair_transcriptome and flair_spliceevents for read correction, gene
+assignment, and isoform filtering.
+"""
+
+from flair.isoform_data import Exon, exons_to_juncs
+
+
+class AnnotData(object):
+    def __init__(self):
+        # map of (transcript_id, gene_id) -> tuple of Exon
+        self.transcript_to_exons = {}
+
+        # list of (transcript_id, gene_id, strand)
+        self.transcripts = []
+
+        # map of junction chain tuple -> (transcript_id, gene_id)
+        self.juncchain_to_transcript = {}
+
+        # map of Junc -> set of (transcript_id, gene_id)
+        self.junc_to_gene = {}
+
+        # single-exon annotations by strand: {'+': [], '-': []}
+        # each entry is Exon(start, end, gene_id), sorted for binary search
+        # FIXME: rename once it is figured out how this works in get_single_exon_gene_overlaps
+        self.all_annot_SE = {'+': [], '-': []}
+
+        # map of strand -> gene_id -> set of Exon
+        # FIXME: why is strand needed here
+        self.spliced_exons = {'+': {}, '-': {}}
+
+        # map of gene_id -> set of Junc
+        self.gene_to_annot_juncs = {}
+
+        # map of gene_id -> strand
+        self.gene_to_strand = {}
+
+        # map of gene_id -> tuple of sorted exon coordinate tuples
+        # union of all exons across all transcripts in the gene
+        self.gene_to_exons = {}
+
+        # map of junction chain tuple -> gene_id
+        self.sjc_to_gene = {}
+
+        # map of "transcript_id_gene_id" -> junction chain tuple
+        # key format matches BED name field from gtf_to_bed --include_gene
+        self.transcript_to_sjc = {}
+
+        # map of Junc -> gene_id (single value, last gene seen wins)
+        # used by spliceevents for simple junction-to-gene lookup
+        self.junc_to_gene_id = {}
+
+        self.gene_to_cds_starts = {}
+
+        self.transcript_to_nmd_except = {}
+
+        self.start_codon_count = 0
+
+
+def annot_data_from_gtf(gtf_data, region):
+    """Build AnnotData for a region from a pre-partitioned GtfData object."""
+    annots = AnnotData()
+    if gtf_data is None:
+        return annots
+    region_map = {region: annots}
+    for trans in gtf_data.transcripts:
+        if len(trans.exons) > 0:
+            _process_transcript(annots, region, region_map, trans)
+    # finalize gene_to_exons as sorted tuples
+    for gene_id in annots.gene_to_exons:
+        annots.gene_to_exons[gene_id] = tuple(sorted(annots.gene_to_exons[gene_id]))
+    return annots
+
+def _process_transcript(annots, region, region_map, trans):
+    trans.gene_id = trans.gene_id.replace('_', '-')
+    exons = [Exon(exon.start, exon.end) for exon in trans.exons]
+    sorted_exons = sorted(exons)
+    t_start = sorted_exons[0].start
+    t_end = sorted_exons[-1].end
+    _save_transcript_annot(trans.transcript_id, trans.gene_id, region,
+                           region_map, t_start, t_end, trans.strand, sorted_exons,
+                           trans.attrs['tag'], trans.start_codon)
+
+def _save_cds_starts(gene_id, start_codon, strand, annots):
+    if gene_id not in annots.gene_to_cds_starts:
+        annots.gene_to_cds_starts[gene_id] = set()
+    if start_codon is not None:
+        annots.start_codon_count += 1
+        if strand == '+':
+            annots.gene_to_cds_starts[gene_id].add(start_codon.start)
+        else:
+            annots.gene_to_cds_starts[gene_id].add(start_codon.end)
+
+def _save_spliced_transcript_info(gene_id, t_exons, juncs, transcript_id, strand, annots):
+    if gene_id not in annots.spliced_exons[strand]:
+        annots.spliced_exons[strand][gene_id] = set()
+    annots.spliced_exons[strand][gene_id].update(set(t_exons))
+    annots.juncchain_to_transcript[tuple(juncs)] = (transcript_id, gene_id)
+    annots.sjc_to_gene[tuple(juncs)] = gene_id
+    annots.transcript_to_sjc[f"{transcript_id}_{gene_id}"] = tuple(juncs)
+    if gene_id not in annots.gene_to_annot_juncs:
+        annots.gene_to_annot_juncs[gene_id] = set()
+    for j in juncs:
+        if j not in annots.junc_to_gene:
+            annots.junc_to_gene[j] = set()
+        annots.junc_to_gene[j].add((transcript_id, gene_id))
+        annots.junc_to_gene_id[j] = gene_id
+        annots.gene_to_annot_juncs[gene_id].add(j)
+
+def _save_transcript_annot(transcript_id, gene_id, region, region_map, t_start, t_end,
+                           strand, t_exons, transcript_tags, start_codon):
+    # region is a SeqRegion object
+    annots = region_map[region]
+
+    _save_cds_starts(gene_id, start_codon, strand, annots)
+    annots.transcript_to_nmd_except[transcript_id] = False
+    if 'NMD_exception' in transcript_tags:
+        annots.transcript_to_nmd_except[transcript_id] = True
+
+    annots.transcript_to_exons[(transcript_id, gene_id)] = tuple(t_exons)
+    juncs = exons_to_juncs(t_exons)
+    annots.transcripts.append((transcript_id, gene_id, strand))
+    if gene_id not in annots.gene_to_strand:
+        annots.gene_to_strand[gene_id] = strand
+    # accumulate exons per gene (as coordinate tuples for spliceevents compatibility)
+    exon_coords = set((e.start, e.end) for e in t_exons)
+    if gene_id not in annots.gene_to_exons:
+        annots.gene_to_exons[gene_id] = exon_coords
+    else:
+        annots.gene_to_exons[gene_id].update(exon_coords)
+    if len(juncs) == 0:
+        annots.all_annot_SE[strand].append(Exon(t_start, t_end, gene_id))
+    else:
+        _save_spliced_transcript_info(gene_id, t_exons, juncs, transcript_id, strand, annots)
+    for strand in ['+', '-']:
+        annots.all_annot_SE[strand] = sorted(annots.all_annot_SE[strand])  # FIXME: make set? Colette note: needs to be sorted for binary search later