diff --git a/.github/workflows/check-apis.yml b/.github/workflows/check-apis.yml new file mode 100644 index 0000000..c44cdd0 --- /dev/null +++ b/.github/workflows/check-apis.yml @@ -0,0 +1,42 @@ +name: check-apis + +# A weekly canary against the live IMGT and OGRDB APIs. It runs the same code a +# user's download would, using the endpoint constants the sources define, so a +# red run here is the early warning that an upstream API changed shape. It never +# blocks a merge; it only watches. workflow_dispatch allows an on-demand check. +on: + schedule: + - cron: "0 6 * * 1" # Mondays, 06:00 UTC + workflow_dispatch: + +concurrency: + group: check-apis + cancel-in-progress: true + +jobs: + probe: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: + python-version: '3.12' + - name: Install package + run: | + python -m pip install --upgrade pip + pip install . + - name: Probe the reference APIs + env: + SOURCERER_LIVE: '1' + run: python -m unittest tests.test_live -v + # On failure the run goes red and GitHub notifies the watchers. Opening an + # issue automatically is intentionally left off; enable the step below if a + # tracked issue is wanted instead of (or as well as) the email. + # + # - name: Open an issue on failure + # if: failure() + # uses: JasonEtco/create-an-issue@v2 + # env: + # GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + # with: + # filename: .github/api-breakage-issue.md diff --git a/NEWS.rst b/NEWS.rst index 3727502..7117350 100644 --- a/NEWS.rst +++ b/NEWS.rst @@ -26,4 +26,24 @@ General: + Added a download provenance record (what was fetched, from where, when, and its hash) written alongside every download. +Germline references: + ++ Added germline reference sources IMGT_ (``sourcerer imgt``) and OGRDB_ + (``sourcerer ogrdb``, also reachable as ``sourcerer airrc``), and an + ``airrc-imgt`` blend that takes immunoglobulin V, D and J from OGRDB's AIRR-C + sets and the T-cell receptor and remaining constants from IMGT. ++ Added ``sourcerer download ``, which writes the germline + ``reference_base`` in the `nf-core/airrflow`_ layout, and ``--igblast`` to also + build the IgBLAST databases (``makeblastdb`` plus the NCBI internal_data and + optional_file trees). ++ Added ``sourcerer reference build``, to validate a germline reference folder + -- in the ``reference_base`` layout or a flat folder of FASTAs -- and build + its IgBLAST databases, with ``--check`` to validate without building. ++ nf-core/airrflow can fetch germlines through ``sourcerer`` for its ``imgt`` + and ``airrc-imgt`` database types, or consume a ``sourcerer``-built reference + passed to ``--reference_fasta`` / ``--reference_igblast``. + .. _OAS: https://opig.stats.ox.ac.uk/webapps/oas/ +.. _IMGT: https://www.imgt.org/genedb/ +.. _OGRDB: https://ogrdb.airr-community.org/ +.. _nf-core/airrflow: https://nf-co.re/airrflow diff --git a/README.rst b/README.rst index b25b3c0..f042b1f 100644 --- a/README.rst +++ b/README.rst @@ -2,11 +2,21 @@ sourcerer ================================================================================ ``sourcerer`` downloads data from online immune repertoire databases and formats -it for use with the Immcantation_ framework. Each external source is a module; -the first is OAS_ (Observed Antibody Space). +it for use with the Immcantation_ framework and `nf-core/airrflow`_. Each +external source is a module, and sources come in two kinds: + +- *dataset* sources such as OAS_ (Observed Antibody Space) download sequencing + data and write an airrflow samplesheet; +- *germline reference* sources -- IMGT_, OGRDB_, and an ``airrc-imgt`` blend of + the two -- download germline sets and build the ``reference_base`` and IgBLAST + databases airrflow consumes. ``sourcerer reference`` can also validate and + build those databases from a reference folder you already have. .. _Immcantation: https://immcantation.readthedocs.io +.. _nf-core/airrflow: https://nf-co.re/airrflow .. _OAS: https://opig.stats.ox.ac.uk/webapps/oas/ +.. _IMGT: https://www.imgt.org/genedb/ +.. _OGRDB: https://ogrdb.airr-community.org/ Why -------------------------------------------------------------------------------- @@ -25,6 +35,8 @@ a reviewable diff and a failing test, not as a silently wrong download. Usage -------------------------------------------------------------------------------- +Datasets (OAS), producing an airrflow samplesheet: + .. code-block:: bash sourcerer --version @@ -38,6 +50,25 @@ Usage --outdir airrflow_out -c ../airrflow.config \ --clonal_threshold 0.2 -resume +Germline references, producing the ``reference_base`` and IgBLAST databases: + +.. code-block:: bash + + # IMGT germline for a species, and (with --igblast) the IgBLAST databases + sourcerer imgt download human --outdir ref --igblast + + # the AIRR-C sets blended with IMGT (immunoglobulin from OGRDB, TR and the + # remaining constants from IMGT) -- the airrflow airrc-imgt reference + sourcerer airrc-imgt download human --outdir ref --igblast + + # validate a germline folder someone provided and build its databases; + # --check validates only, without makeblastdb + sourcerer reference build ref/reference_base --out igblast_base --check + +nf-core/airrflow uses the result either way: point ``--reference_fasta`` and +``--reference_igblast`` at ``ref/reference_base`` and ``igblast_base`` with +``--fetch_germlines none``. + License -------------------------------------------------------------------------------- diff --git a/docs/api.rst b/docs/api.rst index 4e97ed4..cb808a1 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -13,9 +13,14 @@ API modules/Catalog modules/Convert modules/Airrflow + modules/Reference modules/Provenance modules/Gzip modules/Exceptions modules/Sources modules/SourcesBase + modules/SourcesGermline modules/SourcesOas + modules/SourcesImgt + modules/SourcesOgrdb + modules/SourcesAirrcImgt diff --git a/docs/info.rst b/docs/info.rst index 0e00c4b..d3d9b7c 100644 --- a/docs/info.rst +++ b/docs/info.rst @@ -37,6 +37,32 @@ asks that both of the following be cited: *Protein Sci*. 2022;31(1):141-146. doi:`10.1002/pro.4205 `__ +IMGT +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Data from `IMGT `__ is governed by the `IMGT terms +of use `__: it is free for academic +research on condition that IMGT is cited. IMGT asks that the following be cited: + +- Lefranc MP, Giudicelli V, Duroux P, et al. IMGT, the international + ImMunoGeneTics information system 25 years on. *Nucleic Acids Res*. + 2015;43(Database issue):D413-D422. + doi:`10.1093/nar/gku1056 `__ + +The ``airrc-imgt`` blend uses IMGT data too, so its downloads carry this +obligation as well as the OGRDB one below. + +OGRDB +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Data from `OGRDB `__ is distributed under a +`CC BY 4.0 `__ license. In +exchange, OGRDB asks that the following be cited: + +- Lees WD, Busse CE, Corcoran M, et al. OGRDB: a reference database of inferred + immune receptor genes. *Nucleic Acids Res*. 2020;48(D1):D964-D970. + doi:`10.1093/nar/gkz822 `__ + License -------------------------------------------------------------------------------- diff --git a/docs/modules/Reference.rst b/docs/modules/Reference.rst new file mode 100644 index 0000000..11e90f8 --- /dev/null +++ b/docs/modules/Reference.rst @@ -0,0 +1,7 @@ +sourcerer.Reference +------------------- + +.. automodule:: sourcerer.Reference + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/modules/SourcesAirrcImgt.rst b/docs/modules/SourcesAirrcImgt.rst new file mode 100644 index 0000000..ecf6819 --- /dev/null +++ b/docs/modules/SourcesAirrcImgt.rst @@ -0,0 +1,7 @@ +sourcerer.Sources.AirrcImgt +--------------------------- + +.. automodule:: sourcerer.Sources.AirrcImgt + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/modules/SourcesGermline.rst b/docs/modules/SourcesGermline.rst new file mode 100644 index 0000000..773b5d3 --- /dev/null +++ b/docs/modules/SourcesGermline.rst @@ -0,0 +1,7 @@ +sourcerer.Sources.Germline +-------------------------- + +.. automodule:: sourcerer.Sources.Germline + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/modules/SourcesImgt.rst b/docs/modules/SourcesImgt.rst new file mode 100644 index 0000000..7209228 --- /dev/null +++ b/docs/modules/SourcesImgt.rst @@ -0,0 +1,7 @@ +sourcerer.Sources.Imgt +---------------------- + +.. automodule:: sourcerer.Sources.Imgt + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/modules/SourcesOgrdb.rst b/docs/modules/SourcesOgrdb.rst new file mode 100644 index 0000000..f9ffd94 --- /dev/null +++ b/docs/modules/SourcesOgrdb.rst @@ -0,0 +1,7 @@ +sourcerer.Sources.Ogrdb +----------------------- + +.. automodule:: sourcerer.Sources.Ogrdb + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/usage/airrc-imgt.rst b/docs/usage/airrc-imgt.rst new file mode 100644 index 0000000..504ca8c --- /dev/null +++ b/docs/usage/airrc-imgt.rst @@ -0,0 +1,15 @@ +.. _UsageAirrcImgt: + +sourcerer airrc-imgt +================================================================================ + +The AIRR-C germline sets blended with IMGT: immunoglobulin V, D and J from +OGRDB, and everything OGRDB does not cover -- all of the T-cell receptor, and +the immunoglobulin constants without a published set -- from IMGT. Offers +``human`` and ``mouse`` collections. ``download`` writes an airrflow +``reference_base`` mixing ``airrc_`` and ``imgt_`` files; ``--igblast`` +additionally builds the IgBLAST databases. + +.. autoprogram:: sourcerer.Cli:getArgParser() + :prog: sourcerer + :start_command: airrc-imgt diff --git a/docs/usage/imgt.rst b/docs/usage/imgt.rst new file mode 100644 index 0000000..c585e57 --- /dev/null +++ b/docs/usage/imgt.rst @@ -0,0 +1,14 @@ +.. _UsageImgt: + +sourcerer imgt +================================================================================ + +`IMGT/GENE-DB `__: germline V, D, J and C +reference sequences. Offers ``human`` and ``mouse`` collections, each narrowed +with the ``--locus`` and ``--segment`` filters below. ``download`` writes an +airrflow ``reference_base``; ``--igblast`` additionally builds the IgBLAST +databases, which needs ``makeblastdb`` on the path. + +.. autoprogram:: sourcerer.Cli:getArgParser() + :prog: sourcerer + :start_command: imgt diff --git a/docs/usage/index.rst b/docs/usage/index.rst index 36fd112..b821d4a 100644 --- a/docs/usage/index.rst +++ b/docs/usage/index.rst @@ -4,17 +4,26 @@ Commandline Usage ================================================================================ ``sourcerer`` is a single command with a subcommand tree: one subcommand per -external source (``oas`` today), a ``schema`` subcommand to inspect and -re-harvest the stored snapshot each source is built from, and a ``sources`` -subcommand that lists what is registered. +external source (``oas``, ``imgt``, ``ogrdb`` and ``airrc-imgt`` today), a +``schema`` subcommand +to inspect and re-harvest the stored snapshot each source is built from, a +``reference`` subcommand to validate a germline reference folder and build its +IgBLAST databases, and a ``sources`` subcommand that lists what is registered. + +Sources come in two kinds. A repertoire source such as ``oas`` downloads +sequencing data and writes an airrflow samplesheet; a germline reference source +such as ``imgt`` and ``ogrdb`` downloads germline sets and writes an airrflow +``reference_base``, optionally building the IgBLAST databases with ``--igblast``. +The ``reference`` subcommand does that same build for a reference folder supplied +by hand, in either the ``reference_base`` layout or a flat folder of FASTAs. Every source subcommand exposes the same two actions, ``search`` and -``download``, each taking a collection (for example ``paired`` or -``unpaired``) as a further subcommand. The filter flags under a collection — -``--species``, ``--disease``, and so on — are not hardcoded: they are -generated at parser-construction time from the checked-in schema snapshot -described in :ref:`API`, which is also why they appear below exactly as they -would in ``--help`` on the machine building these docs. +``download``, each taking a collection (for example ``paired`` and ``unpaired`` +for OAS, or a species for the germline sources) as a further subcommand. The +filter flags under a collection — ``--species``, ``--locus``, and so on — are +not hardcoded: they are generated at parser-construction time from the checked-in +schema snapshot described in :ref:`API`, which is also why they appear below +exactly as they would in ``--help`` on the machine building these docs. Commands are documented one page per top level subcommand, mirroring how the commandline itself groups them: :doc:`sources` and :doc:`schema` apply @@ -27,4 +36,8 @@ one page here alongside it. sources schema + reference oas + imgt + ogrdb + airrc-imgt diff --git a/docs/usage/ogrdb.rst b/docs/usage/ogrdb.rst new file mode 100644 index 0000000..8d1cec6 --- /dev/null +++ b/docs/usage/ogrdb.rst @@ -0,0 +1,19 @@ +.. _UsageOgrdb: + +sourcerer ogrdb +================================================================================ + +`OGRDB `__: AIRR Community curated +immunoglobulin germline sets. Offers ``human`` and ``mouse`` collections, +narrowed with the ``--locus`` filter below. ``download`` writes an airrflow +``reference_base``; ``--igblast`` additionally builds the IgBLAST databases, +which needs ``makeblastdb`` on the path. + +OGRDB is the AIRR Community database, so ``ogrdb`` also answers to the alias +``airrc`` (``sourcerer airrc download ...``). For a reference that additionally +fills in the T-cell receptor and the remaining constants from IMGT, use the +``airrc-imgt`` source instead. + +.. autoprogram:: sourcerer.Cli:getArgParser() + :prog: sourcerer + :start_command: ogrdb diff --git a/docs/usage/reference.rst b/docs/usage/reference.rst new file mode 100644 index 0000000..b107664 --- /dev/null +++ b/docs/usage/reference.rst @@ -0,0 +1,16 @@ +.. _UsageReference: + +sourcerer reference +================================================================================ + +Validate a folder of germline FASTAs and build the IgBLAST databases from it, +for a reference someone supplies rather than one sourcerer downloaded. Files are +recognised by name in any directory layout -- +``[_][aa_]_.fasta``, for example ``human_IGHV.fasta`` or +``imgt_human_IGHV.fasta`` -- so a nested ``reference_base`` and a flat folder both +work. ``--check`` validates and reports what would build without building +anything, and needs no ``makeblastdb``. + +.. autoprogram:: sourcerer.Cli:getArgParser() + :prog: sourcerer + :start_command: reference diff --git a/src/sourcerer/Cli.py b/src/sourcerer/Cli.py index 255b702..12d2bc4 100644 --- a/src/sourcerer/Cli.py +++ b/src/sourcerer/Cli.py @@ -15,13 +15,13 @@ from pathlib import Path # Sourcerer imports -from sourcerer import Catalog, Convert, Provenance +from sourcerer import Catalog, Convert, Provenance, Reference from sourcerer.Airrflow import buildSamplesheet from sourcerer.Commandline import CommonHelpFormatter, setupLogging from sourcerer.Exceptions import SourcererError from sourcerer.Http import HttpClient from sourcerer.Schema import loadSchema, saveSchema -from sourcerer.Sources import REGISTRY, getSource +from sourcerer.Sources import ALIASES, REGISTRY, canonicalName, getSource from sourcerer.Version import __date__, __version__ log = logging.getLogger('sourcerer') @@ -144,12 +144,49 @@ def getArgParser(): formatter_class=CommonHelpFormatter) _addSchemaParser(commands) + _addReferenceParser(commands) for name, source in sorted(REGISTRY.items()): _addSourceParser(commands, name, source) return parser +def _addReferenceParser(commands): + """Add the reference subcommand: validate and build from a reference folder.""" + reference = commands.add_parser( + 'reference', help='validate a germline reference folder and build ' + 'IgBLAST databases from it', + description='Check that a folder of germline FASTAs is in a format ' + 'airrflow can use, and build the IgBLAST databases from it. ' + 'Files are recognised by name, in any directory layout: the ' + 'species and chain, with an optional source prefix and an ' + 'optional aa marker for translated V, as in human_IGHV.fasta ' + 'or imgt_human_IGHV.fasta.', + formatter_class=CommonHelpFormatter) + actions = reference.add_subparsers(dest='action', metavar='ACTION', + required=True) + + build = actions.add_parser( + 'build', help='validate a reference folder and build IgBLAST databases', + description='Validate the reference folder and build the IgBLAST ' + 'databases from it. With --check, only validate and report, ' + 'building nothing (and needing no makeblastdb).', + formatter_class=CommonHelpFormatter) + build.add_argument('folder', type=Path, + help='a reference_base tree or a flat folder of germline ' + 'FASTAs named _.fasta') + build.add_argument('--out', type=Path, default=None, + help='directory to write igblast_base into; required ' + 'unless --check') + build.add_argument('--check', action='store_true', + help='validate the folder and report what would build, ' + 'without building anything') + build.add_argument('--species', nargs='+', choices=list(Reference.SPECIES), + default=None, + help='limit to these species; default is every species ' + 'found in the folder') + + def _addSchemaParser(commands): """Add the schema subcommand tree.""" schema = commands.add_parser( @@ -167,8 +204,8 @@ def _addSchemaParser(commands): 'collection\'s fields and how many values each accepts, ' 'or every value a single field accepts.', formatter_class=CommonHelpFormatter) - show.add_argument('--source', required=True, choices=sorted(REGISTRY), - help='which source to read the snapshot of') + show.add_argument('--source', required=True, choices=sorted(REGISTRY) + sorted(ALIASES), + help='which source (or alias) to read the snapshot of') show.add_argument('--collection', default=None, help='list this collection\'s fields; without it, print ' 'one summary line per collection') @@ -184,7 +221,7 @@ def _addSchemaParser(commands): 'fetched detail-page enrichment is carried forward ' 'unless --refresh-details asks to redo it.', formatter_class=CommonHelpFormatter) - refresh.add_argument('--source', required=True, choices=sorted(REGISTRY), + refresh.add_argument('--source', required=True, choices=sorted(REGISTRY) + sorted(ALIASES), help='which source to contact and re-harvest') refresh.add_argument('--out', default=None, type=Path, help='directory to write into; without it the packaged ' @@ -205,7 +242,7 @@ def _addSourceParser(commands, name, source): schema = loadSchemaQuietly(name) parser = commands.add_parser( - name, help=source.description, + name, aliases=list(source.aliases), help=source.description, description='%s\n\nHomepage: %s' % (source.description, source.homepage), formatter_class=CommonHelpFormatter) actions = parser.add_subparsers(dest='action', metavar='ACTION', @@ -254,23 +291,34 @@ def _addSourceParser(commands, name, source): help='write the hits to a TSV file') else: leaf.add_argument('--outdir', type=Path, required=True, - help='directory to write into; one ' - 'subdirectory per format, plus a ' - 'samplesheet for each converted format') - leaf.add_argument('--format', action='append', dest='formats', - choices=FORMATS, - help='what to write, repeatable to write ' - 'several; raw mirrors the source files ' - 'untouched and is always written because ' - 'the others are converted from it, so ' - 'omitting this writes raw alone') + help='directory to write into') leaf.add_argument('--dry-run', action='store_true', help='report what would be fetched, then stop') leaf.add_argument('--no-resume', action='store_true', help='re-download in full rather than ' 'continuing a partly fetched file') - leaf.add_argument('--strict-airr', action='store_true', - help='drop columns the AIRR schema does not define') + if source.output == 'reference': + # Reference sources build a germline reference_base rather + # than converting to AIRR, so they take the igblast options + # instead of the format and AIRR-strictness ones. + leaf.add_argument('--igblast', action='store_true', + help='also build the IgBLAST databases ' + 'from the reference; needs makeblastdb ' + 'on PATH') + leaf.add_argument('--igblast-out', type=Path, default=None, + help='where to write igblast_base; ' + 'defaults to /igblast_base') + else: + leaf.add_argument('--format', action='append', dest='formats', + choices=FORMATS, + help='what to write, repeatable to write ' + 'several; raw mirrors the source files ' + 'untouched and is always written because ' + 'the others are converted from it, so ' + 'omitting this writes raw alone') + leaf.add_argument('--strict-airr', action='store_true', + help='drop columns the AIRR schema does ' + 'not define') def makeClient(args): @@ -282,6 +330,8 @@ def handleSources(args): """List registered sources.""" for name, source in sorted(REGISTRY.items()): print('%-10s %s' % (name, source.description)) + if source.aliases: + print('%-10s alias: %s' % ('', ', '.join(source.aliases))) print('%-10s %s' % ('', source.homepage)) if source.license: print('%-10s license: %s' % ('', source.license)) @@ -293,6 +343,7 @@ def handleSources(args): def handleSchemaShow(args): """Print a stored snapshot.""" + args.source = canonicalName(args.source) schema = loadSchema(args.source) if args.collection is None: @@ -324,6 +375,7 @@ def handleSchemaShow(args): def handleSchemaRefresh(args): """Re-harvest a snapshot and its catalogs.""" + args.source = canonicalName(args.source) client = makeClient(args) source = getSource(args.source, client) @@ -391,11 +443,88 @@ def handleSearch(args): return 0 +def handleReference(args): + """Validate a reference folder and, unless --check, build its IgBLAST base.""" + if not args.folder.is_dir(): + raise SourcererError('no such reference folder: %s' % args.folder) + + plan = Reference.planReference(args.folder, species=args.species) + print(plan.summary()) + + if not plan.ok: + raise SourcererError('no databases can be built from %s; check the file ' + 'names against _.fasta' % args.folder) + + if args.check: + return 0 + + if args.out is None: + raise SourcererError('--out is required to build; pass --check to only ' + 'validate the folder') + + Reference.buildFromPlan(plan, args.out, makeClient(args)) + log.info('wrote %s', args.out) + + return 0 + + +def handleReferenceDownload(args, source): + """Download germline sets and build an airrflow reference_base.""" + query = source.validateQuery(args.collection, collectFilters(args)) + if args.limit is not None: + query = type(query)(collection=query.collection, filters=query.filters, + limit=args.limit) + + units = source.searchUnits(query) + log.info('%d germline files for %s %s', + len(units), args.source, args.collection) + + if args.dry_run: + for unit in units: + print('%-48s %s' % (unit.unit_id, unit.url)) + log.info('dry run: nothing downloaded') + return 0 + + outdir = Path(args.outdir) + raw_dir = outdir / 'raw' + + entries, provenance = [], [] + for unit in units: + result = source.fetchUnit(unit, raw_dir, resume=not args.no_resume) + entries.append((unit, result.path)) + provenance.append(Provenance.buildUnitRecord(unit, result, outdir, {})) + + reference_dir = outdir / 'reference_base' + source.buildReference(entries, reference_dir).logSummary() + log.info('wrote %s', reference_dir) + + formats = ['reference'] + if args.igblast: + igblast_out = args.igblast_out or (outdir / 'igblast_base') + Reference.buildIgblastBase(reference_dir, igblast_out, source.client, + species=[args.collection]) + log.info('wrote %s', igblast_out) + formats.append('igblast') + + record = Provenance.writeDownloadMetadata( + outdir, args.source, args.collection, collectFilters(args), args.limit, + formats, provenance, schema=source.schema, license=source.license, + citation=source.citation) + log.info('wrote %s', record) + + return 0 + + def handleDownload(args): """Download and optionally convert matching data units.""" client = makeClient(args) source = getSource(args.source, client) + # Reference sources (germline sets) build a reference_base instead of + # converting repertoires to AIRR and writing a samplesheet. + if source.output == 'reference': + return handleReferenceDownload(args, source) + query = source.validateQuery(args.collection, collectFilters(args)) if args.limit is not None: query = type(query)(collection=query.collection, filters=query.filters, @@ -499,10 +628,16 @@ def main(): return handleSchemaRefresh(args) parser.parse_args([args.command, '--help']) - if args.command in REGISTRY: + if args.command == 'reference': + return handleReference(args) + + source_name = canonicalName(args.command) if args.command else None + if source_name in REGISTRY: # The action and collection levels are required subparsers, so - # argparse has already rejected a commandline missing either. - args.source = args.command + # argparse has already rejected a commandline missing either. The + # command may be an alias (e.g. 'airrc'); resolve it to the canonical + # source so schema and provenance use one name. + args.source = source_name if args.action == 'search': return handleSearch(args) if args.action == 'download': diff --git a/src/sourcerer/Exceptions.py b/src/sourcerer/Exceptions.py index 5920ba2..2f36382 100644 --- a/src/sourcerer/Exceptions.py +++ b/src/sourcerer/Exceptions.py @@ -40,6 +40,16 @@ class OasParseError(ParseError): pass +class ImgtParseError(ParseError): + """IMGT content did not match the expected structure.""" + pass + + +class OgrdbParseError(ParseError): + """OGRDB content did not match the expected structure.""" + pass + + class SchemaError(SourcererError): """A stored schema snapshot is missing, malformed or too new to understand.""" pass diff --git a/src/sourcerer/Reference.py b/src/sourcerer/Reference.py new file mode 100644 index 0000000..cc2a079 --- /dev/null +++ b/src/sourcerer/Reference.py @@ -0,0 +1,620 @@ +""" +Germline reference output + +The airrflow artifact for a germline source is not a samplesheet: it is a +reference tree that IgBLAST and Change-O read. This module is to the reference +sources what Airrflow.py is to the dataset sources -- the one place that knows +the shape of the output nf-core/airrflow expects, kept out of the sources +themselves so a second reference source inherits it unchanged. + +This module holds the builders that turn downloaded germline FASTAs into the two +directory layouts airrflow consumes: a reference_base of per-chain FASTAs, and -- +only when asked, because it needs the BLAST+ binary -- an igblast_base of BLAST +databases plus the internal_data and optional_file trees mirrored from NCBI. The +base class germline sources extend, ReferenceSource, lives in +sourcerer.Sources.Germline, kept there rather than here so this module never +imports from sourcerer.Sources and the two stay free of an import cycle. + +The reference_base keeps the source's own FASTA verbatim, gaps and all, exactly +as airrflow's bin/fetch_references.sh leaves it. Cleaning (gap removal, dedup) +happens only when the BLAST database is built, so nothing that reads the +reference for its IMGT numbering, such as Change-O's germline reconstruction, +loses it. +""" + +# Info +__author__ = 'Ayelet Peres' + +# Imports +import logging +import shutil +import subprocess +import tarfile +from dataclasses import dataclass, field +from pathlib import Path +from urllib.parse import urljoin, urlparse + +from bs4 import BeautifulSoup + +# Sourcerer imports +from sourcerer.Exceptions import SourcererError + +log = logging.getLogger(__name__) + +#: Species airrflow builds references for, and the leading directory in the +#: reference tree. New species are added here and in each source's SETS/CHAINS. +SPECIES = ('human', 'mouse') + +#: Receptor classes, matching airrflow's canonical database basenames. +LOCI = ('ig', 'tr') + +#: Gene segments, in the order IgBLAST names its databases. +SEGMENTS = ('v', 'd', 'j', 'c') + +#: The chains that make up each canonical (locus, segment) database. This is the +#: aggregation airrflow's ref2igblast.sh performs: one BLAST database per class +#: and segment, built from every locus in that class. +LOCUS_CHAINS = { + ('ig', 'v'): ('IGHV', 'IGKV', 'IGLV'), + ('ig', 'd'): ('IGHD',), + ('ig', 'j'): ('IGHJ', 'IGKJ', 'IGLJ'), + ('ig', 'c'): ('IGHC', 'IGKC', 'IGLC'), + ('tr', 'v'): ('TRAV', 'TRBV', 'TRDV', 'TRGV'), + ('tr', 'd'): ('TRBD', 'TRDD'), + ('tr', 'j'): ('TRAJ', 'TRBJ', 'TRDJ', 'TRGJ'), + ('tr', 'c'): ('TRAC', 'TRBC', 'TRDC', 'TRGC'), +} + +#: Every chain a reference FASTA may be named for, flattened from LOCUS_CHAINS. +KNOWN_CHAINS = frozenset(chain for chains in LOCUS_CHAINS.values() + for chain in chains) + +#: Which reference_base subdirectory a chain's FASTA lives in. Constant regions +#: are kept apart from V/D/J because airrflow's tree does, and amino acid V has +#: its own directory because it becomes a protein database rather than a +#: nucleotide one. +KIND_VDJ = 'vdj' +KIND_CONSTANT = 'constant' +KIND_AA = 'vdj_aa' + +#: NCBI's IgBLAST release trees, mirrored into igblast_base so that igblastn has +#: the auxiliary data it cannot derive from the germline FASTAs alone. The +#: old_* directories are the layout airrflow's fetch_igblastdb.sh already tracks. +NCBI_IGBLAST_ROOT = ('https://ftp.ncbi.nlm.nih.gov/blast/executables/igblast/' + 'release/') +NCBI_DATABASE_URL = urljoin(NCBI_IGBLAST_ROOT, 'database/') +NCBI_INTERNAL_URL = urljoin(NCBI_IGBLAST_ROOT, 'old_internal_data/') +NCBI_OPTIONAL_URL = urljoin(NCBI_IGBLAST_ROOT, 'old_optional_file/') + +#: Archives NCBI ships inside database/ that have to be unpacked in place for the +#: mirrored tree to be usable. +NCBI_TAR_ARCHIVES = ('mouse_gl_VDJ.tar', 'rhesus_monkey_VJ.tar') + + +@dataclass +class ReferenceReport: + """ + Summary of a reference build, in the same spirit as OAS's conversion report. + + Arguments: + written (list): reference_base FASTAs written, as (chain, path) or basename. + built (list): canonical BLAST database basenames created. + skipped_empty (list): canonical databases skipped because no chain in them + had any sequence. This is the normal outcome for what a source does not + cover, such as TR from OGRDB. + """ + written: list = field(default_factory=list) + built: list = field(default_factory=list) + skipped_empty: list = field(default_factory=list) + + def logSummary(self): + """Log a one-line summary of what the build produced.""" + log.info('reference: %d files written, %d databases built, %d skipped', + len(self.written), len(self.built), len(self.skipped_empty)) + + +@dataclass +class ReferencePlan: + """ + What building an IgBLAST base from a reference folder would produce. + + Computed without running makeblastdb, so it doubles as the format check: it + says which databases would build, which come up empty, which files were not + recognised, and where duplicate allele names were dropped. + + Arguments: + found_species (list): species seen in the folder. + databases (list): (basename, dbtype, records) that will build. + empty (list): canonical basenames with no sequence to build from. + unrecognized (list): FASTA paths whose names are not in the reference format. + empty_files (list): recognised FASTA paths that held no sequence. + duplicates (dict): basename to the number of duplicate names dropped. + """ + found_species: list = field(default_factory=list) + databases: list = field(default_factory=list) + empty: list = field(default_factory=list) + unrecognized: list = field(default_factory=list) + empty_files: list = field(default_factory=list) + duplicates: dict = field(default_factory=dict) + + @property + def ok(self): + """bool: True if at least one database can be built.""" + return bool(self.databases) + + def summary(self): + """ + Render the plan as a human-readable report. + + Returns: + str: the report. + """ + lines = ['species found: %s' % (', '.join(self.found_species) or 'none')] + if self.databases: + lines.append('databases to build (%d):' % len(self.databases)) + for basename, dbtype, records in self.databases: + dropped = self.duplicates.get(basename) + note = ' (%d duplicate name(s) dropped)' % dropped if dropped else '' + lines.append(' %-16s %5d seq %s%s' + % (basename, len(records), dbtype, note)) + if self.empty: + lines.append('empty, nothing to build: %s' % ', '.join(self.empty)) + for path in self.empty_files: + lines.append('warning: %s held no sequence' % path) + for path in self.unrecognized: + lines.append('warning: %s is not in the reference naming format, ' + 'skipped' % path.name) + + return '\n'.join(lines) + + +# --------------------------------------------------------------------------- +# FASTA helpers +# --------------------------------------------------------------------------- + +def parseFasta(text): + """ + Read FASTA text into (header, sequence) pairs, in source order. + + Whitespace inside a sequence is collapsed; the header keeps everything after + the '>' verbatim, because a source's own header, IMGT's pipe-delimited line + for one and OGRDB's allele name for another, is what identifies the allele. + + Arguments: + text (str): FASTA text. + + Returns: + list: (header, sequence) tuples. + """ + records = [] + header, seq = None, [] + for line in text.splitlines(): + line = line.strip() + if not line: + continue + if line.startswith('>'): + if header is not None: + records.append((header, ''.join(seq))) + header = line[1:] + seq = [] + else: + seq.append(''.join(line.split())) + if header is not None: + records.append((header, ''.join(seq))) + + return records + + +def alleleName(header): + """ + Take the allele name from a FASTA header. + + IMGT headers are pipe-delimited and put the name in the second field + (``>X02897|IGHV1-2*02|Homo sapiens|F|...``); OGRDB writes the bare name. The + first token is used when there is no pipe, so both are handled by one rule. + + Arguments: + header (str): the header line without its leading '>'. + + Returns: + str: the allele name. + """ + if '|' in header: + return header.split('|')[1].strip() + + return header.split()[0] + + +def writeFastaText(path, records): + """ + Write (header, sequence) pairs to a FASTA file verbatim, one line per part. + + Arguments: + path (Path): output path. Parent directories are created. + records (iterable): (header, sequence) tuples. + + Returns: + int: the number of records written. + """ + path = Path(path) + path.parent.mkdir(parents=True, exist_ok=True) + + written = 0 + with open(path, 'w') as handle: + for header, sequence in records: + handle.write('>%s\n%s\n' % (header, sequence)) + written += 1 + + return written + + +def cleanForBlast(records): + """ + Prepare germline records for makeblastdb. + + Gaps are removed, sequences upper-cased and duplicate names dropped, keeping + the first. Deduplication is not cosmetic: makeblastdb -parse_seqids refuses a + database with a repeated identifier, so a duplicate allele name is a hard + failure rather than a warning. The name is taken with alleleName so an IMGT + pipe header collapses to just the allele, which is what IgBLAST reports. + + Arguments: + records (iterable): (header, sequence) tuples. + + Returns: + list: (name, sequence) tuples, cleaned and de-duplicated. + """ + seen = set() + cleaned = [] + for header, sequence in records: + name = alleleName(header) + if name in seen: + continue + seen.add(name) + cleaned.append((name, sequence.replace('.', '').upper())) + + return cleaned + + +# --------------------------------------------------------------------------- +# reference_base +# --------------------------------------------------------------------------- + +def referenceFastaPath(reference_dir, prefix, species, kind, chain): + """ + Locate one chain's FASTA in the reference tree. + + The layout matches airrflow's: ``//__``, + with amino acid V spelled ``_aa__`` so a nucleotide + and a protein file for the same chain do not collide. + + Arguments: + reference_dir (Path): the reference_base root. + prefix (str): the source tag, 'imgt' or 'airrc'. + species (str): the species. + kind (str): the subdirectory, one of KIND_VDJ, KIND_CONSTANT, KIND_AA. + chain (str): the chain, e.g. 'IGHV'. + + Returns: + Path: where the chain's FASTA belongs. + """ + if kind == KIND_AA: + name = '%s_aa_%s_%s.fasta' % (prefix, species, chain) + else: + name = '%s_%s_%s.fasta' % (prefix, species, chain) + + return Path(reference_dir) / species / kind / name + + +def parseReferenceName(filename): + """ + Read (species, chain, is_aa) from a reference FASTA's name. + + The accepted form is ``[_][aa_]_.fasta``: an optional + source prefix (``imgt_``, ``airrc_``, ...) that is ignored, an optional + ``aa_`` marking translated V, the species, and a known chain. Only the name is + read, never the directory, so a file nested in a reference_base and a file in a + flat folder are recognised the same way -- which is what lets both layouts + build. + + Arguments: + filename (str): a FASTA file's basename. + + Returns: + tuple: (species, chain, is_aa), or None if the name does not match. + """ + if not filename.endswith('.fasta'): + return None + + tokens = filename[:-len('.fasta')].split('_') + for index, token in enumerate(tokens): + if token in SPECIES and index + 1 < len(tokens): + chain = tokens[index + 1] + if chain in KNOWN_CHAINS: + return token, chain, 'aa' in tokens[:index] + + return None + + +def discoverReference(reference_dir): + """ + Find every reference FASTA under a folder, by filename, in any layout. + + The folder is searched recursively, so a nested reference_base and a flat + folder of FASTAs are both handled; classification is by name alone. + + Arguments: + reference_dir (Path): a reference_base tree or a flat folder of FASTAs. + + Returns: + tuple: (files, unrecognized) where files is a list of + (species, chain, is_aa, Path), and unrecognized is the list of .fasta paths + whose names are not in the reference format. + """ + files, unrecognized = [], [] + for path in sorted(Path(reference_dir).rglob('*.fasta')): + parsed = parseReferenceName(path.name) + if parsed is None: + unrecognized.append(path) + else: + species, chain, is_aa = parsed + files.append((species, chain, is_aa, path)) + + return files, unrecognized + + +# --------------------------------------------------------------------------- +# igblast_base +# --------------------------------------------------------------------------- + +def runMakeblastdb(fasta, out_base, dbtype): + """ + Build one BLAST database from a cleaned FASTA. + + Arguments: + fasta (Path): the input FASTA. + out_base (Path): the database basename, without an extension. + dbtype (str): 'nucl' or 'prot'. + + Raises: + SourcererError: if makeblastdb is not on PATH or exits non-zero. + """ + if shutil.which('makeblastdb') is None: + raise SourcererError( + 'makeblastdb not found on PATH; install NCBI BLAST+ (for example ' + 'conda install -c bioconda blast) or drop --igblast to write only ' + 'the reference FASTAs') + + result = subprocess.run( + ['makeblastdb', '-parse_seqids', '-dbtype', dbtype, + '-in', str(fasta), '-out', str(out_base)], + capture_output=True, text=True) + if result.returncode != 0: + raise SourcererError('makeblastdb failed for %s: %s' + % (Path(fasta).name, result.stderr.strip())) + + +def planReference(reference_dir, species=None): + """ + Work out which IgBLAST databases a reference folder would produce. + + Files are grouped by name into the canonical (species, locus, segment) + databases airrflow expects, whatever layout they came in and whatever prefix + wrote them, so a nested reference_base and a flat folder both plan. No + makeblastdb is run, so this is also the format check: the returned plan says + what would build, what is empty, what was not recognised, and where duplicate + names were dropped. + + Arguments: + reference_dir (Path): a reference_base tree or a flat folder of FASTAs. + species (iterable): limit to these species, or None for every species found. + + Returns: + ReferencePlan: the databases that would build and the diagnostics. + """ + files, unrecognized = discoverReference(reference_dir) + found = sorted({item[0] for item in files}) + wanted = list(species) if species else found + + contents = {path: parseFasta(path.read_text()) + for _sp, _chain, _aa, path in files} + empty_files = [path for path, records in contents.items() if not records] + + def collect(sp, chains, is_aa): + records = [] + for f_sp, f_chain, f_aa, f_path in files: + if f_sp == sp and f_aa == is_aa and f_chain in chains: + records.extend(contents[f_path]) + return records + + plan = ReferencePlan(found_species=found, unrecognized=unrecognized, + empty_files=empty_files) + for sp in wanted: + for locus in LOCI: + for segment in SEGMENTS: + _addToPlan(plan, '%s_%s_%s' % (sp, locus, segment), 'nucl', + collect(sp, LOCUS_CHAINS[(locus, segment)], False)) + # Amino acid V is a protein database, built only when the folder + # actually carries translated V (OGRDB, for one, does not), so an + # absent one is not reported as an empty gap. + _addToPlan(plan, 'aa_%s_%s_v' % (sp, locus), 'prot', + collect(sp, LOCUS_CHAINS[(locus, 'v')], True), + keep_empty=False) + + return plan + + +def _addToPlan(plan, basename, dbtype, records, keep_empty=True): + """ + Clean one canonical database's records and record it on the plan. + + Arguments: + plan (ReferencePlan): the plan to add to. + basename (str): the canonical database basename. + dbtype (str): 'nucl' or 'prot'. + records (list): the raw (header, sequence) tuples gathered for it. + keep_empty (bool): whether an empty database is worth reporting as a gap. + """ + if not records: + if keep_empty: + plan.empty.append(basename) + return + + cleaned = cleanForBlast(records) + dropped = len(records) - len(cleaned) + if dropped: + plan.duplicates[basename] = dropped + plan.databases.append((basename, dbtype, cleaned)) + + +def buildIgblastBase(reference_dir, out_dir, client, species=None): + """ + Build the IgBLAST database tree airrflow expects from a reference folder. + + A thin wrapper over planReference and buildFromPlan, kept so the download + path and the standalone `reference build` command share one code path. + + Arguments: + reference_dir (Path): a reference_base tree or a flat folder of FASTAs. + out_dir (Path): the igblast_base to write. + client (HttpClient): used to mirror the NCBI support trees. + species (iterable): limit to these species, or None for every species found. + + Returns: + ReferenceReport: what was built and what was skipped. + + Raises: + SourcererError: if makeblastdb is unavailable. + """ + plan = planReference(reference_dir, species=species) + if plan.unrecognized: + log.warning('%d file(s) skipped: names not in the reference format', + len(plan.unrecognized)) + + return buildFromPlan(plan, out_dir, client) + + +def buildFromPlan(plan, out_dir, client): + """ + Write the databases a plan describes, then mirror the NCBI support trees. + + Arguments: + plan (ReferencePlan): the databases to build, from planReference. + out_dir (Path): the igblast_base to write; fasta/ and database/ are created + inside it, alongside the mirrored internal_data/ and optional_file/. + client (HttpClient): used to mirror the NCBI support trees. + + Returns: + ReferenceReport: what was built and what was skipped. + + Raises: + SourcererError: if makeblastdb is unavailable. + """ + out_dir = Path(out_dir) + fasta_out = out_dir / 'fasta' + db_out = out_dir / 'database' + fasta_out.mkdir(parents=True, exist_ok=True) + db_out.mkdir(parents=True, exist_ok=True) + + report = ReferenceReport(skipped_empty=list(plan.empty)) + for basename, dbtype, records in plan.databases: + fasta = fasta_out / ('%s.fasta' % basename) + with open(fasta, 'w') as handle: + for name, sequence in records: + handle.write('>%s\n%s\n' % (name, sequence)) + runMakeblastdb(fasta, db_out / basename, dbtype) + report.built.append(basename) + + mirrorSupport(out_dir, client) + report.logSummary() + + return report + + +def mirrorSupport(out_dir, client): + """ + Mirror the NCBI IgBLAST support trees into an igblast_base. + + database/, internal_data/ and optional_file/ are copied from NCBI's release + directory, and the tar archives NCBI ships inside database/ are unpacked in + place, so the result matches what airrflow's fetch_igblastdb.sh produces. + + Arguments: + out_dir (Path): the igblast_base root. + client (HttpClient): the shared HTTP client. + """ + out_dir = Path(out_dir) + database_dir = out_dir / 'database' + mirrorTree(NCBI_DATABASE_URL, database_dir, client) + for name in NCBI_TAR_ARCHIVES: + archive = database_dir / name + if archive.exists(): + extractTar(archive, database_dir) + + mirrorTree(NCBI_INTERNAL_URL, out_dir / 'internal_data', client) + mirrorTree(NCBI_OPTIONAL_URL, out_dir / 'optional_file', client) + + +def mirrorTree(url, dest_dir, client, seen=None): + """ + Recursively mirror an Apache/NCBI directory index into a local tree. + + Only links that stay under the starting URL are followed, so a parent link + or an absolute link elsewhere on the host cannot walk the mirror out of the + subtree it was pointed at. + + Arguments: + url (str): the directory index URL, ending in '/'. + dest_dir (Path): where to mirror it. + client (HttpClient): the shared HTTP client. + seen (set): URLs already visited, to guard against a self-referential index. + """ + seen = seen if seen is not None else set() + if url in seen: + return + seen.add(url) + + dest_dir = Path(dest_dir) + dest_dir.mkdir(parents=True, exist_ok=True) + + soup = BeautifulSoup(client.get(url).text, 'html.parser') + for anchor in soup.find_all('a'): + href = anchor.get('href') + if not href or href in ('../', './', '/'): + continue + + child = urljoin(url, href) + if urlparse(child).scheme not in ('http', 'https'): + continue + if not child.startswith(url): + continue + + name = Path(urlparse(child).path).name + if not name: + continue + + if href.endswith('/'): + mirrorTree(child, dest_dir / name, client, seen) + else: + client.fetch(child, dest_dir / name, progress=False) + + +def extractTar(archive, dest_dir): + """ + Unpack a tar archive, refusing any member that would escape the destination. + + Arguments: + archive (Path): the tar file. + dest_dir (Path): where to extract it. + + Raises: + SourcererError: if a member path points outside dest_dir. + """ + dest_dir = Path(dest_dir).resolve() + with tarfile.open(archive) as tar: + for member in tar.getmembers(): + target = (dest_dir / member.name).resolve() + if target != dest_dir and dest_dir not in target.parents: + raise SourcererError('unsafe path %s in %s' + % (member.name, Path(archive).name)) + tar.extractall(dest_dir) diff --git a/src/sourcerer/Sources/AirrcImgt.py b/src/sourcerer/Sources/AirrcImgt.py new file mode 100644 index 0000000..782db23 --- /dev/null +++ b/src/sourcerer/Sources/AirrcImgt.py @@ -0,0 +1,172 @@ +""" +AIRR-C germline sets blended with IMGT + +A germline reference that takes each locus from the source that curates it best: +the immunoglobulin V, D and J from OGRDB's AIRR-C sets, and everything OGRDB does +not cover -- all of the T-cell receptor, and the immunoglobulin constants that +are not in a published set -- from IMGT. It is the reference nf-core/airrflow +builds for its ``airrc-imgt`` database type. + +Rather than reimplement either source, this composes them: it asks the OGRDB +source for the immunoglobulin sets and the IMGT source for the gap, tags each +download with which one it came from, and lets each build its own files back into +one reference tree, so an OGRDB allele lands as ``airrc_...`` and an IMGT allele +as ``imgt_...`` exactly as they would from the sources alone. Amino acid V is not +included, matching what airrflow's airrc-imgt build uses. +""" + +# Info +__author__ = 'Ayelet Peres' + +# Imports +import logging +from datetime import UTC + +# Sourcerer imports +from sourcerer.Reference import KIND_AA, ReferenceReport +from sourcerer.Sources.Base import Query +from sourcerer.Sources.Germline import ReferenceSource +from sourcerer.Sources.Imgt import ImgtSource +from sourcerer.Sources.Ogrdb import OgrdbSource + +log = logging.getLogger(__name__) + +#: The IMGT immunoglobulin constants to take per species: the ones OGRDB has no +#: set for. Human IGHC comes from OGRDB, so only the light constants are taken; +#: mouse has no constant set at all, so all three heavy and light come from IMGT. +IMGT_IG_CONSTANTS = {'human': ('IGKC', 'IGLC'), + 'mouse': ('IGHC', 'IGKC', 'IGLC')} + +#: Which download a unit came from, recorded so buildReference can hand each unit +#: back to the source that knows how to read it. +VIA = 'via' + + +class AirrcImgtSource(ReferenceSource): + """ + OGRDB immunoglobulin sets blended with IMGT for TR and the IG constants. + """ + + name = 'airrc-imgt' + description = ('AIRR-C immunoglobulin sets blended with IMGT for TR and the ' + 'remaining constants') + homepage = 'https://ogrdb.airr-community.org/' + collections = ('human', 'mouse') + collection_help = {'human': 'Homo sapiens blended reference', + 'mouse': 'Mus musculus blended reference'} + license = ('OGRDB data under CC BY 4.0 and IMGT data under the IMGT terms of ' + 'use; cite both OGRDB and IMGT') + citation = OgrdbSource.citation + ImgtSource.citation + + def __init__(self, client, schema=None): + """ + Arguments: + client (HttpClient): the shared HTTP client, passed to both sources. + schema (SourceSchema): the loaded snapshot, or None to load on demand. + """ + super().__init__(client, schema) + self._ogrdb = OgrdbSource(client) + self._imgt = ImgtSource(client) + + def harvestSchema(self): + """ + Build a snapshot for the blended source. + + The blend takes a whole species at a time and has no filters of its own, + so the snapshot is just its collections; the drift checks that matter live + on the imgt and ogrdb snapshots this composes. + + Returns: + SourceSchema: the snapshot. + """ + from datetime import datetime + + from sourcerer.Schema import Collection, SourceSchema + from sourcerer.Version import __version__ + + return SourceSchema( + source=self.name, + harvested=datetime.now(UTC).strftime('%Y-%m-%dT%H:%M:%SZ'), + harvested_by='sourcerer %s' % __version__, + source_urls={'ogrdb': self._ogrdb.name, 'imgt': self._imgt.name}, + collections={sp: Collection(name=sp) for sp in self.collections}) + + def searchUnits(self, query): + """ + Resolve a query to the OGRDB and IMGT files the blend needs. + + Arguments: + query (Query): the validated request; collection is the species. + + Returns: + list: DataUnit objects, each tagged with which source produced it. + """ + species = query.collection + + units = [] + for unit in self._ogrdb.searchUnits( + Query(collection=species, filters={'locus': '*'})): + unit.metadata[VIA] = self._ogrdb.name + units.append(unit) + + for unit in self._imgtGapUnits(species): + unit.metadata[VIA] = self._imgt.name + units.append(unit) + + if query.limit is not None: + units = units[:query.limit] + + return units + + def _imgtGapUnits(self, species): + """ + Pick the IMGT files that fill what OGRDB does not cover. + + That is every T-cell receptor chain, and the immunoglobulin constants + without an OGRDB set; the immunoglobulin V, D and J come from OGRDB, and + amino acid V is not part of this blend. + + Arguments: + species (str): the species. + + Returns: + list: DataUnit objects from the IMGT source. + """ + constants = IMGT_IG_CONSTANTS.get(species, ()) + gap = [] + for unit in self._imgt.searchUnits( + Query(collection=species, + filters={'locus': '*', 'segment': '*'})): + meta = unit.metadata + if meta['kind'] == KIND_AA: + continue + if meta['locus'].startswith('TR'): + gap.append(unit) + elif meta['kind'] == 'constant' and meta['chain'] in constants: + gap.append(unit) + + return gap + + def buildReference(self, entries, reference_dir): + """ + Let each source build its own files back into one reference tree. + + Arguments: + entries (list): (DataUnit, Path) pairs from the fetch step. + reference_dir (Path): the reference_base root. + + Returns: + ReferenceReport: the files written by both sources. + """ + ogrdb_entries = [pair for pair in entries + if pair[0].metadata.get(VIA) == self._ogrdb.name] + imgt_entries = [pair for pair in entries + if pair[0].metadata.get(VIA) == self._imgt.name] + + report = ReferenceReport() + report.written.extend( + self._ogrdb.buildReference(ogrdb_entries, reference_dir).written) + report.written.extend( + self._imgt.buildReference(imgt_entries, reference_dir).written) + + return report diff --git a/src/sourcerer/Sources/Base.py b/src/sourcerer/Sources/Base.py index adba6f6..c149d62 100644 --- a/src/sourcerer/Sources/Base.py +++ b/src/sourcerer/Sources/Base.py @@ -101,6 +101,10 @@ class SourceBase(ABC): #: Short name used on the commandline and as the schema directory name. name = None + #: Alternative commandline names for the same source, e.g. ('airrc',) for + #: OGRDB. They share the source's subcommand, flags and schema; the canonical + #: ``name`` is what schema and provenance are keyed on. + aliases = () #: One line description for `sourcerer sources list`. description = '' #: Where a human can read about the source. @@ -118,6 +122,12 @@ class SourceBase(ABC): #: the record of what was downloaded travels with a reminder of how to #: give the source credit for it. citation = () + #: What the source produces, and therefore which output path `download` + #: drives. 'dataset' sources are repertoires: they convert to AIRR/FASTA and + #: write an airrflow samplesheet. 'reference' sources are germline sets: they + #: build an airrflow germline reference_base instead, and never touch the + #: rearrangement conversion path. See sourcerer.Reference.ReferenceSource. + output = 'dataset' def __init__(self, client, schema=None): """ diff --git a/src/sourcerer/Sources/Germline.py b/src/sourcerer/Sources/Germline.py new file mode 100644 index 0000000..4d06b3f --- /dev/null +++ b/src/sourcerer/Sources/Germline.py @@ -0,0 +1,76 @@ +""" +Germline reference sources + +A germline reference is a set of allele sequences. IMGT and OGRDB provide them, +and ReferenceSource turns a download into the airrflow reference tree through +buildReference, reusing SourceBase's download machinery: searchUnits, fetchUnit, +the schema and the HTTP client. The output shaping itself lives in +sourcerer.Reference. + +SourceBase also declares readUnit and normalizeChunk, the conversion step for +sources that emit AIRR records; a germline source has nothing to convert, so +those two are filled in here with a clear error rather than left for each source +to repeat. +""" + +# Info +__author__ = 'Ayelet Peres' + +# Sourcerer imports +from sourcerer.Exceptions import SourcererError +from sourcerer.Reference import referenceFastaPath, writeFastaText +from sourcerer.Sources.Base import SourceBase + + +class ReferenceSource(SourceBase): + """ + Base class for germline reference sources such as IMGT and OGRDB. + """ + + output = 'reference' + + #: The tag written into reference_base filenames, e.g. 'imgt' gives + #: imgt_human_IGHV.fasta. Subclasses set this. + prefix = None + + def readUnit(self, path, unit): + """Unused: a germline source has nothing to convert.""" + raise SourcererError('%s builds a germline reference; there is nothing ' + 'to convert' % self.name) + + def normalizeChunk(self, metadata, chunk, unit, offset, report): + """Unused: a germline source has nothing to convert.""" + raise SourcererError('%s builds a germline reference; there is nothing ' + 'to convert' % self.name) + + def buildReference(self, entries, reference_dir): + """ + Turn downloaded germline files into an airrflow reference_base. + + Arguments: + entries (list): (DataUnit, Path) pairs, one per fetched file. + reference_dir (Path): the reference_base root to write into. + + Returns: + ReferenceReport: the files written. + """ + raise NotImplementedError + + def writeChain(self, reference_dir, species, kind, chain, records): + """ + Write one chain's FASTA into the reference tree. + + Arguments: + reference_dir (Path): the reference_base root. + species (str): the species. + kind (str): the subdirectory (KIND_VDJ, KIND_CONSTANT or KIND_AA). + chain (str): the chain, e.g. 'IGHV'. + records (iterable): (header, sequence) tuples to write verbatim. + + Returns: + Path: the file written. + """ + path = referenceFastaPath(reference_dir, self.prefix, species, kind, chain) + writeFastaText(path, records) + + return path diff --git a/src/sourcerer/Sources/Imgt.py b/src/sourcerer/Sources/Imgt.py new file mode 100644 index 0000000..d19a259 --- /dev/null +++ b/src/sourcerer/Sources/Imgt.py @@ -0,0 +1,306 @@ +""" +IMGT/GENE-DB + +IMGT exposes no data API. Its GENE-DB GENElect page takes a numbered query and a +chain and returns an HTML page with the FASTA embedded in its second ``
``
+block; the first holds the query echo. A query can fail and still return HTTP
+200, so a valid answer is not the status alone but the presence of that second
+block with sequence in it, which is what isValidResponse checks and the weekly
+API canary relies on.
+
+One germline file per (species, chain) is fetched, matching how airrflow's
+bin/fetch_references.sh drives GENElect: query 7.14 for V, D and J nucleotide,
+14.1 for constant (7.5 for the mouse kappa and lambda constants, which 14.1 does
+not serve), and 7.3 for the translated V. The files are written into the
+reference_base with their IMGT headers and gaps intact; see sourcerer.Reference
+for why cleaning is deferred to the database build.
+"""
+
+# Info
+__author__ = 'Ayelet Peres'
+
+# Imports
+import logging
+from datetime import UTC
+from urllib.parse import quote
+
+from bs4 import BeautifulSoup
+
+# Sourcerer imports
+from sourcerer.Exceptions import ImgtParseError
+from sourcerer.Reference import (
+    KIND_AA,
+    KIND_CONSTANT,
+    KIND_VDJ,
+    ReferenceReport,
+    parseFasta,
+)
+from sourcerer.Sources.Base import DataUnit
+from sourcerer.Sources.Germline import ReferenceSource
+
+log = logging.getLogger(__name__)
+
+#: Endpoints.
+GENELECT = 'https://www.imgt.org/genedb/GENElect'
+RELEASE_URL = 'https://www.imgt.org/download/GENE-DB/RELEASE'
+
+#: GENElect query numbers, per chain kind.
+Q_VDJ = '7.14'          # V, D and J nucleotide
+Q_CONSTANT = '14.1'     # constant nucleotide
+Q_CONSTANT_MOUSE = '7.5'  # mouse IGKC and IGLC, which 14.1 does not serve
+Q_AA = '7.3'            # translated V
+
+#: Species as GENElect wants them in the query string, and as they appear in the
+#: FASTA headers. The query form is pre-encoded so it is not double-escaped.
+SPECIES_QUERY = {'human': 'Homo%20sapiens', 'mouse': 'Mus%20musculus'}
+SPECIES_LABEL = {'human': 'Homo sapiens', 'mouse': 'Mus musculus'}
+
+#: Chains fetched as V/D/J nucleotide.
+VDJ_CHAINS = ('IGHV', 'IGHD', 'IGHJ', 'IGKV', 'IGKJ', 'IGLV', 'IGLJ',
+              'TRAV', 'TRAJ', 'TRBV', 'TRBD', 'TRBJ',
+              'TRDV', 'TRDD', 'TRDJ', 'TRGV', 'TRGJ')
+
+#: Chains fetched as constant nucleotide.
+CONSTANT_CHAINS = ('IGHC', 'IGKC', 'IGLC', 'TRAC', 'TRBC', 'TRGC', 'TRDC')
+
+#: Chains fetched as translated V.
+AA_CHAINS = ('IGHV', 'IGKV', 'IGLV', 'TRAV', 'TRBV', 'TRDV', 'TRGV')
+
+#: Loci and segments a search can be narrowed to, offered as filter flags.
+LOCI = ('IGH', 'IGK', 'IGL', 'TRA', 'TRB', 'TRG', 'TRD')
+SEGMENTS = ('V', 'D', 'J', 'C')
+
+
+def buildQueryUrl(species, query, chain, label=None):
+    """
+    Build a GENElect query URL.
+
+    Arguments:
+      species (str): the species key, e.g. 'human'.
+      query (str): the GENElect query number, e.g. '7.14'.
+      chain (str): the chain, e.g. 'IGHV'.
+      label (str): an optional IMGTlabel qualifier.
+
+    Returns:
+      str: the absolute query URL.
+    """
+    url = '%s?query=%s+%s&species=%s' % (GENELECT, quote(query), chain,
+                                         SPECIES_QUERY[species])
+    if label:
+        url += '&IMGTlabel=%s' % label
+
+    return url
+
+
+def isValidResponse(html):
+    """
+    Report whether a GENElect reply actually carries a germline FASTA.
+
+    GENElect answers a failed query with HTTP 200 and an error page, so a live
+    check cannot trust the status code. A real answer has a second ``
``
+    block, and that block has sequence in it. Both conditions are required.
+
+    Arguments:
+      html (str): the GENElect reply body.
+
+    Returns:
+      bool: True if the reply contains a non-empty germline FASTA.
+    """
+    blocks = BeautifulSoup(html, 'html.parser').find_all('pre')
+    if len(blocks) < 2:
+        return False
+
+    return '>' in blocks[1].get_text()
+
+
+def extractFasta(html, species):
+    """
+    Pull the germline FASTA out of a GENElect reply.
+
+    The FASTA is the second ``
`` block. The species name in the headers has
+    its spaces replaced with underscores, as airrflow does, so a header stays one
+    whitespace-delimited field.
+
+    Arguments:
+      html (str): the GENElect reply body.
+      species (str): the species key, for the header rewrite.
+
+    Returns:
+      str: the FASTA text.
+
+    Raises:
+      ImgtParseError: if the reply has no second ``
`` block, which means the
+        query failed or the page layout changed.
+    """
+    blocks = BeautifulSoup(html, 'html.parser').find_all('pre')
+    if len(blocks) < 2:
+        raise ImgtParseError(
+            'GENElect reply has fewer than two 
 blocks; the query failed or '
+            'the page layout changed')
+
+    text = blocks[1].get_text()
+    label = SPECIES_LABEL.get(species)
+    if label:
+        text = text.replace(label, label.replace(' ', '_'))
+
+    return text
+
+
+def _chainPlan(species):
+    """
+    Enumerate every (chain, kind, query, locus, segment) this source fetches.
+
+    Arguments:
+      species (str): the species key, used to route the mouse constant queries.
+
+    Returns:
+      list: dicts describing one germline file each.
+    """
+    plan = []
+    for chain in VDJ_CHAINS:
+        plan.append({'chain': chain, 'kind': KIND_VDJ, 'query': Q_VDJ,
+                     'locus': chain[:3], 'segment': chain[3]})
+    for chain in CONSTANT_CHAINS:
+        query = Q_CONSTANT
+        if species == 'mouse' and chain in ('IGKC', 'IGLC'):
+            query = Q_CONSTANT_MOUSE
+        plan.append({'chain': chain, 'kind': KIND_CONSTANT, 'query': query,
+                     'locus': chain[:3], 'segment': 'C'})
+    for chain in AA_CHAINS:
+        plan.append({'chain': chain, 'kind': KIND_AA, 'query': Q_AA,
+                     'locus': chain[:3], 'segment': 'V'})
+
+    return plan
+
+
+class ImgtSource(ReferenceSource):
+    """
+    The IMGT/GENE-DB germline reference source.
+    """
+
+    name = 'imgt'
+    prefix = 'imgt'
+    description = 'IMGT/GENE-DB: germline V, D, J and C reference sequences'
+    homepage = 'https://www.imgt.org/genedb/'
+    collections = ('human', 'mouse')
+    collection_help = {'human': 'Homo sapiens germline reference',
+                       'mouse': 'Mus musculus germline reference'}
+
+    #: IMGT's reuse terms are not an open-data licence; germline data may be used
+    #: for research on condition IMGT is cited. Recorded so a reader of a download
+    #: directory sees the obligation without having to consult IMGT separately.
+    license = ('IMGT terms of use (https://www.imgt.org/about/termsofuse.php); '
+               'cite IMGT, the international ImMunoGeneTics information system')
+    citation = (
+        'Lefranc MP, Giudicelli V, Duroux P, et al. IMGT, the international '
+        'ImMunoGeneTics information system 25 years on. Nucleic Acids Res. '
+        '2015;43(Database issue):D413-D422. doi:10.1093/nar/gku1056',
+    )
+
+    def harvestSchema(self):
+        """
+        Contact IMGT for its current release and build a fresh snapshot.
+
+        GENElect has no field-listing endpoint, so the searchable vocabulary is
+        the fixed set of loci and segments this source knows how to query. The
+        network call to the release file is what turns a refresh into a genuine
+        liveness check rather than a rewrite of a constant.
+
+        Returns:
+          SourceSchema: the harvested snapshot.
+        """
+        from datetime import datetime
+
+        from sourcerer.Schema import Collection, Field, SourceSchema
+        from sourcerer.Version import __version__
+
+        # A liveness check, not stored: the release tag changes with every IMGT
+        # build, and keeping it in the snapshot would make a monthly refresh
+        # rewrite a tracked file with no change in the vocabulary.
+        log.info('IMGT release: %s', self.fetchRelease() or 'unknown')
+
+        fields = (Field(name='locus', values=LOCI),
+                  Field(name='segment', values=SEGMENTS))
+        collections = {sp: Collection(name=sp, fields=fields)
+                       for sp in self.collections}
+
+        return SourceSchema(
+            source=self.name,
+            harvested=datetime.now(UTC).strftime('%Y-%m-%dT%H:%M:%SZ'),
+            harvested_by='sourcerer %s' % __version__,
+            source_urls={'genelect': GENELECT, 'release': RELEASE_URL},
+            parse_contracts={'fasta_block': 'second 
 element',
+                             'header': 'pipe-delimited, allele in field 2'},
+            collections=collections)
+
+    def fetchRelease(self):
+        """
+        Return IMGT's current GENE-DB release tag.
+
+        Returns:
+          str: the release tag, e.g. '202619-7', or '' if it cannot be read.
+        """
+        try:
+            return self.client.get(RELEASE_URL).text.strip()
+        except Exception as error:
+            log.warning('could not read the IMGT release tag: %s', error)
+            return ''
+
+    def searchUnits(self, query):
+        """
+        Resolve a query to the germline files to fetch.
+
+        Arguments:
+          query (Query): the validated request; collection is the species, and
+            the locus and segment filters narrow which chains are fetched.
+
+        Returns:
+          list: DataUnit objects, one per germline file.
+        """
+        species = query.collection
+        locus = query.filters.get('locus', '*')
+        segment = query.filters.get('segment', '*')
+
+        units = []
+        for item in _chainPlan(species):
+            if locus not in ('*', item['locus']):
+                continue
+            if segment not in ('*', item['segment']):
+                continue
+
+            url = buildQueryUrl(species, item['query'], item['chain'])
+            unit_id = '%s/%s.html' % (item['kind'], item['chain'])
+            units.append(DataUnit(
+                unit_id=unit_id, collection=species, url=url,
+                metadata={'species': species, 'chain': item['chain'],
+                          'kind': item['kind'], 'locus': item['locus'],
+                          'segment': item['segment'], 'query': item['query']}))
+
+        if query.limit is not None:
+            units = units[:query.limit]
+
+        return units
+
+    def buildReference(self, entries, reference_dir):
+        """
+        Extract each downloaded GENElect page into the reference tree.
+
+        Arguments:
+          entries (list): (DataUnit, Path) pairs from the fetch step.
+          reference_dir (Path): the reference_base root.
+
+        Returns:
+          ReferenceReport: the files written.
+        """
+        report = ReferenceReport()
+        for unit, path in entries:
+            html = path.read_text()
+            fasta = extractFasta(html, unit.metadata['species'])
+            records = parseFasta(fasta)
+            written = self.writeChain(reference_dir, unit.metadata['species'],
+                                      unit.metadata['kind'],
+                                      unit.metadata['chain'], records)
+            report.written.append((unit.metadata['chain'], written))
+            log.info('%s: %d sequences', written.name, len(records))
+
+        return report
diff --git a/src/sourcerer/Sources/Ogrdb.py b/src/sourcerer/Sources/Ogrdb.py
new file mode 100644
index 0000000..57f1558
--- /dev/null
+++ b/src/sourcerer/Sources/Ogrdb.py
@@ -0,0 +1,404 @@
+"""
+OGRDB (AIRR Community germline sets)
+
+OGRDB publishes curated germline sets through a small REST API. Resolving a set
+to a download takes three calls -- species to a numeric id, id to the sets it
+holds, set to its latest release -- after which the FASTA is fetched twice, once
+ungapped and once IMGT-gapped, because a set carries V, D, J and C together and
+each segment is taken from the form that suits it.
+
+The segment split is the load-bearing piece and is kept verbatim from airrdb:
+V is taken gapped, so Change-O keeps the IMGT numbering it needs; D and J are
+taken ungapped; and constant regions are taken gapped. The one ambiguity is the
+delta locus, where the diversity segment IGHD and the constant IGHD share a
+name; they are told apart by length, since the constant is far longer. Getting
+this wrong silently files an allele under the wrong segment, so it is covered by
+fixtures.
+
+OGRDB serves only immunoglobulin sets, and only for the species it has curated,
+so a TR request or an uncovered locus resolves to nothing rather than an error.
+"""
+
+# Info
+__author__ = 'Ayelet Peres'
+
+# Imports
+import logging
+import re
+from datetime import UTC
+from urllib.parse import quote
+
+# Sourcerer imports
+from sourcerer.Exceptions import OgrdbParseError
+from sourcerer.Reference import (
+    KIND_CONSTANT,
+    KIND_VDJ,
+    ReferenceReport,
+    parseFasta,
+)
+from sourcerer.Sources.Base import DataUnit
+from sourcerer.Sources.Germline import ReferenceSource
+
+log = logging.getLogger(__name__)
+
+#: Endpoint.
+API = 'https://ogrdb.airr-community.org/api_v2'
+
+#: Species as OGRDB labels them.
+SPECIES_LABEL = {'human': 'Homo sapiens', 'mouse': 'Mus musculus'}
+
+#: The germline sets to fetch per species and locus, and the chains each covers.
+#: A locus can need more than one set: mouse splits V and J across strain-specific
+#: and all-strain sets. Kept as data so a new set is one line, not new code.
+SETS = {
+    ('human', 'IGH'): (('IGH_VDJ', ('IGHV', 'IGHD', 'IGHJ')), ('IGHC', ('IGHC',))),
+    ('human', 'IGK'): (('IGKappa_VJ', ('IGKV', 'IGKJ')),),
+    ('human', 'IGL'): (('IGLambda_VJ', ('IGLV', 'IGLJ')),),
+    ('mouse', 'IGH'): (('C57BL/6 IGH', ('IGHV', 'IGHD', 'IGHJ')),),
+    ('mouse', 'IGK'): (('C57BL/6J IGKV', ('IGKV',)),
+                       ('IGKJ (all strains)', ('IGKJ',))),
+    ('mouse', 'IGL'): (('C57BL/6J IGLV', ('IGLV',)),
+                       ('IGLJ (all strains)', ('IGLJ',))),
+}
+
+#: Loci OGRDB covers, offered as a filter flag.
+LOCI = ('IGH', 'IGK', 'IGL')
+
+#: The two forms each set is fetched in.
+FORMATS = ('ungapped', 'gapped')
+
+#: A constant region under 100 nucleotides is really the delta diversity segment
+#: wearing the same name; see the module docstring.
+CONSTANT_MIN_LENGTH = 100
+
+
+def normalizeVersion(value):
+    """
+    Render a release version without a trailing '.0'.
+
+    OGRDB reports the version as a number, so an integer release arrives as
+    '3.0' where the download URL wants '3'.
+
+    Arguments:
+      value: the reported version.
+
+    Returns:
+      str: the version as it appears in a download URL.
+    """
+    text = str(value)
+
+    return text[:-2] if text.endswith('.0') else text
+
+
+def safeSetName(set_name):
+    """
+    Make a filesystem-safe token from a set name.
+
+    Set names carry spaces and slashes (``C57BL/6J IGKV``) that must not become
+    directory separators in the raw mirror, but the name is never parsed back:
+    the real set name travels in the unit metadata.
+
+    Arguments:
+      set_name (str): the OGRDB set name.
+
+    Returns:
+      str: an identifier-safe token.
+    """
+    return re.sub(r'[^A-Za-z0-9]+', '_', set_name).strip('_')
+
+
+def bucketChain(name, sequence):
+    """
+    Decide which reference chain a germline allele belongs to.
+
+    V, D and J are filed under their four-character chain (``IGHV``); a constant
+    allele is filed under its locus constant (``IGHM`` -> ``IGHC``) so every
+    isotype of a locus lands in one file, as airrflow expects. Returns None for a
+    name too short to classify.
+
+    Arguments:
+      name (str): the allele name.
+      sequence (str): its sequence, used only to tell the delta segments apart.
+
+    Returns:
+      tuple: (chain, kind) or None.
+    """
+    if len(name) < 4:
+        return None
+
+    segment = name[3]
+    if segment == 'V':
+        return name[:4], KIND_VDJ
+    if segment == 'J':
+        return name[:4], KIND_VDJ
+    if segment == 'D':
+        # IGHD is both the diversity segment (short) and the delta constant
+        # (long); length is the only thing that separates them.
+        if len(sequence.replace('.', '')) < CONSTANT_MIN_LENGTH:
+            return name[:4], KIND_VDJ
+        return name[:3] + 'C', KIND_CONSTANT
+
+    return name[:3] + 'C', KIND_CONSTANT
+
+
+def _splitSegments(forms):
+    """
+    Sort a set's alleles into reference chains, form by form.
+
+    V is taken from the gapped alleles, so its IMGT numbering survives, along
+    with the constant regions; D and J are taken from the ungapped alleles. An
+    allele that appears in both forms is therefore filed once, from the form its
+    segment is read from.
+
+    Arguments:
+      forms (dict): 'ungapped' and 'gapped' each mapping allele name to sequence.
+
+    Returns:
+      dict: (chain, kind) to a list of (name, sequence) tuples.
+    """
+    chains = {}
+    for name, sequence in forms.get('gapped', {}).items():
+        target = bucketChain(name, sequence)
+        if target is None:
+            continue
+        chain, kind = target
+        if kind == KIND_CONSTANT or chain[3] == 'V':
+            chains.setdefault(target, []).append((name, sequence))
+
+    for name, sequence in forms.get('ungapped', {}).items():
+        target = bucketChain(name, sequence)
+        if target is None:
+            continue
+        chain, kind = target
+        if kind == KIND_VDJ and chain[3] in ('D', 'J'):
+            chains.setdefault(target, []).append((name, sequence))
+
+    return chains
+
+
+class OgrdbSource(ReferenceSource):
+    """
+    The OGRDB (AIRR Community) germline reference source.
+    """
+
+    name = 'ogrdb'
+    aliases = ('airrc',)
+    prefix = 'airrc'
+    description = 'OGRDB: AIRR Community curated immunoglobulin germline sets'
+    homepage = 'https://ogrdb.airr-community.org/'
+    collections = ('human', 'mouse')
+    collection_help = {'human': 'Homo sapiens curated IG sets',
+                       'mouse': 'Mus musculus curated IG sets'}
+
+    license = 'CC BY 4.0 (https://creativecommons.org/licenses/by/4.0/)'
+    citation = (
+        'Lees WD, Busse CE, Corcoran M, et al. OGRDB: a reference database of '
+        'inferred immune receptor genes. Nucleic Acids Res. '
+        '2020;48(D1):D964-D970. doi:10.1093/nar/gkz822',
+    )
+
+    # -- API client -------------------------------------------------------
+
+    def speciesId(self, species_label):
+        """
+        Resolve a species label to its OGRDB id.
+
+        Arguments:
+          species_label (str): the label, e.g. 'Homo sapiens'.
+
+        Returns:
+          str: the species id.
+
+        Raises:
+          OgrdbParseError: if the species is not listed.
+        """
+        payload = self.client.get('%s/germline/species' % API).json()
+        for item in payload.get('species', []):
+            if item.get('label') == species_label:
+                return item['id']
+
+        raise OgrdbParseError('OGRDB does not list species %r; the species '
+                              'endpoint changed or the species was withdrawn'
+                              % species_label)
+
+    def resolveSetId(self, species_id, locus, set_name):
+        """
+        Resolve a set name to its germline set id.
+
+        Arguments:
+          species_id (str): the OGRDB species id.
+          locus (str): the locus, e.g. 'IGH'.
+          set_name (str): the set name.
+
+        Returns:
+          str: the germline set id.
+
+        Raises:
+          OgrdbParseError: if the set is not found for the species and locus.
+        """
+        payload = self.client.get('%s/germline/sets/%s' % (API, species_id)).json()
+        for item in payload.get('germline_species', []):
+            if (item.get('germline_set_name') == set_name
+                    and item.get('locus') == locus):
+                return item['germline_set_id']
+
+        raise OgrdbParseError("OGRDB has no set %r for locus %s; the set was "
+                              'renamed or withdrawn' % (set_name, locus))
+
+    def latestRelease(self, set_id):
+        """
+        Read the latest release version and date of a set.
+
+        Arguments:
+          set_id (str): the germline set id.
+
+        Returns:
+          tuple: (version, release_date) with the date truncated to YYYY-MM-DD.
+
+        Raises:
+          OgrdbParseError: if the release payload is not shaped as expected.
+        """
+        safe = quote(set_id, safe='.')
+        payload = self.client.get('%s/germline/set/%s/latest' % (API, safe)).json()
+        try:
+            record = payload['GermlineSet'][0]
+
+            return normalizeVersion(record['release_version']), \
+                record['release_date'][:10]
+        except (KeyError, IndexError, TypeError) as error:
+            raise OgrdbParseError('OGRDB latest-release payload for %s is not '
+                                  'shaped as expected (%s)' % (set_id, error))
+
+    def fastaUrl(self, set_id, version, fmt, human):
+        """
+        Build a set's FASTA download URL.
+
+        Arguments:
+          set_id (str): the germline set id.
+          version (str): the release version.
+          fmt (str): 'ungapped' or 'gapped'.
+          human (bool): whether the species is human, which takes the ``_ex``
+            endpoint variant.
+
+        Returns:
+          str: the absolute download URL.
+        """
+        safe = quote(set_id, safe='.')
+        suffix = '_ex' if human else ''
+
+        return '%s/germline/set/%s/%s/%s%s' % (API, safe, version, fmt, suffix)
+
+    # -- schema -----------------------------------------------------------
+
+    def harvestSchema(self):
+        """
+        Contact OGRDB and build a fresh snapshot of the loci it curates.
+
+        The species and sets endpoints are queried, so a refresh both verifies
+        the API is answering and records which of the loci this source consumes
+        are actually available -- the drift signal that matters for OGRDB.
+
+        Returns:
+          SourceSchema: the harvested snapshot.
+        """
+        from datetime import datetime
+
+        from sourcerer.Schema import Collection, Field, SourceSchema
+        from sourcerer.Version import __version__
+
+        collections = {}
+        for sp in self.collections:
+            species_id = self.speciesId(SPECIES_LABEL[sp])
+            payload = self.client.get('%s/germline/sets/%s'
+                                      % (API, species_id)).json()
+            available = {x.get('locus') for x in payload.get('germline_species', [])}
+            loci = tuple(x for x in LOCI if x in available)
+            collections[sp] = Collection(name=sp,
+                                         fields=(Field(name='locus', values=loci),))
+
+        return SourceSchema(
+            source=self.name,
+            harvested=datetime.now(UTC).strftime('%Y-%m-%dT%H:%M:%SZ'),
+            harvested_by='sourcerer %s' % __version__,
+            source_urls={'api': API},
+            parse_contracts={'segment_split': 'V,C gapped; D,J ungapped; '
+                                              'delta D vs C by length'},
+            collections=collections)
+
+    # -- search and build -------------------------------------------------
+
+    def searchUnits(self, query):
+        """
+        Resolve a query to the germline files to fetch.
+
+        Each set is fetched twice, ungapped and gapped, so both are present when
+        the segments are split in buildReference. The set id and latest version
+        are resolved here so the download URLs are concrete.
+
+        Arguments:
+          query (Query): the validated request; collection is the species and the
+            locus filter narrows which sets are fetched.
+
+        Returns:
+          list: DataUnit objects, two per set.
+        """
+        species = query.collection
+        label = SPECIES_LABEL[species]
+        human = species == 'human'
+        locus_filter = query.filters.get('locus', '*')
+
+        species_id = self.speciesId(label)
+        units = []
+        for (set_species, locus), sets in SETS.items():
+            if set_species != species:
+                continue
+            if locus_filter not in ('*', locus):
+                continue
+
+            for set_name, chains in sets:
+                set_id = self.resolveSetId(species_id, locus, set_name)
+                version, _date = self.latestRelease(set_id)
+                for fmt in FORMATS:
+                    units.append(DataUnit(
+                        unit_id='%s.%s.fasta' % (safeSetName(set_name), fmt),
+                        collection=species,
+                        url=self.fastaUrl(set_id, version, fmt, human),
+                        metadata={'species': species, 'locus': locus,
+                                  'set_name': set_name, 'format': fmt,
+                                  'set_id': set_id, 'version': version,
+                                  'chains': list(chains)}))
+
+        if query.limit is not None:
+            units = units[:query.limit]
+
+        return units
+
+    def buildReference(self, entries, reference_dir):
+        """
+        Split the downloaded sets into per-chain reference FASTAs.
+
+        Arguments:
+          entries (list): (DataUnit, Path) pairs from the fetch step.
+          reference_dir (Path): the reference_base root.
+
+        Returns:
+          ReferenceReport: the files written.
+        """
+        report = ReferenceReport()
+        for species in sorted({unit.metadata['species'] for unit, _ in entries}):
+            forms = {fmt: {} for fmt in FORMATS}
+            for unit, path in entries:
+                if unit.metadata['species'] != species:
+                    continue
+                fmt = unit.metadata['format']
+                for header, sequence in parseFasta(path.read_text()):
+                    forms[fmt][header.split()[0]] = sequence
+
+            chains = _splitSegments(forms)
+            for (chain, kind), records in sorted(chains.items()):
+                written = self.writeChain(reference_dir, species, kind, chain,
+                                          records)
+                report.written.append((chain, written))
+                log.info('%s: %d sequences', written.name, len(records))
+
+        return report
diff --git a/src/sourcerer/Sources/__init__.py b/src/sourcerer/Sources/__init__.py
index 6dbe609..b9c3542 100644
--- a/src/sourcerer/Sources/__init__.py
+++ b/src/sourcerer/Sources/__init__.py
@@ -9,18 +9,40 @@
 __author__ = 'Susanna Marquez'
 
 # Sourcerer imports
+from sourcerer.Sources.AirrcImgt import AirrcImgtSource
+from sourcerer.Sources.Imgt import ImgtSource
 from sourcerer.Sources.Oas import OasSource
+from sourcerer.Sources.Ogrdb import OgrdbSource
 
-#: Every source sourcerer knows about, by commandline name.
-REGISTRY = {OasSource.name: OasSource}
+#: Every source sourcerer knows about, by canonical commandline name.
+REGISTRY = {source.name: source
+            for source in (OasSource, ImgtSource, OgrdbSource, AirrcImgtSource)}
+
+#: Alternative names that resolve to a canonical source, e.g. 'airrc' -> 'ogrdb'.
+ALIASES = {alias: source.name
+           for source in REGISTRY.values()
+           for alias in source.aliases}
+
+
+def canonicalName(name):
+    """
+    Resolve an alias to the canonical source name, or return it unchanged.
+
+    Arguments:
+      name (str): a source name or alias.
+
+    Returns:
+      str: the canonical source name.
+    """
+    return ALIASES.get(name, name)
 
 
 def getSource(name, client, schema=None):
     """
-    Instantiate a source by name.
+    Instantiate a source by name or alias.
 
     Arguments:
-      name (str): the source name.
+      name (str): the source name or alias.
       client (HttpClient): the shared HTTP client.
       schema (SourceSchema): a preloaded snapshot, or None to load on demand.
 
@@ -28,8 +50,9 @@ def getSource(name, client, schema=None):
       SourceBase: the source.
 
     Raises:
-      KeyError: if the name is not registered.
+      KeyError: if the name is not a known source or alias.
     """
+    name = canonicalName(name)
     if name not in REGISTRY:
         raise KeyError("unknown source '%s'; known sources: %s"
                        % (name, ', '.join(sorted(REGISTRY))))
diff --git a/src/sourcerer/data/schemas/airrc-imgt/schema.yaml b/src/sourcerer/data/schemas/airrc-imgt/schema.yaml
new file mode 100644
index 0000000..c5c13f7
--- /dev/null
+++ b/src/sourcerer/data/schemas/airrc-imgt/schema.yaml
@@ -0,0 +1,17 @@
+collections:
+  human:
+    fields: []
+    reported_totals: {}
+  mouse:
+    fields: []
+    reported_totals: {}
+field_aliases: {}
+harvested: '2026-08-06T00:00:00Z'
+harvested_by: sourcerer 0.1.0
+parse_contracts: {}
+schema_version: 1
+source: airrc-imgt
+source_urls:
+  imgt: imgt
+  ogrdb: ogrdb
+url_rules: {}
diff --git a/src/sourcerer/data/schemas/imgt/schema.yaml b/src/sourcerer/data/schemas/imgt/schema.yaml
new file mode 100644
index 0000000..1ebbfc9
--- /dev/null
+++ b/src/sourcerer/data/schemas/imgt/schema.yaml
@@ -0,0 +1,57 @@
+collections:
+  human:
+    fields:
+    - name: locus
+      pseudo_values: false
+      values:
+      - IGH
+      - IGK
+      - IGL
+      - TRA
+      - TRB
+      - TRG
+      - TRD
+      wildcard: '*'
+    - name: segment
+      pseudo_values: false
+      values:
+      - V
+      - D
+      - J
+      - C
+      wildcard: '*'
+    reported_totals: {}
+  mouse:
+    fields:
+    - name: locus
+      pseudo_values: false
+      values:
+      - IGH
+      - IGK
+      - IGL
+      - TRA
+      - TRB
+      - TRG
+      - TRD
+      wildcard: '*'
+    - name: segment
+      pseudo_values: false
+      values:
+      - V
+      - D
+      - J
+      - C
+      wildcard: '*'
+    reported_totals: {}
+field_aliases: {}
+harvested: '2026-08-06T00:00:00Z'
+harvested_by: sourcerer 0.1.0
+parse_contracts:
+  fasta_block: second 
 element
+  header: pipe-delimited, allele in field 2
+schema_version: 1
+source: imgt
+source_urls:
+  genelect: https://www.imgt.org/genedb/GENElect
+  release: https://www.imgt.org/download/GENE-DB/RELEASE
+url_rules: {}
diff --git a/src/sourcerer/data/schemas/ogrdb/schema.yaml b/src/sourcerer/data/schemas/ogrdb/schema.yaml
new file mode 100644
index 0000000..0171045
--- /dev/null
+++ b/src/sourcerer/data/schemas/ogrdb/schema.yaml
@@ -0,0 +1,31 @@
+collections:
+  human:
+    fields:
+    - name: locus
+      pseudo_values: false
+      values:
+      - IGH
+      - IGK
+      - IGL
+      wildcard: '*'
+    reported_totals: {}
+  mouse:
+    fields:
+    - name: locus
+      pseudo_values: false
+      values:
+      - IGH
+      - IGK
+      - IGL
+      wildcard: '*'
+    reported_totals: {}
+field_aliases: {}
+harvested: '2026-08-06T00:00:00Z'
+harvested_by: sourcerer 0.1.0
+parse_contracts:
+  segment_split: V,C gapped; D,J ungapped; delta D vs C by length
+schema_version: 1
+source: ogrdb
+source_urls:
+  api: https://ogrdb.airr-community.org/api_v2
+url_rules: {}
diff --git a/tests/data/README.md b/tests/data/README.md
index 8de0a7c..daff1c1 100644
--- a/tests/data/README.md
+++ b/tests/data/README.md
@@ -43,3 +43,28 @@ The two paired fixtures are both required. They are **not** the same schema:
 only the 180 column file would leave the common case uncovered. The
 `csv_paired/` fixture also has no run accession in its filename, which is what
 pins the rule that unit identifiers are opaque.
+
+## IMGT
+
+Derived from the live IMGT/GENE-DB GENElect service on **2026-08-06**. IMGT data
+is subject to the IMGT terms of use ()
+and its use requires citing IMGT, the international ImMunoGeneTics information
+system (Lefranc MP et al., *Nucleic Acids Res.* 2015). The excerpts are trimmed to
+the smallest form that exercises the parser.
+
+| File | Content |
+|---|---|
+| `imgt_ighd.html` | A GENElect reply reduced to its two `
` blocks — the query echo and three real human IGHD records — so `extractFasta` reads the second block. |
+| `imgt_error.html` | A hand-written stand-in for an IMGT error page: HTTP 200 with a single `
` and no FASTA, which is why validity cannot be the status code alone. |
+
+## OGRDB
+
+Trimmed from the live OGRDB `api_v2` human `IGKappa_VJ` set on **2026-08-06**.
+OGRDB data is distributed under **CC BY 4.0**; cite Lees WD et al., *Nucleic Acids
+Res.* 2020. Both forms of the same set are kept because the segment split reads V
+from one and J from the other.
+
+| File | Content |
+|---|---|
+| `ogrdb_igk_ungapped.fasta` | Two IGKV and two IGKJ alleles, ungapped. J is taken from here. |
+| `ogrdb_igk_gapped.fasta` | The same alleles IMGT-gapped; the IGKV records carry `.` gaps. V is taken from here, which is what keeps its numbering. |
diff --git a/tests/data/imgt_error.html b/tests/data/imgt_error.html
new file mode 100644
index 0000000..34d9fc1
--- /dev/null
+++ b/tests/data/imgt_error.html
@@ -0,0 +1,4 @@
+
+

IMGT/GENE-DB

+
No result for your query.
+ diff --git a/tests/data/imgt_ighd.html b/tests/data/imgt_ighd.html new file mode 100644 index 0000000..d2da1bd --- /dev/null +++ b/tests/data/imgt_ighd.html @@ -0,0 +1,14 @@ + +IMGT/GENE-DB +

IMGT/GENE-DB reference sequences

+
Homo sapiens IGHD: query 7.14
+

Result

+
+>X97051|IGHD1-1*01|Homo sapiens|F|D-REGION|33714..33730|17 nt|1| | | | |17+0=17| | |
+ggtacaactggaacgac
+>X13972|IGHD1-14*01|Homo sapiens|ORF|D-REGION|14518..14534|17 nt|1| | | | |17+0=17| | |
+ggtataaccggaaccac
+>X97051|IGHD1-20*01|Homo sapiens|F|D-REGION|62015..62031|17 nt|1| | | | |17+0=17| | |
+ggtataactggaacgac
+
+ diff --git a/tests/data/ogrdb_igk_gapped.fasta b/tests/data/ogrdb_igk_gapped.fasta new file mode 100644 index 0000000..4074415 --- /dev/null +++ b/tests/data/ogrdb_igk_gapped.fasta @@ -0,0 +1,8 @@ +>IGKV1-12*01 +GACATCCAGATGACCCAGTCTCCATCTTCCGTGTCTGCATCTGTAGGAGACAGAGTCACCATCACTTGTCGGGCGAGTCAGGGTATT..................AGCAGCTGGTTAGCCTGGTATCAGCAGAAACCAGGGAAAGCCCCTAAGCTCCTGATCTATGCTGCA.....................TCCAGTTTGCAAAGTGGGGTCCCA...TCAAGGTTCAGCGGCAGTGGA......TCTGGGACAGATTTCACTCTCACCATCAGCAGCCTGCAGCCTGAAGATTTTGCAACTTACTATTGTCAACAGGCTAACAGTTTCCCTCC +>IGKV1-13*01 +GCCATCCAGTTGACCCAGTCTCCATCCTCCCTGTCTGCATCTGTAGGAGACAGAGTCACCATCACTTGCCGGGCAAGTCAGGGCATT..................AGCAGTGCTTTAGCCTGATATCAGCAGAAACCAGGGAAAGCTCCTAAGCTCCTGATCTATGATGCC.....................TCCAGTTTGGAAAGTGGGGTCCCA...TCAAGGTTCAGCGGCAGTGGA......TCTGGGACAGATTTCACTCTCACCATCAGCAGCCTGCAGCCTGAAGATTTTGCAACTTATTACTGTCAACAGTTTAATAATTACCCTCA +>IGKJ1*01 +GTGGACGTTCGGCCAAGGGACCAAGGTGGAAATCAAAC +>IGKJ2*01 +TGTACACTTTTGGCCAGGGGACCAAGCTGGAGATCAAAC diff --git a/tests/data/ogrdb_igk_ungapped.fasta b/tests/data/ogrdb_igk_ungapped.fasta new file mode 100644 index 0000000..2453dff --- /dev/null +++ b/tests/data/ogrdb_igk_ungapped.fasta @@ -0,0 +1,8 @@ +>IGKV1-12*01 +GACATCCAGATGACCCAGTCTCCATCTTCCGTGTCTGCATCTGTAGGAGACAGAGTCACCATCACTTGTCGGGCGAGTCAGGGTATTAGCAGCTGGTTAGCCTGGTATCAGCAGAAACCAGGGAAAGCCCCTAAGCTCCTGATCTATGCTGCATCCAGTTTGCAAAGTGGGGTCCCATCAAGGTTCAGCGGCAGTGGATCTGGGACAGATTTCACTCTCACCATCAGCAGCCTGCAGCCTGAAGATTTTGCAACTTACTATTGTCAACAGGCTAACAGTTTCCCTCC +>IGKV1-13*01 +GCCATCCAGTTGACCCAGTCTCCATCCTCCCTGTCTGCATCTGTAGGAGACAGAGTCACCATCACTTGCCGGGCAAGTCAGGGCATTAGCAGTGCTTTAGCCTGATATCAGCAGAAACCAGGGAAAGCTCCTAAGCTCCTGATCTATGATGCCTCCAGTTTGGAAAGTGGGGTCCCATCAAGGTTCAGCGGCAGTGGATCTGGGACAGATTTCACTCTCACCATCAGCAGCCTGCAGCCTGAAGATTTTGCAACTTATTACTGTCAACAGTTTAATAATTACCCTCA +>IGKJ1*01 +GTGGACGTTCGGCCAAGGGACCAAGGTGGAAATCAAAC +>IGKJ2*01 +TGTACACTTTTGGCCAGGGGACCAAGCTGGAGATCAAAC diff --git a/tests/test_AirrcImgt.py b/tests/test_AirrcImgt.py new file mode 100644 index 0000000..5d792d1 --- /dev/null +++ b/tests/test_AirrcImgt.py @@ -0,0 +1,119 @@ +""" +Unit tests for the airrc-imgt blended source +""" + +# Info +__author__ = 'Ayelet Peres' + +# Imports +import os +import tempfile +import unittest +from pathlib import Path + +# Sourcerer imports +from sourcerer.Sources.AirrcImgt import AirrcImgtSource +from sourcerer.Sources.Base import DataUnit, Query +from tests.test_ogrdb import StubClient + +test_path = os.path.dirname(os.path.realpath(__file__)) +data_path = os.path.join(test_path, 'data') + + +def readFixture(name): + """Read a captured fixture from tests/data.""" + with open(os.path.join(data_path, name)) as handle: + return handle.read() + + +class TestSearchUnits(unittest.TestCase): + """ + Tests for composing the OGRDB and IMGT halves of the blend + """ + + def setUp(self): + self.source = AirrcImgtSource(client=StubClient()) + self.units = self.source.searchUnits(Query(collection='human')) + self.imgt = [u for u in self.units if u.metadata['via'] == 'imgt'] + self.ogrdb = [u for u in self.units if u.metadata['via'] == 'ogrdb'] + + def test_both_sources_contribute(self): + """The blend draws from OGRDB and IMGT, each tagged with its origin.""" + self.assertTrue(self.ogrdb) + self.assertTrue(self.imgt) + self.assertEqual({u.metadata['via'] for u in self.units}, + {'ogrdb', 'imgt'}) + + def test_immunoglobulin_vdj_comes_only_from_ogrdb(self): + """No IMGT unit supplies an immunoglobulin V, D or J: those are OGRDB's.""" + ig_vdj = [u for u in self.imgt + if u.metadata['kind'] == 'vdj' + and not u.metadata['locus'].startswith('TR')] + self.assertEqual(ig_vdj, []) + + def test_imgt_fills_tr_and_the_light_constants(self): + """IMGT supplies all TR, and the IG constants OGRDB has no set for.""" + loci = {u.metadata['locus'] for u in self.imgt} + self.assertEqual(loci & {'TRA', 'TRB', 'TRG', 'TRD'}, + {'TRA', 'TRB', 'TRG', 'TRD'}) + ig_constants = {u.metadata['chain'] for u in self.imgt + if u.metadata['kind'] == 'constant' + and not u.metadata['locus'].startswith('TR')} + # Human IGHC is OGRDB's; IMGT provides only the light constants. + self.assertEqual(ig_constants, {'IGKC', 'IGLC'}) + + def test_no_amino_acid_in_the_blend(self): + """Amino acid V is not part of the airrc-imgt blend.""" + self.assertNotIn('vdj_aa', {u.metadata['kind'] for u in self.imgt}) + + def test_mouse_takes_all_ig_constants_from_imgt(self): + """Mouse has no OGRDB constant set, so all three come from IMGT.""" + # The IMGT gap is pure IMGT selection, so it needs no OGRDB call. + gap = self.source._imgtGapUnits('mouse') + ig_constants = {u.metadata['chain'] for u in gap + if u.metadata['kind'] == 'constant' + and not u.metadata['locus'].startswith('TR')} + self.assertEqual(ig_constants, {'IGHC', 'IGKC', 'IGLC'}) + + +class TestBuildReference(unittest.TestCase): + """ + Tests for routing each unit back to the source that produced it + """ + + def test_each_source_writes_its_own_prefix(self): + """OGRDB units land as airrc_ files and IMGT units as imgt_ files.""" + with tempfile.TemporaryDirectory() as tmp: + tmp = Path(tmp) + source = AirrcImgtSource(client=StubClient()) + + ogrdb_paths = {} + for fmt in ('ungapped', 'gapped'): + path = tmp / ('igk_%s.fasta' % fmt) + path.write_text(readFixture('ogrdb_igk_%s.fasta' % fmt)) + ogrdb_paths[fmt] = path + ogrdb_entries = [ + (DataUnit(unit_id='IGKappa_VJ.%s.fasta' % fmt, collection='human', + url='x', metadata={'species': 'human', 'locus': 'IGK', + 'set_name': 'IGKappa_VJ', 'format': fmt, + 'via': 'ogrdb'}), path) + for fmt, path in ogrdb_paths.items()] + + imgt_page = tmp / 'IGHD.html' + imgt_page.write_text(readFixture('imgt_ighd.html')) + imgt_entries = [ + (DataUnit(unit_id='vdj/IGHD.html', collection='human', url='x', + metadata={'species': 'human', 'chain': 'IGHD', + 'kind': 'vdj', 'locus': 'IGH', 'segment': 'D', + 'via': 'imgt'}), imgt_page)] + + source.buildReference(ogrdb_entries + imgt_entries, + tmp / 'reference_base') + + vdj = tmp / 'reference_base' / 'human' / 'vdj' + self.assertTrue((vdj / 'airrc_human_IGKV.fasta').exists()) + self.assertTrue((vdj / 'imgt_human_IGHD.fasta').exists()) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/test_Imgt.py b/tests/test_Imgt.py new file mode 100644 index 0000000..5967dd3 --- /dev/null +++ b/tests/test_Imgt.py @@ -0,0 +1,143 @@ +""" +Unit tests for the IMGT source +""" + +# Info +__author__ = 'Ayelet Peres' + +# Imports +import os +import tempfile +import unittest +from pathlib import Path + +# Sourcerer imports +from sourcerer.Exceptions import ImgtParseError +from sourcerer.Sources.Base import DataUnit, Query +from sourcerer.Sources.Imgt import ( + ImgtSource, + buildQueryUrl, + extractFasta, + isValidResponse, +) + +test_path = os.path.dirname(os.path.realpath(__file__)) +data_path = os.path.join(test_path, 'data') + + +def readFixture(name): + """Read a captured fixture from tests/data.""" + with open(os.path.join(data_path, name)) as handle: + return handle.read() + + +class TestQueryUrl(unittest.TestCase): + """ + Tests for GENElect URL construction + """ + + def test_encodes_query_and_species(self): + """The query number is escaped and the species is sent pre-encoded.""" + url = buildQueryUrl('human', '7.14', 'IGHV') + self.assertEqual( + url, + 'https://www.imgt.org/genedb/GENElect?query=7.14+IGHV' + '&species=Homo%20sapiens') + + def test_appends_label(self): + """An IMGTlabel qualifier is appended when given.""" + url = buildQueryUrl('mouse', '8.1', 'IGHV', label='L-PART1+L-PART2') + self.assertTrue(url.endswith('&IMGTlabel=L-PART1+L-PART2')) + + +class TestResponseParsing(unittest.TestCase): + """ + Tests for reading a GENElect reply + """ + + def test_valid_response_has_second_pre_with_fasta(self): + """A real reply has a second
 block carrying a FASTA."""
+        self.assertTrue(isValidResponse(readFixture('imgt_ighd.html')))
+
+    def test_error_page_is_not_valid_despite_http_200(self):
+        """An error page with a single 
 is rejected."""
+        self.assertFalse(isValidResponse(readFixture('imgt_error.html')))
+
+    def test_extract_returns_fasta_with_underscored_species(self):
+        """extractFasta pulls the FASTA and underscores the species name."""
+        fasta = extractFasta(readFixture('imgt_ighd.html'), 'human')
+        self.assertIn('>X97051|IGHD1-1*01|Homo_sapiens|F', fasta)
+        self.assertNotIn('Homo sapiens', fasta)
+        self.assertEqual(fasta.count('>'), 3)
+
+    def test_extract_raises_on_error_page(self):
+        """A page with no second 
 is a parse error, not an empty result."""
+        with self.assertRaises(ImgtParseError):
+            extractFasta(readFixture('imgt_error.html'), 'human')
+
+
+class TestSearchUnits(unittest.TestCase):
+    """
+    Tests for enumerating the germline files to fetch
+    """
+
+    def setUp(self):
+        self.source = ImgtSource(client=None)
+
+    def _query(self, **filters):
+        resolved = {'locus': '*', 'segment': '*'}
+        resolved.update(filters)
+        return Query(collection='human', filters=resolved)
+
+    def test_unfiltered_covers_vdj_constant_and_aa(self):
+        """With no filter every VDJ, constant and AA chain is scheduled."""
+        units = self.source.searchUnits(self._query())
+        # 17 VDJ + 7 constant + 7 AA V.
+        self.assertEqual(len(units), 31)
+        kinds = {u.metadata['kind'] for u in units}
+        self.assertEqual(kinds, {'vdj', 'constant', 'vdj_aa'})
+
+    def test_locus_and_segment_filter(self):
+        """--locus IGH --segment V leaves the nucleotide and amino acid V."""
+        units = self.source.searchUnits(self._query(locus='IGH', segment='V'))
+        self.assertEqual({u.metadata['chain'] for u in units}, {'IGHV'})
+        self.assertEqual({u.metadata['kind'] for u in units}, {'vdj', 'vdj_aa'})
+
+    def test_mouse_light_constant_uses_special_query(self):
+        """Mouse IGKC and IGLC take query 7.5, which 14.1 does not serve."""
+        source = ImgtSource(client=None)
+        units = source.searchUnits(
+            Query(collection='mouse',
+                  filters={'locus': 'IGK', 'segment': 'C'}))
+        self.assertEqual(len(units), 1)
+        self.assertIn('query=7.5+IGKC', units[0].url)
+
+
+class TestBuildReference(unittest.TestCase):
+    """
+    Tests for extracting downloaded pages into the reference tree
+    """
+
+    def test_writes_chain_fasta_from_page(self):
+        """A downloaded GENElect page becomes one per-chain reference FASTA."""
+        with tempfile.TemporaryDirectory() as tmp:
+            tmp = Path(tmp)
+            raw = tmp / 'IGHD.html'
+            raw.write_text(readFixture('imgt_ighd.html'))
+            unit = DataUnit(
+                unit_id='vdj/IGHD.html', collection='human', url='x',
+                metadata={'species': 'human', 'chain': 'IGHD',
+                          'kind': 'vdj', 'locus': 'IGH', 'segment': 'D'})
+
+            source = ImgtSource(client=None)
+            report = source.buildReference([(unit, raw)], tmp / 'reference_base')
+
+            written = (tmp / 'reference_base' / 'human' / 'vdj'
+                       / 'imgt_human_IGHD.fasta')
+            self.assertTrue(written.exists())
+            self.assertEqual(written.read_text().count('>'), 3)
+            self.assertEqual(len(report.written), 1)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/tests/test_Reference.py b/tests/test_Reference.py
new file mode 100644
index 0000000..d8a69ef
--- /dev/null
+++ b/tests/test_Reference.py
@@ -0,0 +1,255 @@
+"""
+Unit tests for the germline reference output layer
+"""
+
+# Info
+__author__ = 'Ayelet Peres'
+
+# Imports
+import io
+import shutil
+import tarfile
+import tempfile
+import unittest
+from pathlib import Path
+
+# Sourcerer imports
+from sourcerer import Reference
+from sourcerer.Exceptions import SourcererError
+
+HAS_MAKEBLASTDB = shutil.which('makeblastdb') is not None
+
+
+class Canned:
+    """A stand-in response exposing just .text for the directory-index mirror."""
+
+    def __init__(self, text=''):
+        self.text = text
+
+
+class EmptyIndexClient:
+    """An HTTP client whose directory listings are empty, so nothing mirrors."""
+
+    def get(self, url):
+        return Canned('')
+
+    def fetch(self, url, dest, **kwargs):
+        raise AssertionError('an empty index should not trigger a fetch')
+
+
+class TestFastaHelpers(unittest.TestCase):
+    """
+    Tests for the FASTA parse, name and clean helpers
+    """
+
+    def test_parse_round_trips_records(self):
+        """parseFasta reads headers and collapses wrapped sequence lines."""
+        text = '>a\nACGT\nACGT\n>b\nTTTT\n'
+        self.assertEqual(Reference.parseFasta(text),
+                         [('a', 'ACGTACGT'), ('b', 'TTTT')])
+
+    def test_parse_ignores_blank_lines(self):
+        """Blank lines between records are not mistaken for sequence."""
+        self.assertEqual(Reference.parseFasta('\n>a\n\nACGT\n\n'),
+                         [('a', 'ACGT')])
+
+    def test_allele_name_from_imgt_pipe_header(self):
+        """The allele name is the second pipe field of an IMGT header."""
+        header = 'X02897|IGHV1-2*02|Homo sapiens|F|V-REGION'
+        self.assertEqual(Reference.alleleName(header), 'IGHV1-2*02')
+
+    def test_allele_name_from_plain_header(self):
+        """A header with no pipe yields its first whitespace-delimited token."""
+        self.assertEqual(Reference.alleleName('IGKV1-12*01 extra'), 'IGKV1-12*01')
+
+    def test_clean_degaps_uppercases_and_dedups(self):
+        """cleanForBlast removes gaps, uppercases, and drops repeated names."""
+        records = [('X1|IGHV1-2*02|H', 'ac.gt'),
+                   ('X2|IGHV1-2*02|H', 'aaaa'),   # duplicate name, dropped
+                   ('IGHV3*01', 'gg..cc')]
+        self.assertEqual(
+            Reference.cleanForBlast(records),
+            [('IGHV1-2*02', 'ACGT'), ('IGHV3*01', 'GGCC')])
+
+
+class TestReferencePaths(unittest.TestCase):
+    """
+    Tests for where a chain's FASTA lands in the reference tree
+    """
+
+    def test_vdj_path(self):
+        """A V/D/J chain lands under vdj/ with the source prefix."""
+        path = Reference.referenceFastaPath('/ref', 'imgt', 'human',
+                                            Reference.KIND_VDJ, 'IGHV')
+        self.assertEqual(path, Path('/ref/human/vdj/imgt_human_IGHV.fasta'))
+
+    def test_constant_path(self):
+        """A constant chain lands under constant/."""
+        path = Reference.referenceFastaPath('/ref', 'airrc', 'human',
+                                            Reference.KIND_CONSTANT, 'IGHC')
+        self.assertEqual(path, Path('/ref/human/constant/airrc_human_IGHC.fasta'))
+
+    def test_amino_acid_path_is_disambiguated(self):
+        """Amino acid V carries an aa_ tag so it cannot collide with nucleotide V."""
+        path = Reference.referenceFastaPath('/ref', 'imgt', 'mouse',
+                                            Reference.KIND_AA, 'IGHV')
+        self.assertEqual(path,
+                         Path('/ref/mouse/vdj_aa/imgt_aa_mouse_IGHV.fasta'))
+
+
+class TestExtractTar(unittest.TestCase):
+    """
+    Tests for the path-traversal guard on archive extraction
+    """
+
+    def test_rejects_member_escaping_destination(self):
+        """A member pointing outside the destination is refused, not written."""
+        with tempfile.TemporaryDirectory() as tmp:
+            tmp = Path(tmp)
+            archive = tmp / 'evil.tar'
+            with tarfile.open(archive, 'w') as tar:
+                info = tarfile.TarInfo('../escaped.txt')
+                info.size = 3
+                tar.addfile(info, io.BytesIO(b'bad'))
+
+            dest = tmp / 'out'
+            dest.mkdir()
+            with self.assertRaises(SourcererError):
+                Reference.extractTar(archive, dest)
+            self.assertFalse((tmp / 'escaped.txt').exists())
+
+
+@unittest.skipUnless(HAS_MAKEBLASTDB, 'makeblastdb not on PATH')
+class TestBuildIgblastBase(unittest.TestCase):
+    """
+    Tests for aggregating a reference_base into BLAST databases
+    """
+
+    def _referenceBase(self, root):
+        vdj = root / 'human' / 'vdj'
+        vdj.mkdir(parents=True)
+        (vdj / 'imgt_human_IGHV.fasta').write_text('>IGHV1-2*02\nACGTACGTACGT\n')
+        (vdj / 'imgt_human_IGHJ.fasta').write_text('>IGHJ1*01\nTTTTGGGGCCCC\n')
+
+    def test_builds_present_and_skips_absent(self):
+        """Only chains with sequences build a database; the rest are skipped."""
+        with tempfile.TemporaryDirectory() as tmp:
+            tmp = Path(tmp)
+            reference = tmp / 'reference_base'
+            self._referenceBase(reference)
+
+            report = Reference.buildIgblastBase(reference, tmp / 'igblast_base',
+                                                EmptyIndexClient(),
+                                                species=['human'])
+
+            self.assertIn('human_ig_v', report.built)
+            self.assertIn('human_ig_j', report.built)
+            self.assertIn('human_ig_d', report.skipped_empty)
+            self.assertIn('human_tr_v', report.skipped_empty)
+            self.assertTrue((tmp / 'igblast_base' / 'database'
+                             / 'human_ig_v.nsq').exists())
+
+    def test_missing_makeblastdb_is_a_clear_error(self):
+        """runMakeblastdb names the missing binary rather than failing obscurely."""
+        original = shutil.which
+        try:
+            shutil.which = lambda name: None
+            with tempfile.TemporaryDirectory() as tmp:
+                fasta = Path(tmp) / 'x.fasta'
+                fasta.write_text('>a\nACGT\n')
+                with self.assertRaises(SourcererError) as caught:
+                    Reference.runMakeblastdb(fasta, Path(tmp) / 'x', 'nucl')
+                self.assertIn('makeblastdb', str(caught.exception))
+        finally:
+            shutil.which = original
+
+
+class TestParseReferenceName(unittest.TestCase):
+    """
+    Tests for reading (species, chain, is_aa) from a filename
+    """
+
+    def test_flat_name(self):
+        """A bare species_chain name is recognised."""
+        self.assertEqual(Reference.parseReferenceName('human_IGHV.fasta'),
+                         ('human', 'IGHV', False))
+
+    def test_prefix_is_ignored(self):
+        """A source prefix such as imgt_ or airrc_ is allowed and ignored."""
+        self.assertEqual(Reference.parseReferenceName('imgt_human_IGHV.fasta'),
+                         ('human', 'IGHV', False))
+        self.assertEqual(Reference.parseReferenceName('airrc_mouse_IGKC.fasta'),
+                         ('mouse', 'IGKC', False))
+
+    def test_aa_marker(self):
+        """aa_ marks a translated V, with or without a prefix."""
+        self.assertEqual(Reference.parseReferenceName('aa_human_IGHV.fasta'),
+                         ('human', 'IGHV', True))
+        self.assertEqual(Reference.parseReferenceName('imgt_aa_human_IGHV.fasta'),
+                         ('human', 'IGHV', True))
+
+    def test_rejects_unknown_species_or_chain_or_extension(self):
+        """A name that is not species_knownchain.fasta is not recognised."""
+        self.assertIsNone(Reference.parseReferenceName('rat_IGHV.fasta'))
+        self.assertIsNone(Reference.parseReferenceName('human_XYZ.fasta'))
+        self.assertIsNone(Reference.parseReferenceName('human_IGHV.txt'))
+        self.assertIsNone(Reference.parseReferenceName('notes.fasta'))
+
+
+class TestPlanReference(unittest.TestCase):
+    """
+    Tests for planning an IgBLAST build from a reference folder
+    """
+
+    def _write(self, path, name, body):
+        path.mkdir(parents=True, exist_ok=True)
+        (path / name).write_text(body)
+
+    def test_flat_and_nested_layouts_plan_the_same(self):
+        """The same files plan identically whether flat or nested."""
+        with tempfile.TemporaryDirectory() as tmp:
+            tmp = Path(tmp)
+            flat, nested = tmp / 'flat', tmp / 'nested'
+            self._write(flat, 'human_IGHV.fasta', '>IGHV1-2*02\nACGT\n')
+            self._write(flat, 'human_IGHJ.fasta', '>IGHJ1*01\nTTGG\n')
+            self._write(nested / 'human' / 'vdj', 'imgt_human_IGHV.fasta',
+                        '>IGHV1-2*02\nACGT\n')
+            self._write(nested / 'human' / 'vdj', 'imgt_human_IGHJ.fasta',
+                        '>IGHJ1*01\nTTGG\n')
+
+            built_flat = {b for b, _t, _r in
+                          Reference.planReference(flat).databases}
+            built_nested = {b for b, _t, _r in
+                            Reference.planReference(nested).databases}
+            self.assertEqual(built_flat, {'human_ig_v', 'human_ig_j'})
+            self.assertEqual(built_flat, built_nested)
+
+    def test_reports_empty_unrecognized_and_duplicates(self):
+        """The plan surfaces gaps, unknown names, and dropped duplicate names."""
+        with tempfile.TemporaryDirectory() as tmp:
+            tmp = Path(tmp)
+            self._write(tmp, 'human_IGHV.fasta',
+                        '>IGHV1-2*02\nACGT\n>IGHV1-2*02\nAAAA\n')  # duplicate name
+            self._write(tmp, 'notes.fasta', '>x\nACGT\n')          # unrecognized
+
+            plan = Reference.planReference(tmp)
+            self.assertEqual(plan.found_species, ['human'])
+            self.assertEqual(plan.duplicates.get('human_ig_v'), 1)
+            self.assertIn('human_ig_d', plan.empty)
+            self.assertEqual([p.name for p in plan.unrecognized], ['notes.fasta'])
+            self.assertTrue(plan.ok)
+
+    def test_species_filter(self):
+        """--species narrows the plan to the requested species."""
+        with tempfile.TemporaryDirectory() as tmp:
+            tmp = Path(tmp)
+            self._write(tmp, 'human_IGHV.fasta', '>IGHV1-2*02\nACGT\n')
+            self._write(tmp, 'mouse_IGHV.fasta', '>IGHV1*01\nACGT\n')
+
+            built = {b for b, _t, _r in
+                     Reference.planReference(tmp, species=['human']).databases}
+            self.assertEqual(built, {'human_ig_v'})
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/tests/test_live.py b/tests/test_live.py
new file mode 100644
index 0000000..91fb88f
--- /dev/null
+++ b/tests/test_live.py
@@ -0,0 +1,86 @@
+"""
+Live API checks for the reference sources
+
+These contact IMGT and OGRDB for real, so they are skipped unless SOURCERER_LIVE
+is set in the environment. The weekly check-apis workflow sets it and runs this
+module on its own; a red run there is the early warning that an upstream API
+changed shape before a user's download hits the same failure.
+
+The endpoint constants come from the source modules, so this tests exactly what
+production calls rather than a second copy of the URLs.
+"""
+
+# Info
+__author__ = 'Ayelet Peres'
+
+# Imports
+import os
+import unittest
+
+# Sourcerer imports
+from sourcerer.Http import HttpClient
+from sourcerer.Sources.Imgt import (
+    ImgtSource,
+    buildQueryUrl,
+    extractFasta,
+    isValidResponse,
+)
+from sourcerer.Sources.Ogrdb import OgrdbSource
+
+LIVE = os.environ.get('SOURCERER_LIVE')
+
+
+@unittest.skipUnless(LIVE, 'set SOURCERER_LIVE=1 to contact IMGT and OGRDB')
+class TestImgtLive(unittest.TestCase):
+    """
+    Live checks against IMGT/GENE-DB
+    """
+
+    def setUp(self):
+        self.source = ImgtSource(client=HttpClient())
+
+    def test_genelect_returns_a_germline_fasta(self):
+        """A GENElect query returns a page with a real FASTA in it."""
+        url = buildQueryUrl('human', '7.14', 'IGHD')
+        html = self.source.client.get(url).text
+        self.assertTrue(isValidResponse(html),
+                        'GENElect no longer returns a second 
 with a FASTA')
+        self.assertIn('>', extractFasta(html, 'human'))
+
+    def test_release_tag_is_readable(self):
+        """The GENE-DB release tag is still published and non-empty."""
+        self.assertTrue(self.source.fetchRelease(),
+                        'the IMGT release tag could not be read')
+
+
+@unittest.skipUnless(LIVE, 'set SOURCERER_LIVE=1 to contact IMGT and OGRDB')
+class TestOgrdbLive(unittest.TestCase):
+    """
+    Live checks against OGRDB
+    """
+
+    def setUp(self):
+        self.source = OgrdbSource(client=HttpClient())
+
+    def test_harvest_schema_sees_the_consumed_loci(self):
+        """species and sets resolve, and human still exposes IGH, IGK and IGL."""
+        schema = self.source.harvestSchema()
+        human = schema.getCollection('human').getField('locus')
+        for locus in ('IGH', 'IGK', 'IGL'):
+            self.assertIn(locus, human.values,
+                          'OGRDB no longer lists %s for human' % locus)
+
+    def test_set_resolves_to_a_downloadable_fasta(self):
+        """A set resolves to a release whose FASTA download is non-empty."""
+        from sourcerer.Sources.Base import Query
+
+        units = self.source.searchUnits(
+            Query(collection='human', filters={'locus': 'IGK'}))
+        self.assertTrue(units, 'no OGRDB units resolved for human IGK')
+
+        body = self.source.client.get(units[0].url).text
+        self.assertIn('>', body, 'the OGRDB FASTA download was empty')
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/tests/test_ogrdb.py b/tests/test_ogrdb.py
new file mode 100644
index 0000000..0ee4d59
--- /dev/null
+++ b/tests/test_ogrdb.py
@@ -0,0 +1,198 @@
+"""
+Unit tests for the OGRDB source
+"""
+
+# Info
+__author__ = 'Ayelet Peres'
+
+# Imports
+import os
+import tempfile
+import unittest
+from pathlib import Path
+
+# Sourcerer imports
+from sourcerer.Reference import KIND_CONSTANT, KIND_VDJ
+from sourcerer.Sources.Base import DataUnit, Query
+from sourcerer.Sources.Ogrdb import (
+    OgrdbSource,
+    bucketChain,
+    normalizeVersion,
+    safeSetName,
+)
+
+test_path = os.path.dirname(os.path.realpath(__file__))
+data_path = os.path.join(test_path, 'data')
+
+
+def readFixture(name):
+    """Read a captured fixture from tests/data."""
+    with open(os.path.join(data_path, name)) as handle:
+        return handle.read()
+
+
+class Canned:
+    """A response exposing the .json() the OGRDB client reads."""
+
+    def __init__(self, payload):
+        self._payload = payload
+
+    def json(self):
+        return self._payload
+
+
+class StubClient:
+    """An OGRDB API client backed by canned payloads, no network."""
+
+    SPECIES = {'species': [{'label': 'Homo sapiens', 'id': '9606'}]}
+    SETS = {'germline_species': [
+        {'germline_set_name': 'IGH_VDJ', 'locus': 'IGH',
+         'germline_set_id': '9606.IGH_VDJ'},
+        {'germline_set_name': 'IGHC', 'locus': 'IGH',
+         'germline_set_id': '9606.IGHC'},
+        {'germline_set_name': 'IGKappa_VJ', 'locus': 'IGK',
+         'germline_set_id': '9606.IGK'},
+        {'germline_set_name': 'IGLambda_VJ', 'locus': 'IGL',
+         'germline_set_id': '9606.IGL'}]}
+    LATEST = {'GermlineSet': [
+        {'release_version': 2.0, 'release_date': '2024-06-01T00:00:00'}]}
+
+    def get(self, url):
+        if '/latest' in url:
+            return Canned(self.LATEST)
+        if '/germline/sets/' in url:
+            return Canned(self.SETS)
+        if '/germline/species' in url:
+            return Canned(self.SPECIES)
+        raise AssertionError('unexpected url %s' % url)
+
+
+class TestHelpers(unittest.TestCase):
+    """
+    Tests for the small pure helpers
+    """
+
+    def test_normalize_version_strips_trailing_zero(self):
+        """An integer release reported as 3.0 becomes 3 for the URL."""
+        self.assertEqual(normalizeVersion(3.0), '3')
+        self.assertEqual(normalizeVersion('2.1'), '2.1')
+
+    def test_safe_set_name_replaces_separators(self):
+        """A set name with spaces and slashes becomes one safe token."""
+        self.assertEqual(safeSetName('C57BL/6J IGKV'), 'C57BL_6J_IGKV')
+
+
+class TestBucketChain(unittest.TestCase):
+    """
+    Tests for classifying an allele into a reference chain
+    """
+
+    def test_v_and_j_use_four_character_chain(self):
+        """V and J file under their own four-character chain."""
+        self.assertEqual(bucketChain('IGKV1-12*01', 'A' * 300),
+                         ('IGKV', KIND_VDJ))
+        self.assertEqual(bucketChain('IGKJ1*01', 'A' * 38), ('IGKJ', KIND_VDJ))
+
+    def test_short_ighd_is_the_diversity_segment(self):
+        """A short IGHD is the D segment and files under vdj."""
+        self.assertEqual(bucketChain('IGHD1-1*01', 'A' * 17), ('IGHD', KIND_VDJ))
+
+    def test_long_ighd_is_the_delta_constant(self):
+        """A long IGHD is the delta constant and files under the locus constant."""
+        self.assertEqual(bucketChain('IGHD*01', 'A' * 400), ('IGHC', KIND_CONSTANT))
+
+    def test_isotype_constant_files_under_locus_constant(self):
+        """A heavy isotype such as IGHM files under IGHC."""
+        self.assertEqual(bucketChain('IGHM*01', 'A' * 400), ('IGHC', KIND_CONSTANT))
+
+
+class TestSearchUnits(unittest.TestCase):
+    """
+    Tests for resolving a query to downloads
+    """
+
+    def test_emits_two_forms_per_set_with_human_ex_endpoint(self):
+        """Each set is fetched ungapped and gapped, human via the _ex endpoint."""
+        source = OgrdbSource(client=StubClient())
+        units = source.searchUnits(
+            Query(collection='human', filters={'locus': 'IGK'}))
+
+        self.assertEqual(len(units), 2)
+        formats = {u.metadata['format'] for u in units}
+        self.assertEqual(formats, {'ungapped', 'gapped'})
+        for unit in units:
+            self.assertTrue(unit.url.endswith('_ex'))
+            self.assertIn('/9606.IGK/2/', unit.url)
+            self.assertEqual(unit.metadata['version'], '2')
+
+    def test_locus_filter_narrows_to_one_locus(self):
+        """A locus filter fetches only that locus's sets."""
+        source = OgrdbSource(client=StubClient())
+        units = source.searchUnits(
+            Query(collection='human', filters={'locus': 'IGK'}))
+        self.assertEqual({u.metadata['locus'] for u in units}, {'IGK'})
+
+    def test_wildcard_covers_every_configured_locus(self):
+        """With no locus filter, every immunoglobulin locus is fetched."""
+        source = OgrdbSource(client=StubClient())
+        units = source.searchUnits(
+            Query(collection='human', filters={'locus': '*'}))
+        self.assertEqual({u.metadata['locus'] for u in units},
+                         {'IGH', 'IGK', 'IGL'})
+
+
+class TestBuildReference(unittest.TestCase):
+    """
+    Tests for splitting downloaded sets into per-chain FASTAs
+    """
+
+    def _entries(self, tmp):
+        entries = []
+        for fmt, fixture in (('ungapped', 'ogrdb_igk_ungapped.fasta'),
+                             ('gapped', 'ogrdb_igk_gapped.fasta')):
+            path = tmp / ('%s.fasta' % fmt)
+            path.write_text(readFixture(fixture))
+            unit = DataUnit(
+                unit_id='IGKappa_VJ.%s.fasta' % fmt, collection='human', url='x',
+                metadata={'species': 'human', 'locus': 'IGK',
+                          'set_name': 'IGKappa_VJ', 'format': fmt,
+                          'chains': ['IGKV', 'IGKJ']})
+            entries.append((unit, path))
+        return entries
+
+    def test_v_from_gapped_and_j_from_ungapped(self):
+        """V keeps its gaps from the gapped form; J comes from the ungapped."""
+        with tempfile.TemporaryDirectory() as tmp:
+            tmp = Path(tmp)
+            source = OgrdbSource(client=None)
+            source.buildReference(self._entries(tmp), tmp / 'reference_base')
+
+            vdj = tmp / 'reference_base' / 'human' / 'vdj'
+            v = (vdj / 'airrc_human_IGKV.fasta').read_text()
+            j = (vdj / 'airrc_human_IGKJ.fasta').read_text()
+
+            self.assertIn('IGKV1-12*01', v)
+            self.assertIn('.', v)                 # gapped V keeps its IMGT gaps
+            self.assertIn('IGKJ1*01', j)
+            self.assertNotIn('.', j)              # J is taken ungapped
+
+
+class TestAlias(unittest.TestCase):
+    """
+    Tests for the airrc alias
+    """
+
+    def test_airrc_resolves_to_ogrdb(self):
+        """'airrc' is an alias that resolves to the ogrdb source."""
+        from sourcerer.Sources import canonicalName, getSource
+
+        self.assertEqual(canonicalName('airrc'), 'ogrdb')
+        self.assertIsInstance(getSource('airrc', client=None), OgrdbSource)
+
+    def test_ogrdb_declares_the_alias(self):
+        """The source lists airrc among its aliases."""
+        self.assertIn('airrc', OgrdbSource.aliases)
+
+
+if __name__ == '__main__':
+    unittest.main()