Genomics and DNA engineering: parameters and notes¶
See the tutorial for purpose, methods and examples. These APIs show progress by default; use verbose=False to disable it.
Genomic relationships and breeding prediction¶
Install numeric dependencies:
Rows are individuals; columns must have identical marker order and allele orientation. Only diploid numeric dosages in [0, 2] are supported. None denotes missing dosage; missing="error" rejects missingness. Training means control imputation and centering, including prediction rows. All-missing and allele-monomorphic training markers are excluded and reported.
The relationship matrix is VanRaden method 1's additive relationship matrix; conventional kinship is G/2. Prediction is fixed-lambda ridge regression/RR-BLUP: ridge_lambda is required, and 1.0 is only an example. Select it using training-only validation and account for relatedness when splitting samples. This interface does not fit REML, polyploid or multi-trait models, environmental covariates, or mating allocation.
Prokaryotic gene prediction¶
Install the real Pyrodigal backend:
The tutorial's artificial example demonstrates the interface; it does not guarantee a gene call. Supply one bacterial/archaeal linear contig per call. meta=True uses upstream pretrained models. meta=False requires an explicit suitable training_sequence; failures do not silently switch modes. translation_table is restricted to training mode.
Outputs include CDS, proteins, model scores, translation tables and backend version. Coordinates are 0-based half-open on the input; negative-strand CDSs are reverse-complemented. Ambiguous IUPAC symbols become N without shifting coordinates. The function does not identify eukaryotic intron/exon structures or assign biological gene functions. Circular input requires deliberate linearization.
Methylation counts¶
The reader accepts modkit's 18-column extended bedMethyl, including gzip/BGZF. Column 11 is a percentage in [0, 100]. Summaries preserve chromosome, strand, modification and motif groups; the fraction is total modified count divided by total valid coverage. Zero coverage produces no fraction, and 5mC/5hmC remain separate. Defaults are max_records=1_000_000, max_bytes=268_435_456 decompressed bytes and max_line_bytes=65_536. Exceeding a limit raises an error without returning partial data; partition larger inputs or select suitable limits.
This performs descriptive summaries, not differential methylation tests or disease prediction. The external modkit-pileup adapter produces bedMethyl from modified BAM files.
DNA storage¶
Install the codec:
This adapts Chamaeleo's fixed-Huffman Goldman codec. Preserve the entire JSON packet: indices, lengths and SHA-256 are JSON metadata, not DNA-encoded. Intact reordered segments can be restored; missing/corrupted segments or digest mismatches fail. The digest detects corruption but provides neither correction nor authentication. This does not implement the original paper's complete overlapping archive protocol or any physical synthesis/sequencing experiment.
DNA nanostructures and oxDNA export¶
Install the geometry dependency:
The design is scadnano JSON suitable for a .sc file. Export returns configuration and topology strings, without writing files or running simulation. Caller-supplied designs are also accepted, with complete canonical DNA assignments on every strand.
Only linear strands are supported: the verified scadnano 0.21.1 export loses circular closure bonds, so circular designs are rejected. Non-finite geometry or unsupported upstream offsets return explicit errors. Geometry construction/export does not establish folding, relaxation or molecular dynamics results.
Genotype PCA¶
genotype_pca() uses NumPy, installed through .[genomics] or .[bio-analysis]. Input rows are samples and columns are aligned diploid biallelic dosages in [0, 2], including fractional dosages. None and NaN denote missing values. Optional sample_ids must be unique, nonempty and match the row count.
Estimate each observed marker mean μ and frequency p=μ/2, remove all-missing or observation-invariant markers, impute other missing entries with μ, and calculate an SVD of (dosage−μ)/sqrt(2p(1−p)). Scores are U*S, retained-marker loadings are V, and covariance eigenvalues are S²/(n_samples−1). Variance ratios use the full spectrum. components must be in 1..min(n_samples−1, n_markers_used).
marker_indices maps loading rows to original columns. Means and frequencies retain the full input width, with null for all-missing markers. Imputation counts cover retained markers only. The largest absolute loading is positive; degenerate eigenspaces may still rotate between numerical libraries. Input arrays are not mutated and output is JSON-compatible.
This dense in-memory PCA does not read BED/VCF, prune LD or exclude related samples. Its eigenvalue scale differs from PLINK's relationship-matrix scale. Harmonize alleles and perform QC first; large cohorts should use specialized backends.
Genome ANI¶
average_nucleotide_identity(query, reference) uses PyFastANI from .[genomics], without a separate FastANI executable. Supply DNA strings or local FASTA paths, including gzip. Strings containing only IUPAC letters are DNA; use ./filename or a Path object for such filenames. All FASTA records describe contigs of ONE genome, with boundaries preserved, rather than separate samples.
The interface retains FastANI defaults, including 3,000 bp query fragments, for related prokaryotic genomes. Short contigs and incomplete terminal fragments cannot contribute complete query fragments. Inputs are loaded into memory; ambiguous bases become N and their count is recorded.
Only query→reference is computed; swap inputs for the reverse direction. ani_percent measures identity in homologous fragments. query_aligned_fraction is matching fragments divided by complete query fragments, not base-level coverage. A missing report returns status="no_match" with null identity and aligned fraction, never a zero-percent identity claim.
Results record versions, parameters, input digests, fragment counts and method citations. Progress shows loading, indexing and querying stages; verbose=False disables progress but does not suppress upstream warnings such as short-reference warnings. Precompiled wheel availability depends on the dependency and platform.
Full parameters¶
dnakit.genomics
¶
Optional genomic prediction, gene finding, methylation and DNA engineering APIs.
genotype_pca(genotypes, *, components=2, sample_ids=None, verbose=True)
¶
Compute standard PCA on frequency-standardized diploid allele dosages.
Rows are samples and columns are aligned biallelic markers. Dosages may be fractional and must be in [0, 2]. None/NaN means missing; strings, booleans, complex values and infinity are rejected. Numeric NumPy matrices are also accepted. At least two samples and one variable marker are required.
For each marker, estimate mean dosage mu using observed samples and set p = mu/2. Drop entirely missing or observation-invariant columns (including all-heterozygous columns); fill remaining missing entries with mu. Calculate X = (dosages - mu) / sqrt(2p(1-p)), then a thin SVD X = USV.T. Scores are US, marker loadings are V, and covariance eigenvalues are S*2/(n_samples-1). Variance ratios use the full spectrum as denominator, even when only a subset of components is requested.
The largest absolute loading in each component is assigned a positive
sign (ties select the first retained marker). Equal eigenvalues may still
rotate their non-unique subspace across numerical libraries. Components
must be in 1..min(n_samples-1, n_markers_used); collinearity can produce
zero-variance components. numerical_rank reports the nonzero rank.
Results are JSON-compatible. Scores have shape samples x components;
loadings have shape retained-markers x components, in marker_indices
order. Indices are zero-based. Means and frequencies cover ALL input
columns, with null for wholly missing markers. Imputation counts concern
retained markers only. Neither dosages nor sample labels are modified.
This dense, in-memory implementation does not read PLINK files, prune LD, filter related samples, or provide PLINK-compatible eigenvalue scaling. Its covariance eigenvalues differ from eigenvalues of X*X.T/n_markers. Large cohorts should use an appropriate specialized analysis backend.
average_nucleotide_identity(query, reference, *, threads=1, verbose=True)
¶
Estimate directional genome ANI with the optional PyFastANI library.
Install with pip install 'dnakit[genomics]'. Each input is a non-empty,
ungapped IUPAC DNA string or a local FASTA path (gzip is supported by DNAKit
I/O). A string consisting only of IUPAC letters is always DNA; use a
Path object or an explicit ./ prefix for such a filename. FASTA
text itself is not an input format. Every FASTA record is a contig of ONE
genome; contigs are never concatenated. Inputs are loaded into memory.
Ambiguous bases are represented as N without removing positions.
FastANI's standard settings are retained: 3,000-base query fragments, 16-mers, 80-percent identity tuning, and a 0.2 shared-genome reporting threshold. Contigs shorter than a query fragment and incomplete terminal fragments cannot contribute query fragments. Reference contigs shorter than the backend's sketching window are skipped by PyFastANI.
ANI is a percentage for the reported homologous fragments, not an identity
across every input base. query_aligned_fraction is matching fragments
divided by all complete query fragments, not base-level genome coverage.
The comparison runs query -> reference only and can be asymmetric; call
again with swapped inputs if both directions are needed. This API targets
related prokaryotic genomes, not short-read alignment or distant genomes.
status='no_match' means the backend reported no estimate (or no query
fragment exists); ANI, matching counts and aligned fraction are then None.
It does not mean zero-percent identity. Invalid or non-finite backend
results raise BackendExecutionError. Output is JSON serializable.
Progress reports loading, indexing and mapping stages, not base-level work.
genomic_relationship_matrix(genotypes, *, ploidy=2, missing='mean', verbose=True)
¶
Calculate the VanRaden method-1 additive genomic relationship matrix.
Rows are samples (at least two), columns are biallelic markers, and entries
are alternate-allele dosages in [0, 2]; fractional dosages are accepted.
NaN/None (JSON null) is mean-imputed per marker with missing='mean' or rejected with
missing='error'. Entirely missing and allele-monomorphic (frequency 0 or
1) columns are excluded and reported by zero-based indices. An all-
heterozygous marker is allele-polymorphic and retained, as in rrBLUP A.mat.
Text and booleans are not dosage measurements and are rejected.
With frequencies estimated from these samples, Z = M - 2p and G = Z Z' / (2 sum(p (1-p))). This is an additive relationship matrix; the coancestry/kinship matrix under the usual convention is G/2. It is relative to this sample's frequencies, not an inferred base population. No shrinkage, pedigree blending, or polyploid generalization is applied.
Returns a JSON-compatible dictionary including relationship_matrix,
full-width frequencies (None for entirely missing markers), and filtering
metadata. Uses O(samples**2) memory for the resulting dense matrix.
References: P. M. VanRaden (2008), doi:10.3168/jds.2007-0980; J. B. Endelman and J.-L. Jannink (2012), doi:10.1534/g3.112.004259.
genomic_prediction(train_genotypes, phenotypes, predict_genotypes, *, ridge_lambda, ploidy=2, missing='mean', verbose=True)
¶
Fit an intercept and fixed-lambda RR-BLUP marker effects, then predict.
Both genotype matrices use sample rows and the SAME marker/allele column order, with diploid alternate-allele dosages in [0, 2]. Arrays carry no marker IDs: callers must harmonize both marker identity and effect alleles. Training requires at least two samples and one finite phenotype per row. Missing phenotypes are rejected; filter them and their genotype rows first.
The objective is sum((y - mu - Z b)2) + ridge_lambda * sum(b2),
with unpenalized mu and Z centered using only training-marker means.
ridge_lambda must be finite and positive. It is the residual/marker
effect variance ratio for this dosage coding, not an estimated heritability.
Unlike R rrBLUP mixed.solve, this function does not estimate lambda by
REML. Choose it with training-only/nested validation outside this function.
Missing dosages are replaced by training means (or rejected), and entirely missing or allele-monomorphic training columns are excluded. Prediction samples never affect filtering, frequencies, centering, or model fitting. A thin SVD solves the ridge objective without forming a marker-square matrix.
Returns JSON-compatible predictions, training_predictions, full-width
marker_effects (zero for excluded markers), and filtering metadata.
intercept is the phenotype mean for CENTERED dosages; raw_intercept
applies to uncentered, imputed dosages. breeding_values are predictions
minus the training phenotype mean. This is an additive research baseline,
without environmental covariates, genotype-by-environment effects, or a
prediction-accuracy claim.
Reference: J. B. Endelman (2011), doi:10.3835/plantgenome2011.08.0024.
predict_prokaryotic_genes(sequence, *, sequence_id='sequence', meta=True, training_sequence=None, translation_table=None, closed=False, mask=True, include_stop=True, verbose=True)
¶
Predict prokaryotic CDSs using an optional, locally installed Pyrodigal.
Install with pip install 'dnakit[genomics-genes]'. Metagenomic mode
uses Prodigal's supplied models and selects the translation table from
the winning model; each gene reports that table. To train a genome model,
explicitly pass meta=False and training_sequence=.... Training
uses table 11 unless translation_table is provided. Pyrodigal checks
training length and rejects unsuitable short inputs; no silent fallback
to metagenomic mode is performed.
Input must be non-empty, ungapped IUPAC DNA. Ambiguous IUPAC letters are
represented as N without changing positions. Returned CDS strings are
oriented 5-prime to 3-prime on the predicted gene strand. Protein strings
use the gene's prediction table and include a terminal * when present
unless include_stop=False. mask=True uses Pyrodigal's default
50-base unknown-region masking threshold. closed controls contig
edge predictions, not circular topology.
Progress shows prediction/training stages; the backend does not expose nucleotide-level progress. No software, weights or sequence data are downloaded by this call. This predicts CDS structure, not gene function.
read_bedmethyl(path, *, max_records=1000000, max_bytes=268435456, max_line_bytes=65536, verbose=True)
¶
Read modkit's 18-column BED, detecting gzip/BGZF by its magic bytes.
All limits apply to decompressed input. Exceeding a limit raises, never silently truncates. Empty lines, # comments, track and browser headers are ignored; malformed records are rejected with a physical line number. Coordinates are 0-based, end-exclusive. Column 11 is percent (0..100); fraction_modified is calculated from counts and is None at zero coverage. bedRMod has a different denominator and is deliberately unsupported.
summarize_methylation(records, *, modification_codes=None, chromosomes=None, strands=None, min_coverage=0, max_records=1000000, verbose=True)
¶
Summarize validated records from read_bedmethyl using sum(mod)/sum(valid).
Groups always preserve chromosome, modification code, strand and motif. In particular, 5mC (m), 5hmC (h), and combined strands (.) remain separate. Duplicated/overlapping records are counted as supplied; record_count is not a count of unique CpGs. No cross-sample tests or disease inference is made. Empty filter sequences explicitly select nothing. Every input is validated before filtering, and undefined zero-coverage fractions remain None.
encode_dna_storage(data_hex, *, segment_bytes=32, verbose=True)
¶
Encode hex-encoded bytes as indexed DNA payloads with an integrity digest.
Segment indices and metadata are stored in the returned JSON packet, not
synthesized into the DNA. Preserve that packet for decoding. No sequencing
error correction, wet-lab storage, or physical indexing is performed.
Input is limited to one MB; install dnakit[dna-storage] explicitly.
decode_dna_storage(packet, *, verbose=True)
¶
Decode an intact Goldman packet, rejecting missing/corrupt segments.
A SHA-256 match checks accidental corruption; it is not authentication or an error-correcting code. DNA-only sequences without packet metadata are insufficient for this interface.
design_dna_duplex(sequence, *, verbose=True)
¶
Create a two-strand scadnano design with explicitly assigned complementary DNA.
Requires dnakit[dna-nanotechnology]. Return a scadnano JSON object that
can be saved as a .sc file or passed to export_oxdna_design. This simple
duplex constructor does not optimize origami folding or mechanical stability.
export_oxdna_design(design, *, verbose=True)
¶
Export caller-supplied scadnano geometry to oxDNA configuration and topology text.
No files are written and no oxDNA process is run. All strands must have complete canonical sequences; wildcard bases cannot become arbitrary bases. Only linear strands are supported: scadnano's exporter does not preserve circular bonds. Non-finite geometry and invalid backend output are rejected. Relaxation, simulation settings and physical interpretation are separate tasks.