基因组与 DNA 工程：参数与补充说明¶
作用、计算方法和示例见功能教程。以下接口默认显示进度，可使用 verbose=False 关闭。
基因组关系矩阵与育种预测¶
安装数值计算依赖：
关系矩阵使用 VanRaden 方法 1 的加性关系定义；常见亲缘系数矩阵为其二分之一。预测使用固定惩罚系数的 RR-BLUP／岭回归，ridge_lambda 必须显式提供；示例值不代表最优值。按育种周期、群体或家系划分验证集，避免亲缘泄漏。
缺失剂量可以使用 None，默认以训练样本均值填补；missing="error" 可要求完整输入。全缺失和等位基因单态训练标记会排除，返回其索引。预测样本不参与频率估计或参数拟合。当前只接受二倍体；不包含多倍体、多性状、环境协变量、REML 估计或亲本配组优化。
原核基因识别¶
安装真实 Pyrodigal 后端：
教程中的人工序列仅演示接口，不保证模型检出基因。真实分析每次传入一条细菌／古菌线性 contig。meta=True 使用上游预训练模型；meta=False 需显式提供适于训练的 training_sequence，短序列训练失败不会悄悄切换模式。translation_table 只用于训练模式。
结果包括 CDS、蛋白、模型评分、翻译表与后端版本。坐标为输入序列上的 0-based、左闭右开；负链 CDS 已反向互补。模糊 IUPAC 字母转为 N 并保留原位置。该接口识别原核 CDS，不预测真核外显子／内含子或赋予基因功能名称。环状 DNA 需要先明确线性化位置。
甲基化计数与汇总¶
接受 modkit 的 18 列扩展 bedMethyl，包括 gzip／BGZF。第 11 列使用百分数单位（0–100）。结果按染色体、链、修饰代码和 motif 分组，比例为修饰读数总和／有效覆盖读数总和；零覆盖返回空比例，不混合 5mC 与 5hmC。默认 max_records=1_000_000、max_bytes=268_435_456（解压后字节数）、max_line_bytes=65_536；超限会报错，不返回截断结果。较大数据应分区或显式调整限制。
该统计不执行差异甲基化检验或疾病预测。从带修饰信息的 BAM 生成 bedMethyl，可使用 外部分析 的 modkit-pileup。
DNA 数据存储编码与解码¶
安装编码器：
通过 Chamaeleo 的固定 Huffman Goldman 编码器转换字节。必须保存完整 JSON packet：索引、长度和 SHA-256 存在 JSON 中，未编码到 DNA 链内。支持完整片段乱序重排；丢段、损坏或摘要不符会报错。摘要用于检错，不提供纠错或认证。该接口不包含原论文完整的重叠冗余存档协议，也不执行合成或测序。
DNA 纳米结构与 oxDNA 导出¶
安装结构设计依赖：
design 是可保存为 .sc 文件的 scadnano JSON。导出结果包含 configuration、topology，可供后续 oxDNA 工作流保存使用。也可传入自行设计的 scadnano JSON；所有链必须完整分配 A/C/G/T。
当前仅导出线性链：验证的 scadnano 0.21.1 会丢失环状链闭合键，因此接口拒绝环状设计。无效或非有限几何参数、上游无法处理的偏移会返回明确错误。结构构造和文件导出不等于折叠优化、能量松弛或分子动力学模拟。
基因型 PCA¶
genotype_pca() 使用 NumPy，安装 .[genomics] 或较小的 .[bio-analysis] 即可。输入是样本×位点的二倍体双等位剂量矩阵，允许 0–2 的小数剂量，以 None 或 NaN 表示缺失。样本名称可通过 sample_ids 传入，必须非空、唯一且数量匹配。
先用有观测样本估计各位点均值 μ 和频率 p=μ/2，删除全缺失或观测值不变的位点，再以 μ 填补其余缺失值，按 (dosage−μ)/sqrt(2p(1−p)) 标准化并计算 SVD。返回样本得分 U*S、保留位点载荷 V 和协方差特征值 S²/(n_samples−1)；方差比例的分母包含全部主成分。components 必须在 1 到 min(n_samples−1, n_markers_used) 之间。
marker_indices 映射载荷行至原始位点；均值和频率保留原始列数，全缺失列为 null。填补计数只统计保留位点。最大绝对载荷取正以固定符号；相等特征值对应的子空间仍可能随数值库旋转。输入不被修改，结果可 JSON 序列化。
这是内存中的稠密矩阵 PCA；不读取 BED/VCF、不剪枝 LD，也不自动排除亲缘样本。它与 PLINK 的关系矩阵特征值采用不同尺度。应先统一位点／等位基因方向并完成 QC；大规模队列继续使用专业后端。
基因组 ANI¶
average_nucleotide_identity(query, reference) 使用 .[genomics] 中的 PyFastANI；无需单独安装 FastANI 可执行程序。它接受 DNA 字符串或本地 FASTA 路径（支持 gzip）。仅由 IUPAC 字母组成的字符串按序列解释；这种名称的文件使用 ./文件名 或 Path 对象。FASTA 内所有记录属于同一个基因组，contig 边界保持不变；不把多条记录当作不同样本。
默认使用 FastANI 的 3,000 bp 查询片段等参数，适合近缘原核基因组；短 contig 与不完整末端片段不能贡献完整查询片段。当前实现将序列读入内存；模糊碱基转换为 N 并记录数量。
只计算 query→reference 方向；需要双向结果时交换输入再调用。ani_percent 是同源片段的百分比一致性；query_aligned_fraction 为匹配片段数／完整查询片段数，不是全基因组碱基覆盖率。无可报告命中时返回 status="no_match"，ANI 和匹配比例为 null，不能解释为 0% 一致性。
结果记录后端版本、参数、输入摘要、片段数和方法引用。默认进度条显示加载、建索引和查询阶段；verbose=False 关闭进度，不屏蔽上游对短参考序列等情况的警告。依赖的预编译 wheel 是否可用取决于平台。
完整参数¶
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.