test: cover annotation coordinate and export contracts

.github/workflows/unit-tests.yml

@@ -12,6 +12,8 @@ jobs:
       - uses: actions/setup-python@v5
         with:
           python-version: '3.10'
+      - name: Install BEDTools
+        run: sudo apt-get update && sudo apt-get install -y bedtools
       - name: Install test dependencies
         run: python -m pip install pytest gffutils pybedtools pandas numpy future tqdm
       - name: Run unit tests

annotator/annotate_bed.py

@@ -38,6 +38,12 @@
 
 
 def do_cprofile(func):
+    """Decorate a callable to profile its execution and print call statistics.
+
+    :param func: Callable accepting arbitrary positional and keyword arguments.
+    :return: Wrapper preserving the callable's result and exceptions, printing
+        profiling statistics in either case.
+    """
     def profiled_func(*args, **kwargs):
         profile = cProfile.Profile()
         try:
@@ -202,19 +208,14 @@ def hash_features(db, chromosomes, append_chr, fuzzy,
 
 def get_all_overlapping_features_from_query_stranded(
         chrom, qstart, qend, strand, features_dict):
-    """
-    Given a query location (chr, start, end), return all features that
-    overlap by at least one base. Functions similarly to gffutils db.region(),
-    but uses the pre-hashed self.chr19_features_dict to greatly speed things up.
-
-    :param chrom : string
-    :param qstart : int
-    :param qend : int
-    :param strand : string
-    :param features_dict : dict
-        dictionary of positions for all features
-    :return features : list
-        list of gffutils.Feature objects.
+    """Find same-strand features overlapping a BED coordinate query.
+
+    :param chrom: Query sequence ID.
+    :param qstart: Zero-based inclusive start.
+    :param qend: Exclusive end.
+    :param strand: Query strand, + or -.
+    :param features_dict: defaultdict(list) indexed by (chrom, int(position / HASH_VAL), strand), containing one-based inclusive GFF features.
+    :return: A list of overlapping features with overlap attributes assigned. Features are visited in bucket order; features in several buckets may occur more than once.
     """
     features = []
     append = features.append
@@ -256,20 +257,14 @@ def get_all_overlapping_features_from_query_stranded(
 
 def get_all_overlapping_features_from_query_unstranded(
         chrom, qstart, qend, strand, features_dict):
-    """
-    Returns overlapping features from a query.
-
-    :param chrom: string
-        query chromosome/source
-    :param qstart: int
-        query start
-    :param qend: int
-        query end
-    :param strand: string
-        query strand
-    :param features_dict: dict
-        dictionary of positions for all features
-    :return:
+    """Find overlaps with optional preferred-strand fallback.
+
+    :param chrom: Query sequence ID.
+    :param qstart: Zero-based inclusive start.
+    :param qend: Exclusive end.
+    :param strand: Preferred strand (+ or -), or None to search both.
+    :param features_dict: Bucket index as used by get_all_overlapping_features_from_query_stranded.
+    :return: Same-strand features if found, otherwise opposite-strand features. With strand None, positive-strand results precede negative-strand results.
     """
     if strand is None:  # If no strand given, return features from both strands
         positive_features = get_all_overlapping_features_from_query_stranded(
@@ -305,20 +300,15 @@ def get_all_overlapping_features_from_query_unstranded(
 
 def get_all_overlapping_features_from_query(
         chrom, qstart, qend, strand, features_dict, stranded=True):
-    """
-    Returns overlapping features from a query.
-    If stranded is True, then return only the features that overlap on
-    the same strand. Otherwise, features returned need not to be with respect
-    to strand (although if strand is given, it will return same-stranded
-    features first before attempting to find those on the opposite strand).
+    """Dispatch a BED overlap query to the selected strand policy.
 
-    :param chrom:
-    :param qstart:
-    :param qend:
-    :param strand:
-    :param features_dict:
-    :param stranded:
-    :return:
+    :param chrom: Query sequence ID.
+    :param qstart: Zero-based inclusive start.
+    :param qend: Exclusive end.
+    :param strand: Query strand, or None for an unstranded search.
+    :param features_dict: Bucket index of GFF features.
+    :param stranded: True to require the requested strand; False to use preferred-strand fallback.
+    :return: The overlapping feature list from the chosen strand policy.
     """
     if stranded:
         return get_all_overlapping_features_from_query_stranded(
@@ -332,11 +322,10 @@ def get_all_overlapping_features_from_query(
 
 
 def parse_annotation_string(features_string):
-    """
-    Splits a feature string into a list of feature strings
+    """Split pipe-delimited feature annotations without changing order.
 
-    :param features_string:
-    :return:
+    :param features_string: Pipe-delimited annotation string.
+    :return: A list of strings, preserving empty fields.
     """
     features = features_string.split('|')
     return features
@@ -359,18 +348,12 @@ def is_protein_coding(transcript_type):
 def classify_transcript_type(
         transcript_type, gene_type, distinct_nc_tx=None
 ):
-    """
-    deprecates is_protein_coding().
-
-    :param transcript_type: string
-    :param gene_type: string
-    :param distinct_nc_tx: list or None
-        We can pull out any distinct noncoding regions using this parameter.
-        For example, we can specify distinct_nc_tx = ['lincRNA', 'pseudogene'],
-        and this function will NOT bin these categories into the catchall
-        'non_coding' category.
-    :return transcript_type_class: string
+    """Classify a transcript for annotation priority.
 
+    :param transcript_type: Transcript biotype string.
+    :param gene_type: Gene biotype string, currently unused.
+    :param distinct_nc_tx: Optional iterable of noncoding biotypes to retain separately.
+    :return: protein_coding or miRNA for those exact biotypes, a matching distinct_nc_tx value, or non_coding.
     """
 
     # if the transcript type is protein coding, easy. return 'protein coding'.
@@ -561,16 +544,11 @@ def prioritize_transcript_then_gene(
 
 
 def split_string(priority_string, delim=','):
-    """
-    Parses/splits a string with an expected format and returns
-    three important fields: region, gene_id, gene_name. Returns
-    'intergenic' for all three fields if intergenic
-
-    :param priority_string: string
-        colon-delimited string containing region geneid and
-        genename in fields 4-6 (0-based open).
-    :param delim: change
-    :return region, gene_id, gene_name:
+    """Extract region, gene ID, and gene name from an annotation.
+
+    :param priority_string: Formatted feature string with fields 4 through 6, or intergenic.
+    :param delim: Field delimiter, comma by default.
+    :return: A tuple of region, gene ID, and gene name; all three equal intergenic for that sentinel.
     """
     if priority_string == 'intergenic':
         return 'intergenic', 'intergenic', 'intergenic'
@@ -711,12 +689,10 @@ def annotate(
 
 
 def split_bed_interval(bed_string):
-    """
-    Takes a bed file as a line (tab delimited) and returns 6 fields
-    correctly typed.
-    :param bed_string: string
-    :return chrom, start, end, name, score, strand: string
-        bed interval
+    """Parse a tab-delimited BED6 line.
+
+    :param bed_string: BED6 text with optional trailing whitespace.
+    :return: A tuple of chromosome, integer start/end, name, score, and strand; name, score, and strand remain strings.
     """
     bed = bed_string.rstrip().split('\t')
     return bed[0], int(bed[1]), int(bed[2]), bed[3], bed[4], bed[5]

annotator/annotation_functions.py

@@ -11,14 +11,10 @@
 
 
 def get_keys(gtf_format):
-    """
-    Describes the language for each kind of gtf file 
-    (usually different between sources).
-    Defaults to gencode standard for hg19/mm9/mm10/hg38
+    """Return the feature and attribute vocabulary for a GTF source.
 
-    :param gtf_format: string
-        either one of: ['gencode'],'ensembl','wormbase','genbank','refseq'
-    :return:
+    :param gtf_format: Source name: gencode, ensembl, ce10, ce11, wormbase, genbank, dm6_genbank, refseq, or ucsc. Other names use gencode.
+    :return: A new dictionary mapping semantic names to source keys; unsupported UTR kinds have value None.
     """
     print("gtf_format (species is deprecated) is {}".format(gtf_format))
     if gtf_format == 'ensembl':
@@ -149,8 +145,11 @@ def get_gene_to_transcript_dict(db, gene_id_key, transcript_id_key):
 
 
 def get_longest_transcripts(genes_dict, transcripts_dict):
-    """
-    Returns a dictionary of genes : longest_transcript
+    """Select the longest genomic transcript span for each gene.
+
+    :param genes_dict: Mapping of gene IDs to ordered transcript ID lists.
+    :param transcripts_dict: Mapping of transcript IDs to start and end coordinates.
+    :return: A defaultdict(dict) mapping each gene to its longest transcript ID. The first transcript wins ties; an empty list yields an empty string.
     """
     longest_genes_dict = defaultdict(dict)
     for gene, transcripts in iteritems(genes_dict):
@@ -166,9 +165,11 @@ def get_longest_transcripts(genes_dict, transcripts_dict):
 
 
 def most_upstream_downstream_positions(genes_dict, transcripts_dict):
-    """
-    Assumes the most downstream start < most downstream end (+, flipped for -) foreach transcript.
-    Returns a dictionary of genes : 
+    """Collect the inclusive outer transcript bounds for each gene.
+
+    :param genes_dict: Mapping of gene IDs to nonempty transcript ID lists.
+    :param transcripts_dict: Mapping of transcript IDs to start/end coordinates below one billion.
+    :return: A defaultdict(dict) mapping genes to minimum start and maximum end; independent of strand.
     """
     d = defaultdict(dict)
     for gene, transcripts in iteritems(genes_dict):
@@ -184,12 +185,25 @@ def most_upstream_downstream_positions(genes_dict, transcripts_dict):
     return d
 
 def get_premrna_lengths(gene_id, upstream_downstream):
+    """Return the inclusive genomic span of one gene.
+
+    :param gene_id: Gene ID to look up.
+    :param upstream_downstream: Mapping of gene IDs to inclusive start/end coordinates.
+    :return: The integer end minus start plus one.
+    """
     start = upstream_downstream[gene_id]['start']
     end = upstream_downstream[gene_id]['end']
     return end - start + 1    # (coords are 1-based inclusive still)
 
 
 def get_mrna_lengths(gene_id, exons_dict, genes_dict):
+    """Measure the union of a gene's exon bases across transcripts.
+
+    :param gene_id: Gene ID to look up.
+    :param exons_dict: Mapping of transcript IDs to exon dictionaries with chrom, inclusive start/end, and strand.
+    :param genes_dict: Mapping of gene IDs to transcript IDs.
+    :return: The integer length of the sorted exon union, merging without strand specificity. Requires the bedtools executable.
+    """
     exons_list = []
     total_mrna_length = 0
     for transcript in genes_dict[gene_id]:
@@ -210,12 +224,11 @@ def get_mrna_lengths(gene_id, exons_dict, genes_dict):
 
 
 def get_all_cds_dict(db, cds_key):
-    """
-    For every cds-annotated transcript id (ENST), return a
-    dictionary containing the lowest and highest
-    cds start and end vals for that transcript.
+    """Collect the inclusive outer CDS bounds for each transcript.
 
-    :return chr19_cds_dict : defaultdict{transcript:{'start':START, 'end':END}}
+    :param db: FeatureDB containing CDS features with transcript_id attributes.
+    :param cds_key: Feature type to treat as CDS.
+    :return: A defaultdict(dict) mapping transcript IDs to low and hi coordinates.
     """
     # chr19_cds_dict = defaultdict(lambda: {'low': MAXVAL, 'hi': MINVAL})
     cds_dict = defaultdict(dict)
@@ -274,10 +287,10 @@ def get_all_transcripts_dict(db, transcript_key, transcript_id_key):
 
 
 def chromosome_set(db):
-    """
-    Returns the set of chromosomes that exist in a database.
+    """Return the distinct sequence IDs present in a feature database.
 
-    :return:
+    :param db: FeatureDB exposing execute for the features table.
+    :return: A set of sequence ID strings, empty when no features exist.
     """
     ret = db.execute("SELECT seqid FROM features").fetchall()
     all_chromosomes = [r['seqid'] for r in ret]
@@ -317,16 +330,11 @@ def gene_id_to_name(db, gene_name_key):
 
 
 def classify_utr(utr_feature, cds_dict):
-    """
-    Given a feature classified as a UTR, return whether or not it is
-    upstream (5') or downstream (3') based on CDS positions
-
-    :param utr_feature: gffutils.Feature
-        feature already classified as UTR (within a coding transcript but
-        outside of CDS regions).
-    :param cds_dict: dict
-        dictionary containing cds positions for every transcript
-    :return:
+    """Classify a UTR relative to the coding span and transcript strand.
+
+    :param utr_feature: GFF feature with inclusive coordinates, strand, and transcript_id attributes.
+    :param cds_dict: Mapping of transcript IDs to inclusive low and hi CDS bounds.
+    :return: 5utr, 3utr, or unclassified_utr when CDS bounds or strand do not permit classification. A 5-prime classification takes precedence across transcripts.
     """
     three_prime_utr = False
     five_prime_utr = False

annotator/annotator.py

@@ -65,12 +65,11 @@
 
 
 def parse_annotation_priority(priority_file, delim=','):
-    """
-    From a comma delimited file, return a list of lists
-    describing the priority with which to annotate stuff.
+    """Read class/region priority pairs in file order.
 
-    :param priority_file:
-    :return:
+    :param priority_file: Text path with one two-field class/region row per line.
+    :param delim: Field delimiter, comma by default.
+    :return: A list of two-element lists preserving file order.
     """
     priority = []
     with open(priority_file, 'r') as f:
@@ -81,6 +80,12 @@ def parse_annotation_priority(priority_file, delim=','):
 
 def main():
     # Setup argument parser
+    """Annotate BED6 queries from command-line input/output and gtfdb paths.
+
+    Reads sys.argv. Optional format, priority files, chromosome limits, append-chr,
+    and unstranded flags select annotation policy. Returns None and writes
+    annotated rows through annotate_bed; normal CLI file errors propagate.
+    """
     parser = ArgumentParser()
 
     parser.add_argument(

annotator/build_gffutils_db.py

@@ -38,6 +38,12 @@ def build_db(
 
 def main():
     # Setup argument parser
+    """Build a gffutils database from command-line annotation_file and db_file paths.
+
+    Reads sys.argv; force permits overwrite and disable_infer_genes/transcripts
+    flags control inference. Returns None after writing the database.
+    Database/file errors propagate.
+    """
     parser = ArgumentParser()
 
     parser.add_argument(

annotator/create_AS_STRUCTURE.py

@@ -9,13 +9,13 @@
 from . import annotation_functions as af
 
 def create_as_structure(db, species, out_file):
-    """
-    
-    :param db: gffutils.FeatureDB
-    :param species: string
-        sets the gtf nomenclature
-    :param out_file: string
-    :return: 
+    """Write gene-level transcript spans and exon-union lengths as GFF rows.
+
+    :param db: FeatureDB containing gene, transcript, and exon features.
+    :param species: Source vocabulary name accepted by get_keys.
+    :param out_file: Destination path, overwritten with tab-delimited GFF rows.
+    :return: None; attributes include gene_id, mrna_length (exon union), and
+        premrna_length (inclusive outer transcript span). Requires bedtools.
     """
     keys = af.get_keys(species)
     genes = af.get_gene_to_transcript_dict(
@@ -69,9 +69,11 @@ def create_as_structure(db, species, out_file):
 
 
 def main():
-    """
-    Main program.
+    """Write gene span and exon-union lengths from command-line database options.
 
+    Reads sys.argv: db_file and out_file are required; species defaults to hg19.
+    Returns None and writes GFF rows. Requires bedtools. File/database errors
+    propagate; argparse reports invalid options.
     """
     parser = argparse.ArgumentParser()
     parser.add_argument(