Genostruct — integrated genome · gene · protein · structure viewer

A pipeline and self-contained web viewer that link protein 3D structure, protein sequence, and the underlying genome through a shared, per-element color code. For each protein that has a structure model, every secondary-structure element (α-helix / β-strand) is detected, assigned a unique color, and that same color is used to paint:

1. the residues in the 3D structure (Mol*),
2. the residues in the protein sequence track, and
3. the exact genomic bases (CDS codons, split across exons) that encode the element, in a purpose-built genome track.

The result: you can see a helix in 3D, find the same-colored stretch in the sequence, and read off the exon(s) and genomic coordinates that encode it — with introns and strand correctly represented.

---

Contents

File	What it is
pipeline/genostruct.py	The data pipeline (generic; no organism/file names hardcoded).
pipeline/build_viewer.py	Assembles the self-contained viewer HTML from pipeline output.
pipeline/genostruct_export.py	GenBank export of any transcript locus or genomic region, annotated with colored structural features.
viewer_src/	Front-end source: app.js, app.css, template.html.
genostruct_viewer.html	The deliverable — one self-contained file, built by the pipeline from your data.
tests/validate_designs.py	End-to-end correctness checks for the design pipeline.

No datasets are bundled

This repository contains the tool only. Pipeline output — viewer/data/, viewer/structures/, the built genostruct_viewer.html, and anything under exports/ — is generated from your own genome, annotation and structures, and is deliberately not committed. Open datasets will be published separately.

---

Quick start

Build it from your data, then open genostruct_viewer.html in any modern browser. No server, no network — Mol*, the decompressor, all structure coordinates and per-transcript data are embedded (gzip + base64; ~15 MB for 188 transcripts with tag-design annotation).

Pick a transcript from the dropdown. The four stacked panels are color-linked: hovering or clicking an element in any panel highlights it everywhere.

Viewer controls

• 3D structure (Mol*) — rotate (drag), zoom (scroll), pan (right-drag). Native Mol* interaction.
  • Insertion residue — selecting a design paints the residue(s) the tag exon splices between in magenta, in the 3D structure and in the sequence track, in every color mode and never dimmed.
  • Color mode — the button cycles the cartoon through three schemes: per-element colors (matched to the sequence and genome tracks), the AF2 pLDDT confidence ramp, and relative solvent accessibility (buried → exposed), which is what decides whether a tag-insertion junction is on the protein surface.
• Protein sequence — zoomable: use the − / reset / + buttons, or Ctrl/⌘ + scroll over the sequence.
  • Drag-select a residue range — dims the rest of the structure, recenters the 3D camera on the selection, and zooms the genome track to the encoding bases. Double-click the sequence to clear.
• Genome track — scroll to zoom (centered on the cursor), drag to pan, double-click to reset to the full locus. A readout shows current zoom factor and window size in bp.
  • Show genome 5'→3' / Show protein N→C — toggle the genome-track orientation. By default the track is drawn in the protein's N→C direction (reversed for minus-strand genes); the toggle switches to native 5'→3'.
• Structural elements (legend or sequence) — click an element (a legend chip, a colored sequence block, or a genome feature) to focus it in 3D (that element stays vivid, the rest dims), recenter the camera on it, and zoom the genome track to its bases. Click again or click empty space to clear. Hovering an element highlights it across all three panels.
• Tag knock-in candidates — a ranked table of insertion sites for this transcript, with a separate uniqueness column per nuclease. Selecting a design shows each guide aligned against its off-target sites: sequences in a fixed-width column with mismatched bases in red and the PAM in amber, and the genomic coordinates, mismatch count and score in a separate left-justified column beside them. Click a row to draw its insertion window, both guide pairs and their cut sites in the genome track, zoom to the intron, and highlight the junction residues in the sequence and 3D. The detail block lists each guide's spacer, PAM, GC, heuristic score, identical-site count, distance to the nearest other genomic site and where those sites are, plus the designed cassette sequence. Hide design lane removes the lane from the genome track.
• Search box — type a gene number to filter the transcript list; the full numerically-ordered dropdown stays available.
• Export GenBank (transcript / viewed region) — downloads an annotated GenBank record; the viewed-region export includes the full sequence of the current genome window (see note below).

---

Demo

A short screen recording of the viewer in use — selecting a transcript, focusing structural elements in 3D, and the linked highlighting across the structure, sequence, and genome panels:

Use_genostruct.mp4

(The video plays inline on GitHub. A copy is also committed at docs/Use_genostruct.mp4 for offline/download use.)

---

Running the pipeline on your own data

The pipeline is generic — nothing about the organism, file names, or gene IDs is hardcoded. It needs three inputs:

• a genome FASTA (*.fa),
• a GFF3 annotation (genes / mRNA / exon / CDS with phase),
• a folder of PDB structure models (e.g. AlphaFold predictions), for any subset of the proteins.

The protein FASTA is derived automatically from the genome + GFF3 (CDS translation, keyed by transcript ID, using the same pure-Python GFF3 parser/translator the pipeline validates its CDS round-trip against — no external tools required) — you no longer need to supply one. If you already have a peptide FASTA you trust (e.g. a curated annotation release), you can still pass it explicitly to skip derivation.

import sys; sys.path.insert(0, "pipeline")
import genostruct

genostruct.build(
    genome_fa   = "path/to/genome.fa",
    gff3        = "path/to/annotation.gff3",
    pdb_dir     = "path/to/pdbs",
    out_dir     = "viewer",
    # pep_fa is optional — omit it to derive protein.fa from genome+gff automatically.
    # pep_fa    = "path/to/pep.fa",
    # id_regex is optional. By default, each PDB is matched to a transcript by
    # SEQUENCE (searching the model's own residues against every candidate
    # protein), not by filename — protein/gene names in a PDB filename
    # frequently don't correspond to the transcript/gene IDs in a GFF3. Only
    # set id_regex if you already know filenames reliably encode the
    # transcript id and want to skip sequence search for speed:
    # id_regex  = r"(.+)_ranked_0",   # example for *_ranked_0.pdb
)

Or from the CLI:

python pipeline/genostruct.py \
  --genome path/to/genome.fa --gff path/to/annotation.gff3 \
  --pdb-dir path/to/pdbs --out viewer

This writes viewer/data/<transcript>.json, viewer/data/index.json, and copies matched PDBs into viewer/structures/.

Then build the standalone HTML:

python pipeline/build_viewer.py \
  --out viewer \
  --molstar molstar.js --molstar-css molstar.css --pako pako_inflate.min.js \
  --template viewer_src/template.html \
  --app-js viewer_src/app.js --app-css viewer_src/app.css \
  --html viewer/genostruct_viewer.html

(molstar.js, molstar.css, pako_inflate.min.js are bundled here; they are the only third-party assets and are embedded into the output.)

---

Web app (upload your own data, no local setup)

app.py is a Shiny for Python app that wraps the pipeline end to end: upload a genome FASTA, a GFF3, and a folder of PDB models, click Build viewer, and the app runs genostruct.build() (deriving protein.fa automatically, no external tools) + build_viewer.build_html() server-side, then renders the resulting self-contained viewer inline (with a download button for the generated HTML).

Run locally:

pip install -r requirements.txt
shiny run app.py

Deploy to shinyapps.io:

pip install rsconnect-python
rsconnect deploy shiny . --name <account> --title genostruct-viewer

Everything the app needs (shiny, numpy, biotite) is pip-installable — there's no system binary dependency to configure on the deploy target.

---

GenBank export

genostruct_export.py produces a GenBank record — with genomic sequence — for either a single transcript's locus or an arbitrary region. Structural elements are written as misc_features whose location is the exact codon ranges (as join(...), complement(...) for minus strand), carrying:

• /label = transcript:element (e.g. g113516.t1-00001:H2),
• /note = element kind + protein residue range,
• /color = the viewer hex color (e.g. #e6194B),
• /colour = Artemis/Geneious-style R G B (for tools that read that convention),
• /structure_type = helix or strand.

# a single transcript's locus, with 200 bp flanks
python pipeline/genostruct_export.py \
  --genome genome.fa --data-dir viewer/data \
  --transcript g113516.t1-00001 --flank 200 \
  --species "Acyrthosiphon pisum" --assembly JIC1_v1.0 \
  --out exports/g113516_locus.gb

# an arbitrary region (all overlapping transcripts included)
python pipeline/genostruct_export.py \
  --genome genome.fa --data-dir viewer/data \
  --region scaffold_2:89950000-90000000 \
  --out exports/region.gb

The Export GenBank button in the viewer produces a lightweight feature-table-only version (no sequence — kept out of the browser to hold file size down). For a complete record with sequence, use genostruct_export.py.

---

Tag knock-in design (ALFA / HiBiT via an artificial exon)

Given the same inputs the viewer uses, the pipeline will nominate introns where an artificial exon encoding an epitope tag can be knocked in by NHEJ, design the cassette, and find the overlapping gRNA pairs needed to cut there — one pair for Cas9 and, independently, one for Cas12a. Everything is annotated in the exported GenBank.

The premise: an artificial exon spliced into an internal intron places the tag at the corresponding exon–exon junction of the protein, which is tolerated when that junction sits in a solvent-exposed, unstructured loop on the outside of the fold.

# rank candidate introns, design cassettes, find overlapping gRNA pairs
python pipeline/tagdesign.py \
  --data-dir viewer/data --structures viewer/structures \
  --tags ALFA,HiBiT --out-tsv exports/tag_designs.tsv \
  --out-json-dir viewer/data_tagged

# then export an annotated GenBank for any transcript
python pipeline/genostruct_export.py \
  --genome genome.fa --data-dir viewer/data_tagged \
  --transcript g119212.t1-00001 --flank 200 --max-designs 2 \
  --out exports/g119212_tag_design.gb

Module	What it does
pipeline/introns.py	Coding introns in transcript orientation; splice dinucleotides; intron phase (from cumulative coding length, cross-checked against the GFF3 phase column); the NHEJ insertion-safe window.
pipeline/surface.py	Per-residue SASA / relative solvent accessibility (Shrake–Rupley in numpy — no new dependency).
pipeline/tagsites.py	Scores each intron's exon–exon junction on loop flanks (below), exposure, clearance from secondary structure and distance from the termini.
pipeline/cassette.py	Designs the artificial exon and validates it by in-silico splicing.
pipeline/grna.py	Cas9 / Cas12a site finding on both strands, and overlapping-pair detection.
pipeline/tagdesign.py	Drives all of the above; writes ranked TSV + augmented JSON.
tests/validate_designs.py	Independent end-to-end correctness checks.

The tag must land inside a loop — both flanks

The artificial exon is spliced between the last residue the upstream exon encodes and the first residue the downstream exon encodes. For the tag to sit inside a surface loop, both of those must be coil, and stay coil for a few residues.

Averaging coil content over a window around the junction is not sufficient — it passes sites where one flank is a long loop and the other is a helix, so the tag lands at the edge of the loop, hard against secondary structure:

          upstream exon | downstream exon          verdict
  I3      aaaacccc | cccaaaaa     4/4  loop both sides      accepted
  I8      accccccc | caaaaaaa     7/2  end of a loop        rejected
  I9      aaaaaaac | cccaaaaa     1/4  start of a loop      rejected
  I14     cccccccc | aaaaaaaa    19/1  loop then helix      rejected

So tagsites.flank_coil_runs() counts the consecutive coil residues running back into the upstream exon and forward into the downstream exon separately. Each run is confined twice over:

• to the residues its own exon encodes, so a loop in a neighbouring exon cannot be borrowed across an intron;
• to the residues the structure model actually covers. residue_sse defaults to coil wherever the model does not reach (a trimmed signal peptide, say), so counting those would manufacture a long "loop" out of absent evidence — precisely where a tag could land in structure nobody looked at.

Both runs must reach --min-flank-coil (default 3); --allow-nonloop-flanks relaxes it.

This is a hard gate, and a strict one: on the bundled dataset it rejects 1724 of 2247 otherwise-viable junctions, leaving 490 designs whose tags are genuinely interior to a loop that was actually modelled. Both runs are reported in the TSV (upstream_coil_run, downstream_coil_run), in the viewer's design table, and drive the unstructured subscore.

The insertion-safe window

A cassette is spliced in using its own acceptor and donor, while the native donor and acceptor keep serving the flanking exons. So the break must avoid the native signals: by default the first 30 bp of the intron (5′ donor) and the last 60 bp (branch point, polypyrimidine tract, 3′ acceptor) are excluded (--donor-pad, --acceptor-pad). Guides are only reported when their whole footprint lies inside that window.

Frame handling — every phase is usable

With intron phase p (nucleotides of the junction codon contributed by the upstream exon), the exon body is pad5 + tag codons + pad3 where pad5 = (3-p) % 3 and pad3 = p % 3, so the body is always a multiple of 3 and the downstream exon stays in frame. For a split junction (p ≠ 0) the padding reconstitutes the junction codon from the flanking exons' own bases, so both resulting codons are guaranteed to be sense codons and the native residue is duplicated rather than mutated. Phase is not used to rank sites.

Introns at or past the stop codon are refused — a cassette spliced in there is never translated.

Nothing is organism-specific

Both species-dependent inputs are learned from the supplied annotation:

• codon usage, counted over every CDS, used to encode the tag;
• splice-site consensus, as 5′/3′ position weight matrices over every annotated intron, used to build the synthetic acceptor and donor.

On the bundled A. pisum data that recovers AAG|GTAAGTA (1.96 bits on the donor GT), a long T-rich polypyrimidine tract before AG (1.99 bits), and strongly AT-biased codon usage — all correct for this genome, with no species assumption anywhere in the code.

Overlapping gRNA pairs

Cutting is much more efficient when two guides target overlapping sites, so the unit of output is a pair whose protospacer+PAM footprints overlap.

Nuclease	PAM	Spacer	Cut
SpCas9	NGG, 3′ of the spacer	20 nt	blunt, 3 bp 5′ of the PAM
Cas12a (As/Lb)	TTTV, 5′ of the spacer	23 nt	staggered, 5-nt 5′ overhang

Pairs are found independently per nuclease, and either nuclease is sufficient — a site with an overlapping Cas9 pair but no Cas12a pair (or vice versa) is a perfectly usable design and is kept and scored on its own merits, not penalised. A design is only rejected when neither nuclease has an overlapping pair, since then there is no way to make the break. nucleases_with_pairs records which apply; --require-both-nucleases opts in to the stricter rule.

Tunable: --min-overlap, --max-cut-separation, --opposite-strands, --min-efficiency.

On-target efficiency is a labelled heuristic. No published trained scorer (Doench Rule Set 2/3, Azimuth, DeepCpf1) is installable as a pure-Python offline dependency — they ship as R/Bioconductor wrappers over Keras — and this project deploys without system binaries. The default scorer therefore combines only well-established, checkable determinants (GC optimum, poly-T terminator, homopolymers, hairpin propensity, seed composition, PAM-proximal base) and reports every term separately. It is a relative ranking aid, not a calibrated activity prediction, and is labelled heuristic_v1 wherever it appears. Register a real trained model with grna.register_scorer() to replace it.

Why pLDDT is reported, not scored

On this dataset corr(window pLDDT, junction RSA) = −0.33: AlphaFold draws low-confidence regions as extended chain, which inflates apparent solvent exposure. Rewarding "flexible" (low pLDDT) and "exposed" together would count one artefact twice and flood the ranking with the least trustworthy geometry.

Low pLDDT is also genuinely ambiguous — a good predictor of intrinsic disorder (a real positive for an internal tag) and a warning that the RSA and secondary-structure calls there are unreliable. So pLDDT becomes a reported site_class (ordered_loop / intermediate / disordered) that you filter on with --site-class, rather than a hidden term in the score.

What gets annotated in the GenBank

Feature	Meaning	Colour
intron	candidate intron, with phase, length, splice class, site score	grey
misc_feature :insertion_window	the NHEJ-safe interval	pale green
misc_feature :insertion_point	preferred insertion coordinate; carries the full cassette sequence in /cassette	dark green
misc_feature :junction_resN	the junction codon(s), as join() of their actual bases	purple
protein_bind	each guide's protospacer, on its own strand, with sequence, PAM, GC, efficiency and pair overlap	blue (Cas9) / orange (Cas12a)
misc_feature _PAM / _cut	the PAM, and the blunt cut or staggered overhang	light / red

Off-target assessment

pipeline/offtarget.py compares the whole protospacer plus its PAM against every PAM-adjacent site in the assembly. The search is exhaustive up to a mismatch budget, by enumeration rather than seeding: every sequence within --max-mismatches of a guide is generated (32,551 variants for a 20-mer at ≤3), 2-bit encoded and sorted once, then the genome is streamed and each site matched with a vectorised searchsorted. Guides are batched so the variant index fits in memory; the bundled 526 Mb assembly takes a couple of minutes, cached by genome size, mtime and format version.

This matters: a seed-indexed search (match the 12-nt PAM-proximal seed exactly, then count mismatches) cannot see an off-target that differs inside the seed. Switching to the exhaustive search on this dataset cut the guides called unique from ~55% to ~25% — the earlier figure was an artefact of an incomplete search, not a property of the guides.

Position-weighted scoring

A mismatch next to the PAM largely abolishes cutting; one at the distal end is often tolerated. So each hit carries a weighted score rather than a bare mismatch count:

w(position) = W_MIN + (W_MAX - W_MIN) * proximity_to_PAM ** W_EXP
hit_score   = PROD (1 - w_i)   over mismatched positions i

A single mismatch beside the PAM gives 0.05 (very unlikely to cut); one at the far end gives 0.80 (likely still cuts). Correctly inverted for Cas12a, whose PAM is 5′. Per guide, specificity = 1 / (1 + Σ hit_score) over all other sites, so 1.0 means nothing else was found.

Heuristic, not a fitted model. This has the same shape as published specificity models (CFD/MIT) but not their fitted values — no published coefficient table is reproduced here from memory. Install a measured table with offtarget.set_position_weights(nuclease, weights) to replace it.

Each guide reports:

field	meaning
additional_identical_sites	other places the protospacer occurs exactly with a valid PAM. 0 for a unique guide — the guide's own locus is not counted
found_own_locus	whether the guide was located at all (false means the spacer is not in this assembly)
nearest_mismatches	mismatches to the closest other site, or none if nothing within the budget
specificity	1 / (1 + Σ hit_score) over every other site; 1.0 = clean
worst_hit_score	the most plausible single off-target
total_other_sites	how many other sites exist within the budget
hits / closest	each site's sequence and PAM, scaffold:start-end(strand, Nmm), the mismatched positions and that site's hit_score, most dangerous first

These appear in the TSV, in the viewer (genome-track tooltips and the design table) and in the GenBank as /additional_identical_sites, /nearest_other_site and /closest_offtargets.

Scope. Sites differing by more than --max-mismatches (default 3) are not reported. nearest_mismatches = none means nothing within that budget, not genome-wide uniqueness. Raising to 4 costs ~13× more variants per guide (425k) and is best run on a shortlist.

Why so few guides are unique here

Not a scoring artefact — the BICYCLE family is paralogous:

• Only 944 of 1762 candidate spacers have no additional identical site. The rest occur perfectly elsewhere — some in repeats, with up to 7205 extra copies. That is fatal regardless of weighting.
• Uniqueness is reported per nuclease, because a site can be unusable with Cas9 and clean with Cas12a. 118 of 490 designs have a fully unique pair for at least one nuclease (Cas12a 96, Cas9 43), against only 31 where every nuclease present is clean — so a single combined flag would have hidden 87 usable designs.

Two different questions are reported separately rather than collapsed into one "unique" flag, because they mean different things:

shown	meaning
additional identical sites	the protospacer occurs again exactly. Fatal.
worst off-target score	no exact duplicate, but the best mismatched site is plausible enough to cut

In the viewer each nuclease column shows duplicates · worst score, green when neither applies, amber when there is no duplicate but a plausible mismatch site, red when an exact duplicate exists.

Use --require-unique to keep only those, and sort on specificity rather than treating uniqueness as binary.

Validation

tests/validate_designs.py re-derives every claim independently of the code that made it. On the bundled dataset, all pass:

• 2396 introns — cumulative-length phase == GFF3 phase column; windows strictly inside their intron
• 2393 flank coil runs — recounted straight off the secondary-structure string, confined to the correct exon and to the modelled region
• 2393 cassettes — frame preserved, padding matches phase, native protein preserved on both sides, no internal stop, both tags in frame
• 147613 gRNA sites — spacer and PAM re-extracted from their genomic coordinates, PAM motif valid, footprint size, overhang length, PAM adjacency
• 1787 pairs — reported overlap exact, guides distinct, cuts clear of the native splice signals
• GenBank round-trip — each annotated guide feature extracts back to its stated spacer
• off-target search — identical-site, per-mismatch and total counts all recompared against brute force on a miniature genome; every located hit re-extracted to confirm it differs at exactly the stated positions; scores bounded; PAM-proximal mismatches verified to score below distal ones

python3 tests/validate_designs.py --data-dir viewer/data \
  --genbank exports/g119212_tag_design.gb

---

Method notes

• Secondary structure is computed from the model's Cα coordinates with biotite's annotate_sse (the P-SEA algorithm), then segmented into contiguous helix/strand elements. Each element gets a unique color from a warm palette (helices) or cool palette (strands), cycled if there are many.
• Residue → genome mapping walks each CDS block with its GFF phase, in translation order, on the correct strand, so every protein residue maps to its three genomic codon positions — correctly handling introns and codons that straddle exon boundaries.
• Model ↔ protein alignment. AlphaFold models here were predicted from signal-peptide-trimmed sequences and are numbered from 1, so each model is aligned to its full protein to recover an N-terminal offset (exact substring → gapless scan → local alignment fallback). Colors are indexed by full-protein residue, and the Mol* color theme applies the offset so 3D coloring lines up with the sequence and genome.

---

Validation

• Coordinate correctness: for all 187 transcripts, the exported CDS translates back to the exact deposited protein (187/187), and extracting + translating every structural-element misc_feature reproduces its protein residues exactly (1113/1113 elements, 187/187 transcripts fully correct — no CDS mismatches, no out-of-range elements, no empty ranges).
• Element residue ranges are clamped to the coding protein length in the pipeline, so models that are longer than their protein (from a gapped local alignment) can no longer map an element past the C-terminus into the stop codon or out-of-range codon positions.
• validation_colorlink.png visually confirms sequence-color → genome-segment correspondence, including multi-exon elements on the minus strand.

Known limitations

• 3 PDB models in the source set correspond to stale isoform IDs absent from the protein FASTA and are excluded (documented in the pipeline output warnings).
• 10 models differ from their deposited protein by internal indels and are placed by local alignment (identity ≥ 0.6); their per-residue mapping is still validated correct by the CDS/element round-trip above.
• The viewer runs entirely offline in a single self-contained HTML file. If it is opened inside a sandboxed iframe with a restrictive Content-Security-Policy (some hosted environments do this), the file-download of a GenBank export may be blocked by the sandbox; a copy-to-clipboard / open-in-new-tab fallback is provided for that case. Opened directly in a browser (double-clicking the file, or served from a normal web server), downloads work normally.