From 90c5abc6b3c8e3f516f89ece8e0413431c0caf21 Mon Sep 17 00:00:00 2001 From: Brian Yee Date: Mon, 12 Oct 2026 22:19:41 -0700 Subject: [PATCH 1/2] test: cover annotation coordinate and export contracts --- .github/workflows/unit-tests.yml | 2 + annotator/annotate_bed.py | 118 ++++++--------- annotator/annotation_functions.py | 68 +++++---- annotator/annotator.py | 15 +- annotator/build_gffutils_db.py | 6 + annotator/create_AS_STRUCTURE.py | 20 +-- annotator/create_region_bedfiles.py | 170 +++++++++------------ annotator/gene_name2id.py | 76 ++++++---- annotator/get_region_lengths.py | 28 ++-- annotator/miRNA_name2id.py | 104 ++++++++----- unit_tests/test_cli_contracts.py | 94 ++++++++++++ unit_tests/test_coordinate_contracts.py | 187 ++++++++++++++++++++++++ unit_tests/test_identifier_contracts.py | 128 ++++++++++++++++ unit_tests/test_region_contracts.py | 173 ++++++++++++++++++++++ 14 files changed, 889 insertions(+), 300 deletions(-) create mode 100644 unit_tests/test_cli_contracts.py create mode 100644 unit_tests/test_coordinate_contracts.py create mode 100644 unit_tests/test_identifier_contracts.py create mode 100644 unit_tests/test_region_contracts.py diff --git a/.github/workflows/unit-tests.yml b/.github/workflows/unit-tests.yml index 48ad792..605cb54 100644 --- a/.github/workflows/unit-tests.yml +++ b/.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 diff --git a/annotator/annotate_bed.py b/annotator/annotate_bed.py index dc0a15b..b4fa839 100755 --- a/annotator/annotate_bed.py +++ b/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] diff --git a/annotator/annotation_functions.py b/annotator/annotation_functions.py index c2068a2..a5394a4 100644 --- a/annotator/annotation_functions.py +++ b/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 diff --git a/annotator/annotator.py b/annotator/annotator.py index a4bae86..49d9d1c 100755 --- a/annotator/annotator.py +++ b/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( diff --git a/annotator/build_gffutils_db.py b/annotator/build_gffutils_db.py index 51f00a6..a00399a 100755 --- a/annotator/build_gffutils_db.py +++ b/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( diff --git a/annotator/create_AS_STRUCTURE.py b/annotator/create_AS_STRUCTURE.py index e5bb8ab..ca06bef 100755 --- a/annotator/create_AS_STRUCTURE.py +++ b/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( diff --git a/annotator/create_region_bedfiles.py b/annotator/create_region_bedfiles.py index 942f1e7..030ff0a 100755 --- a/annotator/create_region_bedfiles.py +++ b/annotator/create_region_bedfiles.py @@ -29,18 +29,12 @@ from tqdm import trange def create_bedtools(features, keys, by_transcript=False): - """ - Given a list of features and chr19_keys dictionary, create a bedtool - containing intervals of features whose name is specified - using chr19_keys['gene_id'] - - :param features: list - list of gffutils features (1-based) for which to convert - to bedtool intervals - :param keys: dict - a set of chr19_keys and values which helps translate different - GTF/GFF nomenclatures (ie. 'cds' - :return: + """Convert inclusive GFF features to half-open BED intervals. + + :param features: Iterable of one-based inclusive GFF features. + :param keys: Source vocabulary from annotation_functions.get_keys. + :param by_transcript: True to name intervals by transcript and preserve individual features; False to name by gene and merge exported overlapping intervals. + :return: A BedTool with chromosome, converted start/end, selected ID, score 0, and strand. Multiple IDs produce separate intervals. """ intervals = [] key = 'transcript_id' if by_transcript else 'gene_id' @@ -59,25 +53,10 @@ def create_bedtools(features, keys, by_transcript=False): def merge_bedtool_by_gene(bedtool): - """ - Takes a bedtool and does a merge, but preserves the distinct names. - This is different than using bedtools merge. - - ie. - - chr1 100 200 GENE1 0 + - chr1 150 230 GENE1 0 + - chr1 150 250 GENE2 0 + + """Merge BED intervals separately by chromosome, gene name, and strand. - -> - - chr1 100 230 GENE1 0 + - chr1 150 250 GENE2 0 + - - :param bedtool: pybedtools.BedTool - - :return merged: pandas.DataFrame - table containing unique merged entries per distinct name. + :param bedtool: BED6 BedTool; intervals within each chromosome and gene are merged strand-specifically. + :return: A DataFrame with chrom/start/end/name/score/strand columns, sorted by chromosome and start. Requires the bedtools executable. """ df = bedtool.to_dataframe() merged = pd.DataFrame( @@ -103,17 +82,13 @@ def merge_bedtool_by_gene(bedtool): return merged def create_cds_region_bedfile(db, keys, cds_out, by_transcript=False): - """ - Creates the CDS region bedfile from a gffutils FeatureDB. - - :param db: gffutils.FeatureDB - gffutils sqlite db object - :param keys: dict - {'region':region defined in GTF}. Sets the naming scheme - for GTF parsing - :param cds_out: str - file name for output bedfile - :return cds: + """Return CDS intervals and optionally write their selected-name BED export. + + :param db: gffutils.FeatureDB containing the requested feature types. + :param keys: Source vocabulary from annotation_functions.get_keys. + :param cds_out: Optional CDS BED output path, or None. + :param by_transcript: True to name intervals by transcript and preserve individual features; False to name by gene and merge exported overlapping intervals. + :return: The original unmerged CDS BedTool; a written export is merged by gene unless by_transcript is True. """ cds_features = [] @@ -135,17 +110,13 @@ def create_cds_region_bedfile(db, keys, cds_out, by_transcript=False): return cds def create_exon_region_bedfile(db, keys, exon_out, by_transcript=False): - """ - Creates the exon region bedfile from a gffutils FeatureDB. - - :param db: gffutils.FeatureDB - gffutils sqlite db object - :param keys: dict - {'region':region defined in GTF}. Sets the naming scheme - for GTF parsing - :param cds_out: str - file name for output bedfile - :return: + """Write exon intervals as a BED export. + + :param db: gffutils.FeatureDB containing the requested feature types. + :param keys: Source vocabulary from annotation_functions.get_keys. + :param exon_out: Exon BED output path. + :param by_transcript: True to name intervals by transcript and preserve individual features; False to name by gene and merge exported overlapping intervals. + :return: None; writes a BED6 export, sorted when by_transcript is True. """ exon_features = [] @@ -166,17 +137,12 @@ def create_exon_region_bedfile(db, keys, exon_out, by_transcript=False): def create_gene_region_bedfile(db, keys, gene_out): - """ - Creates the exon region bedfile from a gffutils FeatureDB. - - :param db: gffutils.FeatureDB - gffutils sqlite db object - :param keys: dict - {'region':region defined in GTF}. Sets the naming scheme - for GTF parsing - :param cds_out: str - file name for output bedfile - :return: + """Write gene features as half-open BED intervals. + + :param db: gffutils.FeatureDB containing the requested feature types. + :param keys: Source vocabulary from annotation_functions.get_keys. + :param gene_out: Gene BED output path. + :return: None; writes the unmerged gene intervals. """ gene_features = [] gene_hash = defaultdict(list) @@ -202,26 +168,15 @@ def create_gene_region_bedfile(db, keys, gene_out): # ) def create_utr_region_bedfiles(db, keys, cds_dict, utr3_out, utr5_out, by_transcript=False): - """ - Creates UTR bedfiles from a gffutils FeatureDB. - - :param db: gffutils.FeatureDB - gffutils sqlite db object - :param keys: dict - {'region':region defined in GTF}. Sets the naming scheme - for GTF parsing - :param cds_dict: dict - dictionary of {transcript:cds_regions} - :param utr3_out: str - file name for output bedfile - :param utr5_out: str - file name for output bedfile - :param by_transcript: bool - if True, return all UTR regions, regardless (every UTR in every transcript is reported) - if False, return all UTR regions as a merged region according to gene - (merged, unoverlapping UTR region is reported for every gene) - - :return: + """Return and optionally export five-prime and three-prime UTR intervals. + + :param db: gffutils.FeatureDB containing the requested feature types. + :param keys: Source vocabulary from annotation_functions.get_keys. + :param cds_dict: Mapping of transcript IDs to inclusive low/hi CDS bounds, used to classify generic UTR features. + :param utr3_out: Optional three-prime UTR output path, or None. + :param utr5_out: Optional five-prime UTR output path, or None. + :param by_transcript: True to name intervals by transcript and preserve individual features; False to name by gene and merge exported overlapping intervals. + :return: A tuple (five_prime, three_prime) of unmerged BedTools. Supports generic UTR or explicit utr5/utr3 feature types. Returns 1 with a diagnostic if neither vocabulary is available. """ five_prime_utr_features = [] three_prime_utr_features = [] @@ -272,23 +227,17 @@ def create_utr_region_bedfiles(db, keys, cds_dict, utr3_out, utr5_out, by_transc def create_intron_region_bedfiles(db, exons_dict, transcripts_dict, keys, proxintron_out, distintron_out, allintron_out, by_transcript): - """ - Creates prox and distal intron bedfiles from a gffutils FeatureDB. - - :param db: gffutils.FeatureDB - gffutils sqlite db object - :param keys: dict - {'region':region defined in GTF}. Sets the naming scheme - for GTF parsing - :param exons_dict: dict - dictionary of {transcript:exon_regions} - :param transcripts_dict: dict - dictionary of {transcript:gffutils.Features} - :param proxintron_out: str - file name for output bedfile - :param distintron_out: str - file name for output bedfile - :return: + """Return and optionally export proximal, distal, and complete intron intervals. + + :param db: gffutils.FeatureDB containing the requested feature types. + :param exons_dict: Mapping of transcript IDs to inclusive exon coordinate dictionaries. + :param transcripts_dict: Mapping of transcript IDs to inclusive start/end dictionaries. + :param keys: Source vocabulary from annotation_functions.get_keys. + :param proxintron_out: Optional proximal intron output path, or None. + :param distintron_out: Optional distal intron output path, or None. + :param allintron_out: Optional complete intron output path, or None. + :param by_transcript: True to name intervals by transcript and preserve individual features; False to name by gene and merge exported overlapping intervals. + :return: A tuple (proximal, distal, all_introns) of unmerged BedTools. Proximal regions lie within 500 bases of an exon; introns of length at most 1000 are entirely proximal. """ allintrons = [] proxintrons = [] @@ -384,6 +333,21 @@ def create_region_bedfiles( ): ### get the parsed regions ### + """Write selected BED exports from a feature database. + + :param db_file: Path to a gffutils database. + :param species: Source vocabulary accepted by get_keys. + :param cds_out: Optional CDS BED destination. + :param utr3_out: Optional three-prime UTR BED destination. + :param utr5_out: Optional five-prime UTR BED destination. + :param proxintron_out: Optional proximal intron BED destination. + :param distintron_out: Optional distal intron BED destination. + :param allintron_out: Optional complete intron BED destination. + :param gene_out: Optional gene BED destination. + :param exon_out: Optional exon BED destination. + :param by_transcript: Preserve transcript names instead of gene merging. + :return: None; only non-None destinations select exports. Requires bedtools. + """ keys = af.get_keys(species) db = gffutils.FeatureDB(db_file) featuretypes = [x.lower() for x in list(db.featuretypes())] @@ -416,6 +380,12 @@ def create_region_bedfiles( def main(): # Setup argument parser + """Export requested BED regions from command-line db_file and species options. + + Reads sys.argv. Optional cds/utr3/utr5/proxintron/distintron/allintron/gene/exon + output paths select exports; by-transcript retains transcript names. + Returns None after writing requested files. Requires bedtools. + """ parser = ArgumentParser() parser.add_argument( diff --git a/annotator/gene_name2id.py b/annotator/gene_name2id.py index 7d3fc23..22e4637 100755 --- a/annotator/gene_name2id.py +++ b/annotator/gene_name2id.py @@ -25,7 +25,13 @@ def get_id_from_namerow(row, name2id, name_col): - """ Uses name2id dict to return the id with the associated name. """ + """Look up ordered gene IDs for one table row. The lookup value is converted to str. + + :param row: Mapping or pandas row holding the lookup column. + :param name2id: Mapping of names to ordered ID lists. + :param name_col: Column name used for lookup. + :return: A pipe-joined ID string, or an empty string for an unknown key. + """ try: return '|'.join(name2id[str(row[name_col])]) except KeyError: @@ -33,26 +39,24 @@ def get_id_from_namerow(row, name2id, name_col): return '' def add_id_to_dataframe(df, name2id_dict, name_col): - """ adds id values to pandas dataframe """ + """Append the 'gene id' lookup column to a DataFrame. Raises AssertionError if the ID column already exists. + + :param df: Pandas DataFrame to update in place. + :param name2id_dict: Mapping of names to ordered ID lists. + :param name_col: Column name used for lookup. + :return: The same DataFrame with joined IDs; unknown keys have empty strings. + """ assert 'gene id' not in df.columns df['gene id'] = df.apply(get_id_from_namerow, args=[name2id_dict, name_col], axis=1) return df def build_name_to_id_dict(gffdb_file, custom_file): - """ - Returns {gene_name:gene_id} dictionary. - Names with more than 1 associated ID will be appended, delimited by |. - - :param gffdb_file: string - path of the gtfdb. - Database must contain: - 'gene' in featuretype field - 'gene_name' in feature.attributes.chr19_keys() - 'gene_id' in feature.attributes.chr19_keys() - - :return name2id_dict: dict - {gene_name:gene_id[]} + """Build an identifier mapping from optional database and custom table sources. + + :param gffdb_file: Optional path to a gffutils database, or None. + :param custom_file: Optional tab-delimited mapping path; comment lines start with #. Gene mappings use ID/name columns; microRNA mappings use precursor ID/name/pipe-delimited mature IDs. + :return: A defaultdict(list) mapping names to ID lists; empty when both sources are None. Source errors are printed and do not propagate. """ # delim = '|' name2id = defaultdict(list) @@ -81,22 +85,15 @@ def build_name_to_id_dict(gffdb_file, custom_file): def read_and_append_id(in_file, sep, name_col, gffdb_file, custom_file, delete_original_column): - """ - Reads in a delimited file with miRNA names in at least one column - Returns dataframe containing matching miRNA precursor ids. - - :param in_file: string - input delimited file - :param sep: string - delimiter - :param name_col: string - column name where miRNA names exist - :param gffdb_file: string - path to gffutils database file - :param delete_original_column: bool - if True, delete the original name column (specified with name_col) - if False, keep it. - :return: + """Read a header-bearing table and append gene IDs. + + :param in_file: Input table path with a header row. + :param sep: Input/output table delimiter. + :param name_col: Column name used for lookup. + :param gffdb_file: Optional path to a gffutils database, or None. + :param custom_file: Optional tab-delimited mapping path; comment lines start with #. Gene mappings use ID/name columns; microRNA mappings use precursor ID/name/pipe-delimited mature IDs. + :param delete_original_column: True to drop the name column and index by ID; False to index by the original name. + :return: A DataFrame indexed by original names or IDs according to delete_original_column. """ df = pd.read_table(in_file, sep=sep) @@ -112,6 +109,17 @@ def read_and_append_id(in_file, sep, name_col, gffdb_file, custom_file, delete_o def convert(in_file, sep, name_col, out_file, gtfdb_file, custom_file, delete_original_column): + """Convert a header-bearing identifier table and write its indexed result. + + :param in_file: Input table path with a header row. + :param sep: Input/output table delimiter. + :param name_col: Column name used for lookup. + :param out_file: Output table path; includes header and index. + :param gtfdb_file: Optional path to a gffutils database, or None. + :param custom_file: Optional tab-delimited mapping path; comment lines start with #. Gene mappings use ID/name columns; microRNA mappings use precursor ID/name/pipe-delimited mature IDs. + :param delete_original_column: True to drop the name column and index by ID; False to index by the original name. + :return: None; the output file contains the converted table. + """ df = read_and_append_id( in_file, sep, name_col, gtfdb_file, custom_file, delete_original_column @@ -121,6 +129,12 @@ def convert(in_file, sep, name_col, out_file, gtfdb_file, custom_file, delete_or def main(): # Setup argument parser + """Convert gene names using command-line input/output paths and mapping options. + + Reads sys.argv: custom and/or gffdb supply ID mappings; name_col identifies + the name column; sep defaults to tab. delete_original_column indexes by ID + instead of name. Returns None after writing the converted table. + """ parser = ArgumentParser() parser.add_argument( diff --git a/annotator/get_region_lengths.py b/annotator/get_region_lengths.py index 7b2d02a..6f3999e 100755 --- a/annotator/get_region_lengths.py +++ b/annotator/get_region_lengths.py @@ -178,14 +178,13 @@ def calculate_avg_utr_lengths(db, cds_dict, keys): def calculate_total_intron_lengths(db, exons_dict, transcripts_dict, keys): - """ - Calculates the total length of distal, proximal, and all introns. - - :param db: - :param exons_dict: - :param transcripts_dict: - :param keys: - :return: + """Print union coverage of proximal, distal, and complete intron regions. + + :param db: FeatureDB with transcript features. + :param exons_dict: Transcript IDs mapped to inclusive exon coordinates. + :param transcripts_dict: Transcript IDs mapped to inclusive start/end bounds. + :param keys: Source vocabulary from get_keys. + :return: None; prints each region's base coverage. Requires bedtools. """ proxintron_out = None # don't write to any file @@ -227,14 +226,11 @@ def calculate_avg_lengths(db_file, species): def calculate_total_lengths(db_file, species): - """" - Prints number of bases labeled as CDS, UTR, and intron in a gtf db file. - Can be used to generate ratios for region distribution pie charts. - - :param db_file: gffutils.FeatureDB - :param species: string - used to generate chr19_keys corresonding to GTF column nomenclature. - :return: + """Print CDS, UTR, and intron union coverage from a feature database. + + :param db_file: Path to a gffutils database containing transcripts and exons. + :param species: Source vocabulary name accepted by get_keys. + :return: None; prints measured base counts for each region. Requires bedtools. """ keys = af.get_keys(species) db = gffutils.FeatureDB(db_file) diff --git a/annotator/miRNA_name2id.py b/annotator/miRNA_name2id.py index eb285df..3901c41 100755 --- a/annotator/miRNA_name2id.py +++ b/annotator/miRNA_name2id.py @@ -27,7 +27,13 @@ def get_precursor_id_from_namerow(row, name2id, name_col): - """ Uses name2id dict to return the id with the associated name. """ + """Look up ordered microRNA precursor IDs for one table row. + + :param row: Mapping or pandas row holding the lookup column. + :param name2id: Mapping of names to ordered ID lists. + :param name_col: Column name used for lookup. + :return: A pipe-joined ID string, or an empty string for an unknown key. + """ try: return '|'.join(name2id[row[name_col]]) except KeyError: @@ -36,7 +42,13 @@ def get_precursor_id_from_namerow(row, name2id, name_col): def get_mature_id_from_idrow(row, precursor2mature, name_col): - """ Uses name2id dict to return the id with the associated name. """ + """Look up ordered mature microRNA IDs for one precursor row. + + :param row: Mapping or pandas row holding the lookup column. + :param precursor2mature: Mapping of precursor IDs to ordered mature ID lists. + :param name_col: Column name used for lookup. + :return: A pipe-joined ID string, or an empty string for an unknown key. + """ try: return '|'.join(precursor2mature[row[name_col]]) except KeyError: @@ -45,24 +57,37 @@ def get_mature_id_from_idrow(row, precursor2mature, name_col): def add_precursor_id_to_dataframe(df, name2id_dict, name_col): - """ adds id values to pandas dataframe """ + """Append the 'miRNA id' lookup column to a DataFrame. Raises AssertionError if the ID column already exists. + + :param df: Pandas DataFrame to update in place. + :param name2id_dict: Mapping of names to ordered ID lists. + :param name_col: Column name used for lookup. + :return: The same DataFrame with joined IDs; unknown keys have empty strings. + """ assert 'miRNA id' not in df.columns df['miRNA id'] = df.apply(get_precursor_id_from_namerow, args=[name2id_dict, name_col], axis=1) return df def add_mature_id_to_dataframe(df, precursor2mature_dict, name_col): + """Append the 'mature' lookup column to a DataFrame. + + :param df: Pandas DataFrame to update in place. + :param precursor2mature_dict: Mapping of precursor IDs to ordered mature ID lists. + :param name_col: Column name used for lookup. + :return: The same DataFrame with joined IDs; unknown keys have empty strings. + """ print(df.head()) df['mature'] = df.apply(get_mature_id_from_idrow, args=[precursor2mature_dict, name_col], axis=1) return df def build_precursor_id_to_mature_id_dict(gffdb_file, custom_file): - """ - Uses a gffdb file and a custom file to - :param gffdb_file: - :param custom_file: - :return: + """Build an identifier mapping from optional database and custom table sources. + + :param gffdb_file: Optional path to a gffutils database, or None. + :param custom_file: Optional tab-delimited mapping path; comment lines start with #. Gene mappings use ID/name columns; microRNA mappings use precursor ID/name/pipe-delimited mature IDs. + :return: A mapping of precursor IDs to mature ID lists. Supply one source when a precursor occurs in both; combining repeated precursor keys is currently unsupported. """ precursor2mature = defaultdict(list) if gffdb_file is not None: @@ -94,19 +119,11 @@ def build_precursor_id_to_mature_id_dict(gffdb_file, custom_file): def build_name_to_precursor_id_dict(gffdb_file, custom_file): - """ - Returns {name:id} dictionary. - Names with more than 1 associated ID will be appended, delimited by |. - - :param gffdb_file: string - path of the gffdb. - Database must contain: - 'miRNA_primary_transcript' in featuretype field - 'Name' in feature.attributes.chr19_keys() - 'ID' in feature.attributes.chr19_keys() - - :return name2id_dict: dict - {name:id[]} + """Build an identifier mapping from optional database and custom table sources. + + :param gffdb_file: Optional path to a gffutils database, or None. + :param custom_file: Optional tab-delimited mapping path; comment lines start with #. Gene mappings use ID/name columns; microRNA mappings use precursor ID/name/pipe-delimited mature IDs. + :return: A defaultdict(list) mapping names to ID lists; empty when both sources are None. Source errors are printed and do not propagate. """ # delim = '|' precursor_name2id = defaultdict(list) @@ -134,22 +151,16 @@ def build_name_to_precursor_id_dict(gffdb_file, custom_file): def read_and_append_id(in_file, sep, name_col, gffdb_file, custom_file, add_mature, delete_original_column): - """ - Reads in a delimited file with miRNA names in at least one column - Returns dataframe containing matching miRNA precursor ids. - - :param in_file: string - input delimited file - :param sep: string - delimiter - :param name_col: string - column name where miRNA names exist - :param gffdb_file: string - path to gffutils database file - :param delete_original_column: bool - if True, delete the original name column (specified with name_col) - if False, keep it. - :return: + """Read a header-bearing table and append microRNA precursor IDs. + + :param in_file: Input table path with a header row. + :param sep: Input/output table delimiter. + :param name_col: Column name used for lookup. + :param gffdb_file: Optional path to a gffutils database, or None. + :param custom_file: Optional tab-delimited mapping path; comment lines start with #. Gene mappings use ID/name columns; microRNA mappings use precursor ID/name/pipe-delimited mature IDs. + :param add_mature: True to append mature IDs using the precursor mapping. + :param delete_original_column: True to drop the name column and index by ID; False to index by the original name. + :return: A DataFrame indexed by original names or IDs according to delete_original_column. """ df = pd.read_table(in_file, sep=sep) @@ -169,7 +180,18 @@ def read_and_append_id(in_file, sep, name_col, gffdb_file, custom_file, add_matu def convert(in_file, sep, name_col, out_file, gtfdb_file, custom_file, add_mature, delete_original_column=True): - """ Main runner script for conversion. """ + """Convert a header-bearing identifier table and write its indexed result. + + :param in_file: Input table path with a header row. + :param sep: Input/output table delimiter. + :param name_col: Column name used for lookup. + :param out_file: Output table path; includes header and index. + :param gtfdb_file: Optional path to a gffutils database, or None. + :param custom_file: Optional tab-delimited mapping path; comment lines start with #. Gene mappings use ID/name columns; microRNA mappings use precursor ID/name/pipe-delimited mature IDs. + :param add_mature: True to append mature IDs using the precursor mapping. + :param delete_original_column: True to drop the name column and index by ID; False to index by the original name. + :return: None; the output file contains the converted table. + """ df = read_and_append_id( in_file, sep, name_col, gtfdb_file, custom_file, add_mature, @@ -181,6 +203,12 @@ def convert(in_file, sep, name_col, out_file, gtfdb_file, custom_file, add_matur def main(): # Setup argument parser + """Convert microRNA names using command-line input/output and mapping options. + + Reads sys.argv: custom and/or gffdb supply precursor IDs, name_col selects + the lookup column, and sep defaults to tab. add_mature appends mature IDs. + Returns None after writing the table indexed by precursor ID. + """ parser = ArgumentParser() parser.add_argument( diff --git a/unit_tests/test_cli_contracts.py b/unit_tests/test_cli_contracts.py new file mode 100644 index 0000000..ab1e02f --- /dev/null +++ b/unit_tests/test_cli_contracts.py @@ -0,0 +1,94 @@ +"""Exercise documented CLI entry points with small local input files.""" + +import gffutils +import pandas as pd +import pybedtools +import pytest + +from annotator import annotator, build_gffutils_db, create_AS_STRUCTURE +from annotator import create_region_bedfiles, gene_name2id, miRNA_name2id + + +@pytest.fixture +def database(tmp_path): + gtf = tmp_path / 'input.gtf' + gtf.write_text(''.join('chr1\tx\t{}\t{}\t{}\t.\t+\t.\tgene_id "g1"; gene_name "Alpha"; transcript_id "tx1"; transcript_type "protein_coding"; gene_type "protein_coding";\n'.format(*r) for r in [ + ('gene',1,100), ('transcript',1,100), ('exon',1,100), ('CDS',20,80), ('UTR',1,19), ('UTR',81,100)])) + path = tmp_path / 'input.db' + build_gffutils_db.build_db(str(gtf), str(path), True, True, True) + return gtf, path + + +def test_database_cli_builds_queryable_feature_database(database, tmp_path, monkeypatch): + gtf, _ = database + out = tmp_path / 'cli.db' + monkeypatch.setattr('sys.argv', ['build_gffutils_db', '--annotation_file',str(gtf),'--db_file',str(out), + '--disable_infer_genes','--disable_infer_transcripts']) + build_gffutils_db.main() + db = gffutils.FeatureDB(str(out)) + assert db.count_features_of_type('gene') == 1 + assert db.count_features_of_type('CDS') == 1 + assert [(f.start,f.end) for f in db.features_of_type('CDS')] == [(20,80)] + + +@pytest.mark.parametrize('by_transcript,name', [(False,'g1'),(True,'tx1')]) +def test_region_cli_writes_only_requested_exports(database, tmp_path, monkeypatch, by_transcript, name): + _, db = database + cds = tmp_path / 'cds.bed'; exons = tmp_path / 'exons.bed' + args = ['create_region_bedfiles','--db_file',str(db),'--species','gencode', + '--cds_out',str(cds),'--exons_out',str(exons)] + if by_transcript: args.append('--by-transcript') + monkeypatch.setattr('sys.argv',args) + create_region_bedfiles.main() + assert [(x.start,x.end,x.name) for x in pybedtools.BedTool(str(cds))] == [(19,80,name)] + assert [(x.start,x.end,x.name) for x in pybedtools.BedTool(str(exons))] == [(0,100,name)] + + +def test_annotation_cli_creates_coding_and_intergenic_rows(database, tmp_path, monkeypatch): + _,db = database + source = tmp_path / 'queries.bed'; out = tmp_path / 'annotated.bed' + source.write_text('chr1\t30\t40\tcoding\t0\t+\nchr1\t200\t210\toutside\t0\t+\n') + monkeypatch.setattr('sys.argv',['annotator','--input',str(source),'--output',str(out),'--gtfdb',str(db)]) + annotator.main() + rows = [line.split('\t') for line in out.read_text().splitlines()] + assert rows[0][:6] == ['chr1','30','40','coding','0','+'] + assert rows[0][6:9] == ['g1','Alpha','CDS'] + assert rows[1][6:9] == ['intergenic'] * 3 + + +def test_gene_cli_converts_custom_mapping_and_preserves_measurement(tmp_path, monkeypatch): + source = tmp_path / 'input.tsv'; source.write_text('symbol\tcount\nAlpha\t7\n') + mapping = tmp_path / 'map.tsv'; mapping.write_text('g1\tAlpha\n') + out = tmp_path / 'output.tsv' + monkeypatch.setattr('sys.argv',['gene_name2id','--input',str(source),'--output',str(out), + '--custom',str(mapping),'--name_col','symbol','--delete_original_column']) + gene_name2id.main() + assert pd.read_table(out).to_dict('list') == {'gene id':['g1'], 'count':[7]} + + +def test_mirna_cli_converts_custom_precursor_and_mature_ids(tmp_path, monkeypatch): + source = tmp_path / 'input.tsv'; source.write_text('symbol\tcount\nmir-one\t7\n') + mapping = tmp_path / 'map.tsv'; mapping.write_text('pre1\tmir-one\tmat1|mat2\n') + out = tmp_path / 'output.tsv' + monkeypatch.setattr('sys.argv',['miRNA_name2id','--input',str(source),'--output',str(out), + '--custom',str(mapping),'--name_col','symbol','--add_mature']) + miRNA_name2id.main() + assert pd.read_table(out).to_dict('list') == {'miRNA id':['pre1'], 'count':[7], 'mature':['mat1|mat2']} + + +def test_as_structure_cli_writes_measured_gene_attributes(database, tmp_path, monkeypatch): + _,db = database + out = tmp_path / 'structure.gff' + monkeypatch.setattr('sys.argv',['create_as_structure','--db_file',str(db),'--out_file',str(out),'--species','gencode']) + create_AS_STRUCTURE.main() + assert out.read_text().split('\t')[-1].strip() == 'gene_id=g1;mrna_length=100;premrna_length=100' + + +@pytest.mark.parametrize('module', [annotator, build_gffutils_db, create_region_bedfiles, + create_AS_STRUCTURE, gene_name2id, miRNA_name2id]) +def test_cli_help_lists_options_and_exits_successfully(module, monkeypatch, capsys): + monkeypatch.setattr('sys.argv', ['tool', '--help']) + with pytest.raises(SystemExit) as exc: + module.main() + assert exc.value.code == 0 + assert 'usage:' in capsys.readouterr().out diff --git a/unit_tests/test_coordinate_contracts.py b/unit_tests/test_coordinate_contracts.py new file mode 100644 index 0000000..02b9fd8 --- /dev/null +++ b/unit_tests/test_coordinate_contracts.py @@ -0,0 +1,187 @@ +"""Small synthetic examples for annotation coordinates and format contracts.""" + +from collections import defaultdict +from types import SimpleNamespace + +import gffutils +import pybedtools +import pytest + +from annotator import annotation_functions as af +from annotator import annotate_bed as ab + + +def feature(start, end, strand='+', transcripts=('tx',), featuretype='CDS'): + return gffutils.Feature(seqid='chr1', start=start, end=end, strand=strand, + featuretype=featuretype, + attributes={'transcript_id': list(transcripts)}) + + +@pytest.mark.parametrize('source,expected', [ + ('ensembl', ('CDS', None, 'three_prime_utr', 'five_prime_utr', 'transcript_biotype', 'gene_biotype')), + ('ce10', ('CDS', 'UTR', None, None, 'transcript_biotype', 'gene_biotype')), + ('wormbase', ('CDS', None, 'three_prime_utr', 'five_prime_utr', 'transcript_biotype', 'gene_biotype')), + ('ce11', ('CDS', None, 'three_prime_utr', 'five_prime_utr', 'transcript_biotype', 'gene_biotype')), + ('gencode', ('CDS', 'UTR', None, None, 'transcript_type', 'gene_type')), + ('unknown', ('CDS', 'UTR', None, None, 'transcript_type', 'gene_type')), + ('genbank', ('Genbank CDS', None, None, None, 'gene_biotype', 'gene_biotype')), + ('dm6_genbank', ('Genbank CDS', None, None, None, 'gene_biotype', 'gene_biotype')), + ('refseq', ('CDS', None, None, None, 'gene_biotype', 'gene_biotype')), + ('ucsc', ('CDS', None, '3UTR', '5UTR', 'gene_biotype', 'gene_type')), +]) +def test_get_keys_selects_source_vocabulary(source, expected): + keys = af.get_keys(source) + assert tuple(keys[k] for k in ('cds', 'utr', 'utr3', 'utr5', 'transcript_type', 'gene_type')) == expected + assert set(keys) == {'cds', 'utr', 'utr3', 'utr5', 'transcript_type', 'gene_type', + 'gene', 'gene_name', 'transcript', 'transcript_id', 'exon', 'gene_id'} + + +def test_get_keys_returns_independent_mappings(): + first = af.get_keys('gencode') + first['cds'] = 'changed' + assert af.get_keys('gencode')['cds'] == 'CDS' + + +def test_get_all_cds_dict_unions_extrema_for_each_transcript(): + db = SimpleNamespace(features_of_type=lambda kind: iter([ + feature(100, 110, transcripts=('tx', 'shared')), + feature(50, 75), feature(105, 150), feature(80, 90), + ]) if kind == 'coding' else iter([])) + assert af.get_all_cds_dict(db, 'coding') == { + 'tx': {'low': 50, 'hi': 150}, 'shared': {'low': 100, 'hi': 110}} + + +def test_get_all_cds_dict_empty_database(): + db = SimpleNamespace(features_of_type=lambda kind: iter([])) + assert af.get_all_cds_dict(db, 'CDS') == {} + + +def test_chromosome_set_deduplicates_database_rows(): + db = gffutils.create_db('chr1\tx\tgene\t1\t3\t.\t+\t.\tgene_id "g1";\n' + 'chr2\tx\tgene\t1\t3\t.\t-\t.\tgene_id "g2";\n' + 'chr1\tx\tgene\t8\t9\t.\t+\t.\tgene_id "g3";\n', + dbfn=':memory:', from_string=True) + assert af.chromosome_set(db) == {'chr1', 'chr2'} + + +@pytest.mark.parametrize('strand,start,end,expected', [ + ('+', 10, 49, '5utr'), ('+', 101, 120, '3utr'), + ('-', 10, 49, '3utr'), ('-', 101, 120, '5utr'), + ('+', 60, 80, 'unclassified_utr'), ('.', 10, 49, 'unclassified_utr'), +]) +def test_classify_utr_respects_coding_span_and_strand(strand, start, end, expected): + assert af.classify_utr(feature(start, end, strand), {'tx': {'low': 50, 'hi': 100}}) == expected + + +def test_classify_utr_without_coding_transcript_is_unclassified(): + assert af.classify_utr(feature(1, 10), {}) == 'unclassified_utr' + + +@pytest.mark.parametrize('length,distance,prox,dist', [ + (100, 50, [(10, 110)], []), + (101, 50, [(10, 60), (61, 111)], [(60, 61)]), + (200, 50, [(10, 60), (160, 210)], [(60, 160)]), +]) +def test_proximal_distal_split_preserves_half_open_coordinates(length, distance, prox, dist): + intron = pybedtools.create_interval_from_list(['chr2', '10', str(10 + length), 'tx', '7', '-']) + result = af.get_proxdist_from_intron(intron, distance) + assert [(x.start, x.end) for x in result['prox']] == prox + assert [(x.start, x.end) for x in result['dist']] == dist + assert sum(len(x) for values in result.values() for x in values) == length + assert all((x.chrom, x.strand, x.score) == ('chr2', '-', '7') for values in result.values() for x in values) + + +@pytest.mark.parametrize('exons,expected', [ + ([{'start': 1, 'end': 20}, {'start': 31, 'end': 50}], [{'start': 21, 'end': 30}]), + ([{'start': 31, 'end': 50}, {'start': 1, 'end': 20}], [{'start': 21, 'end': 30}]), + ([{'start': 1, 'end': 50}], []), + ([], []), +]) +def test_find_introns_uses_inclusive_exon_coordinates(exons, expected): + assert af.find_introns({'start': 1, 'end': 50}, exons) == expected + + +@pytest.mark.parametrize('tx_type,gene_type,distinct,expected', [ + ('protein_coding', 'lncRNA', None, 'protein_coding'), + ('miRNA', 'protein_coding', None, 'miRNA'), + ('retained_intron', 'protein_coding', None, 'non_coding'), + ('lncRNA', 'lncRNA', ['lncRNA'], 'lncRNA'), + ('lncRNA', 'lncRNA', ['snoRNA'], 'non_coding'), +]) +def test_classify_transcript_type_preserves_selected_noncoding_classes(tx_type, gene_type, distinct, expected): + assert ab.classify_transcript_type(tx_type, gene_type, distinct) == expected + + +@pytest.mark.parametrize('value,expected', [('protein_coding', True), ('miRNA', False), ('', False)]) +def test_is_protein_coding_requires_exact_type(value, expected): + assert ab.is_protein_coding(value) is expected + + +def test_split_bed_interval_converts_only_coordinates(): + assert ab.split_bed_interval('chr2\t0\t10\tpeak\t3.5\t-\n') == ('chr2', 0, 10, 'peak', '3.5', '-') + + +@pytest.mark.parametrize('line,delimiter,expected', [ + ('intergenic', ',', ('intergenic',) * 3), + ('tx,10,20,+,CDS,g1,Alpha,protein_coding,protein_coding,overlap', ',', ('CDS', 'g1', 'Alpha')), + ('tx:10:20:+:exon:g2:Beta:type:type:overlap', ':', ('exon', 'g2', 'Beta')), +]) +def test_split_string_extracts_public_annotation_columns(line, delimiter, expected): + assert ab.split_string(line, delimiter) == expected + + +@pytest.mark.parametrize('text,expected', [('a|b', ['a', 'b']), ('a|', ['a', '']), ('', [''])]) +def test_parse_annotation_string_preserves_feature_order(text, expected): + assert ab.parse_annotation_string(text) == expected + + +@pytest.mark.parametrize('strand,query,expected', [ + ('+', (9, 20), ['same']), ('-', (9, 20), ['other']), + ('+', (0, 10), ['same']), ('+', (20, 30), []), +]) +def test_stranded_overlap_observes_bed_to_gtf_boundary(strand, query, expected): + features = defaultdict(list) + a = feature(10, 20, '+'); a.id = 'same' + b = feature(10, 20, '-'); b.id = 'other' + features['chr1', 0, '+'] = [a]; features['chr1', 0, '-'] = [b] + result = ab.get_all_overlapping_features_from_query_stranded('chr1', *query, strand, features) + assert [x.id for x in result] == expected + + +@pytest.mark.parametrize('strand,features_on_plus,expected', [ + (None, True, ['plus', 'minus']), ('+', True, ['plus']), ('+', False, ['minus']), +]) +def test_unstranded_queries_prefer_requested_strand_then_fall_back(strand, features_on_plus, expected): + features = defaultdict(list) + a = feature(10, 20, '+'); a.id = 'plus' + b = feature(10, 20, '-'); b.id = 'minus' + if features_on_plus: features['chr1', 0, '+'] = [a] + features['chr1', 0, '-'] = [b] + result = ab.get_all_overlapping_features_from_query('chr1', 9, 20, strand, features, stranded=False) + assert [x.id for x in result] == expected + + +@pytest.mark.parametrize('delimiter,text', [(',', 'protein_coding,CDS\nmiRNA,exon\n'), + ('\t', 'protein_coding\tCDS\nmiRNA\texon\n')]) +def test_priority_file_preserves_class_and_region_order(tmp_path, delimiter, text): + from annotator.annotator import parse_annotation_priority + path = tmp_path / 'priority.txt' + path.write_text(text) + assert parse_annotation_priority(str(path), delimiter) == [['protein_coding','CDS'], ['miRNA','exon']] + + +def test_profiling_decorator_preserves_arguments_result_and_prints_stats(capsys): + @ab.do_cprofile + def add(left, *, right): + return left + right + assert add(2, right=3) == 5 + assert 'function calls' in capsys.readouterr().out + + +def test_profiling_decorator_propagates_failure_and_prints_stats(capsys): + @ab.do_cprofile + def failing(): + raise ValueError('sentinel') + with pytest.raises(ValueError, match='sentinel'): + failing() + assert 'function calls' in capsys.readouterr().out diff --git a/unit_tests/test_identifier_contracts.py b/unit_tests/test_identifier_contracts.py new file mode 100644 index 0000000..9ce22d4 --- /dev/null +++ b/unit_tests/test_identifier_contracts.py @@ -0,0 +1,128 @@ +"""Offline GFF and tabular fixtures for gene and microRNA identifier conversion.""" + +import gffutils +import pandas as pd +import pytest + +from annotator import gene_name2id as genes +from annotator import miRNA_name2id as mirnas + + +@pytest.fixture +def feature_db(tmp_path): + path = tmp_path / 'features.db' + gffutils.create_db('chr1\tx\tgene\t1\t50\t.\t+\t.\tID=g1;gene_id=g1;gene_name=Alpha\n' + 'chr1\tx\tgene\t70\t90\t.\t+\t.\tID=g2;gene_id=g2;gene_name=Alpha\n' + 'chr1\tx\tmiRNA_primary_transcript\t100\t160\t.\t+\t.\tID=pre1;Name=mir-one\n' + 'chr1\tx\tmiRNA_primary_transcript\t170\t200\t.\t-\t.\tID=pre2;Name=mir-one\n' + 'chr1\tx\tmiRNA\t110\t130\t.\t+\t.\tID=mat1;Derives_from=pre1\n' + 'chr1\tx\tmiRNA\t140\t150\t.\t+\t.\tID=mat2;Derives_from=pre1\n', + dbfn=str(path), from_string=True) + return str(path) + + +def test_gene_database_lookup_accumulates_ids_for_shared_name(feature_db): + assert dict(genes.build_name_to_id_dict(feature_db, None)) == {'Alpha': ['g1', 'g2']} + + +def test_gene_custom_mapping_extends_database_and_ignores_comments(feature_db, tmp_path): + custom = tmp_path / 'names.tsv' + custom.write_text('# ID\tname\ng3\tAlpha\ng4\tBeta\n') + assert dict(genes.build_name_to_id_dict(feature_db, str(custom))) == { + 'Alpha': ['g1', 'g2', 'g3'], 'Beta': ['g4']} + + +def test_precursor_database_lookup_accumulates_shared_names(feature_db): + assert dict(mirnas.build_name_to_precursor_id_dict(feature_db, None)) == {'mir-one': ['pre1', 'pre2']} + + +def test_precursor_custom_mapping_uses_first_two_columns(tmp_path): + custom = tmp_path / 'mirnas.tsv' + custom.write_text('# comment\npre1\tmir-one\textra\npre2\tmir-one\n') + assert dict(mirnas.build_name_to_precursor_id_dict(None, str(custom))) == {'mir-one': ['pre1', 'pre2']} + + +def test_mature_database_lookup_groups_ids_by_parent(feature_db): + assert dict(mirnas.build_precursor_id_to_mature_id_dict(feature_db, None)) == {'pre1': ['mat1', 'mat2']} + + +def test_mature_custom_mapping_splits_pipe_delimited_ids(tmp_path): + custom = tmp_path / 'mirnas.tsv' + custom.write_text('# comment\npre1\tmir-one\tmat1|mat2\n') + assert dict(mirnas.build_precursor_id_to_mature_id_dict(None, str(custom))) == {'pre1': ['mat1', 'mat2']} + + +@pytest.mark.parametrize('module,function', [(genes, 'build_name_to_id_dict'), + (mirnas, 'build_name_to_precursor_id_dict'), + (mirnas, 'build_precursor_id_to_mature_id_dict')]) +def test_identifier_mapping_without_sources_is_empty(module, function): + assert getattr(module, function)(None, None) == {} + + +@pytest.mark.parametrize('module,function', [(genes, 'get_id_from_namerow'), + (mirnas, 'get_precursor_id_from_namerow'), + (mirnas, 'get_mature_id_from_idrow')]) +def test_row_lookup_joins_ids_and_handles_unknown_name(module, function): + lookup = getattr(module, function) + assert lookup({'symbol': 'known'}, {'known': ['a', 'b']}, 'symbol') == 'a|b' + assert lookup({'symbol': 'missing'}, {'known': ['a']}, 'symbol') == '' + + +def test_gene_row_lookup_normalizes_numeric_names_to_strings(): + assert genes.get_id_from_namerow({'symbol': 123}, {'123': ['g1']}, 'symbol') == 'g1' + + +@pytest.mark.parametrize('module,function,column', [(genes, 'add_id_to_dataframe', 'gene id'), + (mirnas, 'add_precursor_id_to_dataframe', 'miRNA id'), + (mirnas, 'add_mature_id_to_dataframe', 'mature')]) +def test_dataframe_lookup_preserves_rows_and_unrelated_columns(module, function, column): + df = pd.DataFrame({'symbol': ['known', 'missing'], 'count': [7, 2]}) + assert getattr(module, function)(df, {'known': ['a', 'b']}, 'symbol') is df + assert df.to_dict('list') == {'symbol': ['known', 'missing'], 'count': [7, 2], column: ['a|b', '']} + + +@pytest.mark.parametrize('module,function,column', [(genes, 'add_id_to_dataframe', 'gene id'), + (mirnas, 'add_precursor_id_to_dataframe', 'miRNA id')]) +def test_dataframe_lookup_refuses_to_replace_existing_id_column(module, function, column): + df = pd.DataFrame({'symbol': ['known'], column: ['original']}) + with pytest.raises(AssertionError): + getattr(module, function)(df, {'known': ['new']}, 'symbol') + assert df[column].tolist() == ['original'] + + +@pytest.mark.parametrize('delete_original', [False, True]) +def test_gene_table_conversion_preserves_measurements_and_selects_index(tmp_path, delete_original): + source = tmp_path / 'counts.tsv'; source.write_text('symbol\tcount\nAlpha\t7\nBeta\t2\n') + mapping = tmp_path / 'names.tsv'; mapping.write_text('g1\tAlpha\ng2\tBeta\n') + df = genes.read_and_append_id(str(source), '\t', 'symbol', None, str(mapping), delete_original) + assert df.index.name == ('gene id' if delete_original else 'symbol') + assert df.index.tolist() == (['g1', 'g2'] if delete_original else ['Alpha', 'Beta']) + assert df['count'].tolist() == [7, 2] + out = tmp_path / 'out.tsv' + genes.convert(str(source), '\t', 'symbol', str(out), None, str(mapping), delete_original) + pd.testing.assert_frame_equal(pd.read_csv(out, sep='\t', index_col=0), df) + + +@pytest.mark.parametrize('delete_original', [False, True]) +def test_mirna_table_conversion_can_include_mature_ids(tmp_path, delete_original): + source = tmp_path / 'counts.tsv'; source.write_text('symbol\tcount\nmir-one\t7\nmir-two\t2\n') + mapping = tmp_path / 'names.tsv'; mapping.write_text('pre1\tmir-one\tmat1|mat2\npre2\tmir-two\tmat3\n') + df = mirnas.read_and_append_id(str(source), '\t', 'symbol', None, str(mapping), True, delete_original) + assert df.index.name == ('miRNA id' if delete_original else 'symbol') + assert df.index.tolist() == (['pre1', 'pre2'] if delete_original else ['mir-one', 'mir-two']) + assert df['mature'].tolist() == ['mat1|mat2', 'mat3'] + assert df['count'].tolist() == [7, 2] + out = tmp_path / 'out.tsv' + mirnas.convert(str(source), '\t', 'symbol', str(out), None, str(mapping), True, delete_original) + pd.testing.assert_frame_equal(pd.read_csv(out, sep='\t', index_col=0), df) + + +@pytest.mark.parametrize('module,function', [(genes, 'build_name_to_id_dict'), + (mirnas, 'build_name_to_precursor_id_dict'), + (mirnas, 'build_precursor_id_to_mature_id_dict')]) +def test_mapping_sources_report_missing_paths_without_losing_available_custom_data(module, function, tmp_path, capsys): + missing = tmp_path / 'missing.db' + assert getattr(module, function)(str(missing), str(tmp_path / 'missing.tsv')) == {} + out = capsys.readouterr().out + assert 'Could not use the gffdb file' in out + assert 'Could not use the custom file' in out diff --git a/unit_tests/test_region_contracts.py b/unit_tests/test_region_contracts.py new file mode 100644 index 0000000..8cfccbc --- /dev/null +++ b/unit_tests/test_region_contracts.py @@ -0,0 +1,173 @@ +"""Small real BED/GTF fixtures for region exports and transcript geometry.""" + +import gffutils +import pybedtools +import pytest + +from annotator import annotation_functions as af +from annotator import create_region_bedfiles as regions +from annotator import create_AS_STRUCTURE as structure + + +@pytest.fixture +def db(): + rows = [('gene', 1, 1600), ('transcript', 1, 1600), + ('exon', 1, 100), ('exon', 1501, 1600), + ('CDS', 20, 100), ('CDS', 1501, 1550), + ('UTR', 1, 19), ('UTR', 1551, 1600)] + text = ''.join('chr1\tx\t{}\t{}\t{}\t.\t+\t.\tgene_id "g1"; gene_name "Alpha"; transcript_id "tx1";\n'.format(*r) for r in rows) + return gffutils.create_db(text, dbfn=':memory:', from_string=True, + disable_infer_genes=True, disable_infer_transcripts=True) + + +@pytest.fixture +def keys(): + return af.get_keys('gencode') + + +def coords(tool): + return [(x.chrom, x.start, x.end, x.name, x.strand) for x in tool] + + +@pytest.mark.parametrize('by_transcript,name', [(False, 'g1'), (True, 'tx1')]) +def test_create_bedtools_converts_inclusive_gtf_to_half_open_bed(db, keys, by_transcript, name): + assert coords(regions.create_bedtools(list(db.features_of_type('exon')), keys, by_transcript)) == [ + ('chr1', 0, 100, name, '+'), ('chr1', 1500, 1600, name, '+')] + + +def test_create_bedtools_retains_each_associated_id(keys): + f = gffutils.Feature(seqid='chr2', start=5, end=5, strand='-', + attributes={'gene_id': ['a', 'b']}) + assert coords(regions.create_bedtools([f], keys)) == [ + ('chr2', 4, 5, 'a', '-'), ('chr2', 4, 5, 'b', '-')] + + +def test_merge_bedtool_by_gene_preserves_gene_and_strand_boundaries(): + bed = pybedtools.BedTool('chr1\t10\t20\ta\t0\t+\nchr1\t15\t30\ta\t0\t+\n' + 'chr1\t15\t40\tb\t0\t+\nchr1\t18\t25\ta\t0\t-\n' + 'chr2\t10\t20\ta\t0\t+\n', from_string=True) + df = regions.merge_bedtool_by_gene(bed) + actual = set(df[['chrom', 'start', 'end', 'name', 'strand']].itertuples(index=False, name=None)) + assert actual == {('chr1', 10, 30, 'a', '+'), ('chr1', 15, 40, 'b', '+'), + ('chr1', 18, 25, 'a', '-'), ('chr2', 10, 20, 'a', '+')} + + +@pytest.mark.parametrize('by_transcript,name', [(False, 'g1'), (True, 'tx1')]) +def test_cds_region_export_and_returned_bed_have_same_coordinates(db, keys, tmp_path, by_transcript, name): + out = tmp_path / 'cds.bed' + result = regions.create_cds_region_bedfile(db, keys, str(out), by_transcript) + expected = [('chr1', 19, 100, name, '+'), ('chr1', 1500, 1550, name, '+')] + assert coords(result) == expected + assert coords(pybedtools.BedTool(str(out))) == expected + + +@pytest.mark.parametrize('by_transcript,name', [(False, 'g1'), (True, 'tx1')]) +def test_exon_region_export_preserves_span_and_selected_name(db, keys, tmp_path, by_transcript, name): + out = tmp_path / 'exon.bed' + assert regions.create_exon_region_bedfile(db, keys, str(out), by_transcript) is None + assert coords(pybedtools.BedTool(str(out))) == [('chr1', 0, 100, name, '+'), ('chr1', 1500, 1600, name, '+')] + + +def test_gene_region_export_uses_gene_identifier(db, keys, tmp_path): + out = tmp_path / 'gene.bed' + regions.create_gene_region_bedfile(db, keys, str(out)) + assert coords(pybedtools.BedTool(str(out))) == [('chr1', 0, 1600, 'g1', '+')] + + +@pytest.mark.parametrize('by_transcript,name', [(False, 'g1'), (True, 'tx1')]) +def test_utr_exports_classify_generic_utr_features(db, keys, tmp_path, by_transcript, name): + five = tmp_path / '5utr.bed'; three = tmp_path / '3utr.bed' + result5, result3 = regions.create_utr_region_bedfiles( + db, keys, af.get_all_cds_dict(db, 'CDS'), str(three), str(five), by_transcript) + assert coords(result5) == [('chr1', 0, 19, name, '+')] + assert coords(result3) == [('chr1', 1550, 1600, name, '+')] + assert coords(pybedtools.BedTool(str(five))) == coords(result5) + assert coords(pybedtools.BedTool(str(three))) == coords(result3) + + +def test_utr_exports_accept_explicit_ensembl_utr_types(keys): + f5 = gffutils.Feature(seqid='chr2', start=1, end=10, strand='-', attributes={'transcript_id':['tx'], 'gene_id':['g']}) + f3 = gffutils.Feature(seqid='chr2', start=80, end=100, strand='-', attributes={'transcript_id':['tx'], 'gene_id':['g']}) + class DB: + def features_of_type(self, kind): + return iter({'five_prime_utr': [f5], 'three_prime_utr': [f3]}.get(kind, [])) + five, three = regions.create_utr_region_bedfiles(DB(), af.get_keys('ensembl'), {}, None, None, True) + assert coords(five) == [('chr2', 0, 10, 'tx', '-')] + assert coords(three) == [('chr2', 79, 100, 'tx', '-')] + + +def test_utr_exports_report_unsupported_vocabulary(capsys): + assert regions.create_utr_region_bedfiles(None, af.get_keys('refseq'), {}, None, None, True) == 1 + assert 'Having trouble guessing' in capsys.readouterr().out + + +@pytest.mark.parametrize('by_transcript,name', [(False, 'g1'), (True, 'tx1')]) +def test_intron_return_values_partition_inclusive_exon_gap(db, keys, by_transcript, name): + exons = af.get_all_exons_dict(db, 'exon', 'transcript_id') + transcripts = af.get_all_transcripts_dict(db, 'transcript', 'transcript_id') + prox, dist, all_introns = regions.create_intron_region_bedfiles( + db, exons, transcripts, keys, None, None, None, by_transcript) + assert coords(prox) == [('chr1', 100, 600, name, '+'), ('chr1', 1000, 1500, name, '+')] + assert coords(dist) == [('chr1', 600, 1000, name, '+')] + assert coords(all_introns) == [('chr1', 100, 1500, name, '+')] + + +def test_gene_intron_files_match_returned_regions(db, keys, tmp_path): + paths = [tmp_path / n for n in ['prox.bed', 'dist.bed', 'all.bed']] + tools = regions.create_intron_region_bedfiles( + db, af.get_all_exons_dict(db, 'exon', 'transcript_id'), + af.get_all_transcripts_dict(db, 'transcript', 'transcript_id'), keys, + *(str(p) for p in paths), False) + assert [coords(pybedtools.BedTool(str(p))) for p in paths] == [coords(t) for t in tools] + + +def test_transcript_geometry_chooses_first_longest_and_outer_bounds(): + genes = {'g': ['short', 'long', 'tied']} + transcripts = {'short': {'start': 10, 'end': 20}, 'long': {'start': 5, 'end': 50}, + 'tied': {'start': 7, 'end': 52}} + assert af.get_longest_transcripts(genes, transcripts) == {'g': 'long'} + outer = af.most_upstream_downstream_positions(genes, transcripts) + assert outer == {'g': {'start': 5, 'end': 52}} + assert af.get_premrna_lengths('g', outer) == 48 + + +def test_mrna_length_counts_exon_union_across_isoforms(): + exon = lambda start,end: {'chrom':'chr1', 'start': start, 'end':end, 'strand':'+'} + assert af.get_mrna_lengths('g', {'tx1':[exon(1,10),exon(31,40)], + 'tx2':[exon(5,15),exon(31,40)]}, {'g':['tx1','tx2']}) == 25 + + +def test_as_structure_export_includes_mrna_and_premrna_lengths(db, tmp_path): + out = tmp_path / 'structure.gff' + structure.create_as_structure(db, 'gencode', str(out)) + fields = out.read_text().strip().split('\t') + assert fields[:8] == ['chr1','AS_STRUCTURE','gene','1','1600','.','+','.'] + assert fields[8] == 'gene_id=g1;mrna_length=200;premrna_length=1600' + + +def test_total_cds_length_merges_overlapping_coding_bases(db, keys, capsys): + from annotator import get_region_lengths as lengths + assert lengths.calculate_total_cds_length(db, keys) == 131 + assert 'cds lengths (total): 131' in capsys.readouterr().out + + +def test_total_intron_length_reports_each_partition(db, keys, capsys): + from annotator import get_region_lengths as lengths + lengths.calculate_total_intron_lengths( + db, af.get_all_exons_dict(db, 'exon', 'transcript_id'), + af.get_all_transcripts_dict(db, 'transcript', 'transcript_id'), keys) + out = capsys.readouterr().out + assert 'proximal intron lengths (total): 1000' in out + assert 'distal intron lengths (total): 400' in out + assert 'all intron lengths (total): 1400' in out + + +def test_total_lengths_wrapper_reports_actual_coding_utr_and_intron_bases(db, tmp_path, capsys): + from annotator import get_region_lengths as lengths + path = tmp_path / 'database.db' + db.conn.backup(__import__('sqlite3').connect(str(path))) + assert lengths.calculate_total_lengths(str(path), 'gencode') is None + out = capsys.readouterr().out + for text in ['cds lengths (total): 131', 'five_prime_utr_lengths (total): 19', + 'three_prime_utr_lengths (total): 50', 'all intron lengths (total): 1400']: + assert text in out From 5bc1437dba92f043ca13cf501d6a67554dc68019 Mon Sep 17 00:00:00 2001 From: Brian Yee Date: Mon, 12 Oct 2026 22:40:22 -0700 Subject: [PATCH 2/2] test: cover advanced annotation CLI and boundary cases --- unit_tests/test_cli_contracts.py | 27 +++++++++++++++++++++++++ unit_tests/test_coordinate_contracts.py | 17 ++++++++++++++++ unit_tests/test_region_contracts.py | 13 ++++++++++++ 3 files changed, 57 insertions(+) diff --git a/unit_tests/test_cli_contracts.py b/unit_tests/test_cli_contracts.py index ab1e02f..b964dcf 100644 --- a/unit_tests/test_cli_contracts.py +++ b/unit_tests/test_cli_contracts.py @@ -92,3 +92,30 @@ def test_cli_help_lists_options_and_exits_successfully(module, monkeypatch, caps module.main() assert exc.value.code == 0 assert 'usage:' in capsys.readouterr().out + + +def test_annotation_cli_applies_prefix_strand_fallback_and_priority_files(tmp_path, monkeypatch): + gtf=tmp_path/'input.gtf' + gtf.write_text(''.join('1\tx\t{}\t{}\t{}\t.\t-\t.\tgene_id "g1"; gene_name "Alpha"; transcript_id "tx1"; transcript_type "protein_coding"; gene_type "protein_coding";\n'.format(*r) for r in [ + ('gene',1,100),('transcript',1,100),('exon',1,100),('CDS',20,80)])) + db=tmp_path/'input.db';build_gffutils_db.build_db(str(gtf),str(db),True,True,True) + source=tmp_path/'queries.bed';source.write_text('chr1\t30\t40\tquery\t0\t+\n') + priority=tmp_path/'priority.txt';priority.write_text('protein_coding,CDS\nprotein_coding,exon\n') + out=tmp_path/'output.bed' + monkeypatch.setattr('sys.argv',['annotator','--input',str(source),'--output',str(out),'--gtfdb',str(db), + '--append-chr','--unstranded','--limit-chroms-to','1', + '--transcript-priority-file',str(priority),'--gene-priority-file',str(priority)]) + annotator.main() + assert out.read_text().strip().split('\t')[6:9]==['g1','Alpha','CDS'] + + +@pytest.mark.parametrize('strand,expected', [('+','miRNA'),('-','miRNA')]) +def test_annotation_cli_identifies_mirna_exons(tmp_path,monkeypatch,strand,expected): + gtf=tmp_path/'input.gtf' + gtf.write_text(''.join('chr1\tx\t{}\t1\t100\t.\t{}\t.\tgene_id "g1"; gene_name "mir-one"; transcript_id "tx1"; transcript_type "miRNA"; gene_type "miRNA";\n'.format(kind,strand) for kind in ['gene','transcript','exon'])) + db=tmp_path/'input.db';build_gffutils_db.build_db(str(gtf),str(db),True,True,True) + source=tmp_path/'queries.bed';source.write_text('chr1\t30\t40\tquery\t0\t'+strand+'\n') + out=tmp_path/'output.bed' + monkeypatch.setattr('sys.argv',['annotator','--input',str(source),'--output',str(out),'--gtfdb',str(db)]) + annotator.main() + assert out.read_text().strip().split('\t')[6:9]==['g1','mir-one',expected] diff --git a/unit_tests/test_coordinate_contracts.py b/unit_tests/test_coordinate_contracts.py index 02b9fd8..01cf20c 100644 --- a/unit_tests/test_coordinate_contracts.py +++ b/unit_tests/test_coordinate_contracts.py @@ -185,3 +185,20 @@ def failing(): with pytest.raises(ValueError, match='sentinel'): failing() assert 'function calls' in capsys.readouterr().out + + +@pytest.mark.parametrize('query', [(0,15),(15,30)]) +def test_stranded_overlap_finds_features_when_only_one_query_edge_overlaps(query): + f=feature(10,20);f.id='partial' + features=defaultdict(list);features['chr1',0,'+']=[f] + result=ab.get_all_overlapping_features_from_query_stranded('chr1',*query,'+',features) + assert [x.id for x in result]==['partial'] + + +def test_gene_name_lookup_reports_missing_required_attribute(capsys): + gene=SimpleNamespace(attributes={'gene_id':['g1']}) + db=SimpleNamespace(features_of_type=lambda kind:iter([gene])) + with pytest.raises(SystemExit) as exc: + af.gene_id_to_name(db,'gene_name') + assert exc.value.code==1 + assert 'Warning. Key not found' in capsys.readouterr().out diff --git a/unit_tests/test_region_contracts.py b/unit_tests/test_region_contracts.py index 8cfccbc..a47d07f 100644 --- a/unit_tests/test_region_contracts.py +++ b/unit_tests/test_region_contracts.py @@ -171,3 +171,16 @@ def test_total_lengths_wrapper_reports_actual_coding_utr_and_intron_bases(db, tm for text in ['cds lengths (total): 131', 'five_prime_utr_lengths (total): 19', 'three_prime_utr_lengths (total): 50', 'all intron lengths (total): 1400']: assert text in out + + +def test_gene_export_preserves_distinct_features_with_repeated_gene_id(keys,tmp_path,capsys): + genes=[gffutils.Feature(seqid='chr1',start=1,end=10,strand='+',attributes={'gene_id':['same']}), + gffutils.Feature(seqid='chr2',start=20,end=30,strand='-',attributes={'gene_id':['same']})] + class DB: + def features_of_type(self,kind): + assert kind=='gene' + return iter(genes) + out=tmp_path/'genes.bed' + regions.create_gene_region_bedfile(DB(),keys,str(out)) + assert coords(pybedtools.BedTool(str(out)))==[('chr1',0,10,'same','+'),('chr2',19,30,'same','-')] + assert 'hashkey exists for gene same' in capsys.readouterr().out