Source: https://docs.nygen.io/CyteArc/api/variants.html Revision: 905308f593408f3d36f3ad7e0202f487e1f34646 Variants API reference ********************** Mitochondrial lineage tracing (mtscATAC-seq, MAESTER) and targeted genotyping (GoT, GIFT, Tapestri) count, for each cell, the reads or UMIs of each allele at a set of sites: every position of the mitochondrial genome or of a targeted panel, or targets and probes without coordinates. Each tool writes those counts in its own files. "ds.variants" reads them into one "variant_counts" artifact that holds the sites and their alleles, and each cell’s depth at each site, its count of each allele and the calls the tool made itself. "ds.variants.import_counts" imports the counts of one file or folder, or several of one format, "ds.variants.load_counts" reads them as a "VariantCounts" record, and "ds.plots.variant_depth" draws the depth at each site and in each cell. "ds.variants.informative" scores every variant across cells and selects the informative ones, and "ds.variants.clones" calls the clones that those variants mark. For targeted sites, "ds.variants.genotypes" calls each cell’s genotype and "ds.variants.genotype_clones" the clones that the mutant calls define. cells = ds.snapshot_cell_selection("I") counts = ds.variants.import_counts("mgatk/final", cell_selection=cells) loaded = ds.variants.load_counts(counts, sites=[310, 3244]) loaded.counts.head() # each cell's reads of each base at the two positions ds.plots.variant_depth(counts, group_by="sample") deep = ds.select_cells(counts, low=10) informative = ds.variants.informative(counts, cells=deep, preset="atac") ds.plots.variant_statistics(informative) clones = ds.variants.clones(counts, variants=informative) ds.plots.allele_frequencies(informative, clones=clones) Import counts ============= "import_counts" recognises each path’s format from its contents: a folder by the names of its files, an HDF5 file by its groups and a table by its columns. "format=" names the format when the contents leave it open. +----------------------+----------------------+----------------------+----------------------+----------------------+ | "format" | Files | Sites and alleles | Unit | Depth counts | |======================|======================|======================|======================|======================| | ""mgatk"" | an mgatk folder from | every position of | reads | every base | | | 0.4.7 on: | the reference but | | | | | ".A.txt.gz" to | "N", with the four | | | | | ".T.txt.gz", " | bases | | | | | _refAllele.t | | | | | | xt" and ".dept | | | | | | hTable.txt" | | | | +----------------------+----------------------+----------------------+----------------------+----------------------+ | ""maegatk"" | a maegatk folder of | as mgatk | molecules (UMI | every base | | | the same files, | | consensus reads) | | | | always with six | | | | | | columns | | | | +----------------------+----------------------+----------------------+----------------------+----------------------+ | ""cellsnp_lite"" | a cellsnp-lite | each record, with | declared | DP plus OTH, or DP, | | | folder: | its reference, its | | the reference and | | | "cellSNP.base.vcf", | alternative and, | | alternative, without | | | "cellSNP.samples.tsv | with OTH, ""other"", | | OTH | | | ", | the other bases | | | | | "cellSNP.tag.AD.mtx" | | | | | | and ".DP.mtx", and | | | | | | ".OTH.mtx" where | | | | | | written | | | | +----------------------+----------------------+----------------------+----------------------+----------------------+ | ""mquad"" | an MQuad folder: | each passed variant, | declared | reference and | | | "passed_ad.mtx", | with its reference | | alternative | | | "passed_dp.mtx", "pa | and alternative | | | | | ssed_variant_names.t | | | | | | xt" and | | | | | | "BIC_params.csv", | | | | | | with the "cellSNP.sa | | | | | | mples.tsv" of the | | | | | | cellsnp-lite run it | | | | | | read | | | | +----------------------+----------------------+----------------------+----------------------+----------------------+ | ""vartrix"" | a folder of | each record with one | declared | reference and | | | VarTrix’s coverage | alternative, with | | alternative | | | mode: the input VCF, | both alleles; one | | | | | the barcode list and | record per position | | | | | the two matrices, | | | | | | the reference one | | | | | | with "ref" in its | | | | | | name | | | | +----------------------+----------------------+----------------------+----------------------+----------------------+ | ""mito"" | MiTo’s allele | each position, with | molecules where | every base, the | | | frequency matrix, an | its reference and | "uns" says so, | largest DP of the | | | ".h5ad" file with | alternatives | otherwise declared | position’s variants | | | "AD" and "DP" layers | | | | +----------------------+----------------------+----------------------+----------------------+----------------------+ | ""giftwrap"" | giftwrap’s | each probe, with the | UMIs | every gapfill | | | "counts.N.h5" or "co | gapfills that cells | | | | | unts.N.filtered.h5", | show or the manifest | | | | | or a folder of GEO’s | names | | | | | "probe_reads.tsv.gz" | | | | | | , "barcodes.tsv.gz" | | | | | | and "manifest.tsv" | | | | | | of one sample | | | | +----------------------+----------------------+----------------------+----------------------+----------------------+ | ""ironthrone"" | IronThrone-GoT’s "

.summTable.txt | "WT", "MUT" and | | | | | " | "ambiguous" | | | +----------------------+----------------------+----------------------+----------------------+----------------------+ | ""got"" | a per-cell GoT table | the target, with | UMIs | "WT" and "MUT" | | | of one target, as | "WT" and "MUT" | | | | | deposited on GEO: | | | | | | "BC" with | | | | | | "num.WT.call" and | | | | | | "num.MUT.call", | | | | | | "WT.calls" and | | | | | | "MUT.calls", or | | | | | | "wtUMIs" and | | | | | | "mutUMIs", and a | | | | | | "Genotype" call | | | | +----------------------+----------------------+----------------------+----------------------+----------------------+ | ""got_multi"" | GoT-Multi-ML’s "iron | each target, and | UMIs | every call | | | throne_out_pp.csv" | each target’s gene- | | | | | and "got_multi_out_r | sequence and probe- | | | | | efined.csv" | barcode channels | | | +----------------------+----------------------+----------------------+----------------------+----------------------+ | ""gotcha"" | Gotcha’s tables | each target, with | reads | "WT" and "MUT" | | | written with R’s | "WT" and "MUT" | | | | | "write.csv": | | | | | | "WTcount", | | | | | | "MUTcount" and | | | | | | "genotype_pred", or | | | | | | "_WTcount" and | | | | | | "_MUTcount" | | | | +----------------------+----------------------+----------------------+----------------------+----------------------+ | ""tapestri"" | a Tapestri Pipeline | each position of the | reads | every read at the | | | 2.0 ".dna+pr | variants that pass | | position, "DP" | | | otein.h5" or | and those that | | | | | ".dna.h5", whose "as | "variants" names, | | | | | says/dna_variants" | with its reference | | | | | holds the layers | and the variants’ | | | | | "AF", "DP", "NGT" | alternatives | | | | | and "GQ" | | | | +----------------------+----------------------+----------------------+----------------------+----------------------+ | ""dabseq"" | a DAb-seq ".hdf5" | each position of its | reads | every read at the | | | file: "AD", "RD", | variants, with its | | position, "DP" | | | "DP", "GQ", "GT", | reference and | | | | | "VARIANTS" and | alternatives | | | | | "CELL_BARCODES" | | | | +----------------------+----------------------+----------------------+----------------------+----------------------+ mgatk with "--emit-base-qualities" and maegatk both write six columns, "position,barcode,forward,forward quality,reverse,reverse quality", and count strands the other way round, so a six-column folder needs "format". maegatk counts reverse-aligned reads as forward, and so did mgatk’s "call" and "bcall" from 0.4.0 to 0.5.8: maegatk’s strands are swapped back on import, and old mgatk files need "swap_strands=True", since their files do not reveal the version. In RNA libraries, such as MAESTER’s, strand follows the transcript, so strand tests apply to DNA libraries such as mtscATAC-seq. Deposits rename and compress files, so files are matched by suffix. mgatk names its reference after the contig, and a deposit that prefixes the base files prefixes the reference alike, so the contig is the reference’s name without the base files’ prefix: "TF1_filtered.chrM_refAllele.txt.gz" beside "TF1_filtered.A.txt.gz" names "chrM", and "NC_012920.1_refAllele.txt" names "NC_012920.1". A name that holds another dot, such as a prefix that the base files do not carry, leaves the contig unclear and raises "ValueError". nf- MiTo’s "A.txt.gz" and "refAllele.txt", which name no contig, read as mgatk files of an unnamed contig. maegatk writes its reference beside its tables, while nf-MiTo’s MAESTER folders drop it; copy the "_refAllele.txt" of the run’s reference into such a folder. A one-target GoT table names its target after the file, without its extensions and "summTable", so "Patient1_TET2.2847_summTable.txt" counts the site "Patient1_TET2.2847". IronThrone’s per-UMI lists are not read; the file’s SHA-256 keeps them traceable. GoT-Multi counts each target in two evidence channels, which become the sites "BCOR (geneseq)" and "BCOR (probebc)", while its genotypes, "genotype_geneseq", "genotype_probebc" and "genotype_merged", and the refined table’s "genotype_final", concern the target "BCOR"; both tables import together. A header with one field fewer than the rows, as R’s "write.table" writes it, leaves out the row names, which name the cells when no column does; a table whose every row ends in its delimiter with nothing after it is refused, as clone tables are. A GIFT probe’s gapfills are read in the probe’s orientation, the reverse complement of the transcript. The manifest’s wild-type gapfill ("original_sequence", or "original_gap_probe_sequence" on GEO) is the reference, its expected mutant ("gap_probe_sequence") an alternative, and other gapfills, among them a control probe’s empty gapfill, have the role ""other"". Manifests often name neither, and some name one gapfill as both; that gapfill is the reference, and the source record lists such probes as "same_gapfills". GEO’s per-UMI tables give each cell’s UMIs, reads (the sum of "pcr_duplicate_count") and support (the mean of "percent_supporting") per probe and gapfill, as giftwrap’s HDF5 file does. Tapestri and DAb-seq -------------------- A Tapestri Pipeline 2.0 file lists the candidate variants of its panel, 27,719 in Mission Bio’s AML example, of which its default filters pass a few dozen ("filtered" 0), 40 there. The import reads the variants that pass and those that "variants" names: a variant as the file names it, such as ""chr1:115256513:G/T"", or a position, such as ""chr1:115256513"", which takes every variant there. Each position is a site, whose alleles are its reference and the alternatives of its chosen variants, so the variants of one position share its depth, "DP". Where their references differ, as a SNV and a deletion of one position do once Tapestri splits them from one record, the longest reference is the site’s and each shorter reference is extended, with its alternative, by the bases it lacks, as VCF writes several alternatives of one record: "T/G" beside "TC/T" reads as the allele "GC" of the reference "TC". A variant’s alternative reads are "AF × DP / 100", a whole number of reads; a frequency and depth that give another number raise "ValueError". Tapestri counts no reference reads, and the frequencies of the alternatives of one position may add up to more than 100%, so the depth bounds each alternative’s reads rather than their sum. "NGT" becomes the call column "NGT", with labels "0" to "3", and "GQ" its "quality", where the cell has reads at the position: a call without reads says nothing of the cell, and Tapestri writes some such cells as wild type at "GQ" 0. Rows with an empty "ALT", or alleles that are not sequences of bases, such as a "REF" of "*" or an "ALT" of "C+T", which no documentation explains, are skipped with a warning and counted in the source record’s "n_skipped_empty_alternative" and "n_skipped_other_alleles", with up to five in "skipped_examples"; a variant that "variants" names exactly but cannot be read raises "ValueError". An "ALT" of "*", VCF’s allele of a deletion that spans the position, is read as that allele. The source record also counts the file’s variants, those that pass, those imported and those that "variants" adds, which it names by the SHA-256 of their names, so names that add the same variants reuse the import. It names the samples, the panel, the genome and the pipeline version, and keeps the thresholds of the pipeline’s genotype filters, "gq_cutoff", "dp_cutoff" and "af_cutoff", 30, 10 and 20 in Mission Bio’s example. A file’s cells are named "__", as "TapestriToZarr" names the cells of the store it writes from the same file, so the cells of several samples stay apart. A DAb-seq file holds each cell’s alternative and reference reads ("AD", "RD"), its depth ("DP") and its genotype ("GT", with "GQ") at every variant, named ":::/", or "::/" by the pipeline’s mode for other species, with missing values written as -1. "GT" is read as a diploid genotype, coded as "NGT"; the pipeline’s haploid mode writes 1 for a mutant, which then reads as "heterozygous". It flags no variant, so every variant whose alleles are bases is read; "variants" may name variants by their full name or as "chrom:pos:ref/alt". The variants of one position count the same reference reads, which are stored once and must agree. The FLT3-ITD variants that the pipeline adds from ITDseek hold ITDseek’s depth in "DP" and "RD", so each is a named site of its own, such as "FLT3-ITD:chr13:28608251:T/TCTC", with alternative reads alone; the pipeline calls FLT3-ITD wild type in every cell where ITDseek found none, and those calls, without reads, are left out like Tapestri’s. Cells are named as "CELL_BARCODES" names them, such as "AACAACCTAGTCTTCTCG-1-all", the barcode, the tube and the sample. Before it reads a file’s candidate variants and cells, which it holds as text while it chooses the variants, the import checks them against the memory budget: about 700 bytes for each Tapestri variant (1,100 for DAb-seq, whose names it splits into more parts) and 300 for each cell, with two bytes for each byte of their longest texts, four for DAb-seq’s names, so Mission Bio’s example of 27,719 variants and 1,313 cells needs about 26 MB. Each block of cells then reads its cells’ IDs with its reads, so the import holds no cell’s ID beyond its block. counts = ds.variants.import_counts( "sample.dna+protein.h5", cell_selection=cells, variants=["chr1:115256513:G/T"], ) genotypes = ds.variants.genotypes(counts, method="source") Join barcodes to cells ---------------------- A cell’s ID is the path’s "barcode_prefix", the barcode and "barcode_suffix", matched exactly, and "barcode_pattern" keeps the part of each barcode that its "barcode" group matches, by the rules that Clones API reference describes. Barcodes that match no selected cell are counted, saved in the artifact’s "n_file_barcodes", "n_unmatched_barcodes" and "unmatched_barcode_examples" attributes, and named in a warning, which a reused import gives again. Cell Ranger ARC writes the gene expression barcode into the ATAC BAM’s "CB" tag, so mgatk run with "-bt CB" names cells as the RNA library does. A multiplexed Flex barcode holds 24 nucleotides, the GEM barcode and the probe barcode, so a giftwrap file of a multiplexed run joins a singleplex expression library with "barcode_pattern=r"^(?P[ACGT]{16})"" and "barcode_suffix="-1"". Several libraries ----------------- Several paths of one format import together, such as one mgatk folder per Cell Ranger ATAC library, each with its own prefix. Their shared sites are one site, so they must lie on the same reference: a contig with another length or sequence raises "ValueError", as do one sequence under two contig names, such as mgatk references renamed "lib1_chrM_refAllele.txt" and "lib2_chrM_refAllele.txt", a named and an unnamed contig, files that give one allele two roles, and, outside GIFT’s open set of gapfills, files that give a site different alleles. The paths must also count in one unit, define depth alike and give their counts the same fields, so four-column and six-column mgatk folders, or cellsnp-lite folders with and without OTH, import apart; a table of calls alone, such as GoT-Multi’s refined table, gives no counts. A cell that two files count at one site raises "ValueError". A call that a source repeats for one cell and allele, as GoT-Multi-ML’s refined table does for a few cells, is kept once and counted in the "n_repeated_calls" attribute, while two different calls raise "ValueError". Saved counts are reused when each file that a source’s reader takes holds the same bytes (reuse of imports), so moved files reuse them. counts = ds.variants.import_counts( ["lib1/mgatk/final", "lib2/mgatk/final"], cell_selection=cells, barcode_prefix=["lib1__", "lib2__"], ) Units and depth --------------- The unit of the counts is the format’s own, except for cellsnp-lite, MQuad and VarTrix, which count UMIs or reads depending on how they ran: their import needs "count_unit="umis"" or "count_unit="reads"". A MiTo matrix counts molecules when its "uns" says so: "uns['molecules']", which MiTo 2 records when it turns mtscATAC reads into molecules, or a "pp_method" of maegatk or ReDeeM. Other MiTo matrices, such as mgatk reads that MiTo 1 assembled, need "count_unit". Tools define depth two ways, which the artifact’s "depth_definition" records: every read or UMI at the site (""all""), or those of the reference and alternative only (""reference_and_alternative""). cellsnp-lite’s reference count is DP less AD, and its depth DP plus OTH. Depth 0 means the cell was not covered at the site, so an allele frequency there is undefined rather than 0. Signac’s "AlleleFreq" turns such cells into frequency 0; CyteArc stores no depth for them, so a reader can tell “not covered” from “covered, without the allele”, which VarTrix’s explicit zeros mark. The import sorts each file’s rows by cell, whatever their order, and rows that outgrow the memory budget go to sorted runs in the system’s temporary folder ("TMPDIR" on Linux and macOS), about 4 bytes per count row of mgatk’s files, removed when the import ends. mgatk’s four base files parse on up to four threads when the budget is 1 GiB or more. Before it reads a row, the import checks the selected cells’ IDs against the budget, as text and as their read, which decodes a stored chunk of IDs at the width of the longest; the parsed chunks, the sort buffer and the bands share what the barcode join leaves. Read counts =========== "load_counts" returns a "VariantCounts" record. Its "sites" table holds one row per site, indexed by label, such as "chrM:310", "310" on an unnamed contig, or "IDH2_Mut1": "contig", "position", "name", "channel", "reference", and "depth" and "cells", the depth summed over all of the artifact’s cells and the cells covered. "alleles" holds each site’s alleles and their roles, "cells" each cell’s "listed", "mean_depth" and "covered_sites", and "contigs" each contig’s name, length and the SHA-256 of its sequence. The long tables hold one row per stored entry: "depth" ("cell", "site", "depth") where a cell was covered, "counts" ("cell", "site", "allele", "count", then the fields the source has: "forward" and "reverse", "quality_forward" and "quality_reverse", NaN where the strand count is 0, "reads" and "support") where a cell showed an allele, and "calls" ("cell", "site", "allele", "source", "call", "quality") where the source made a call, on the mutant allele, with the call’s quality where the source gives one, such as Tapestri’s "GQ", and NaN elsewhere. A frequency is a count over the depth at its cell and site: loaded = ds.variants.load_counts(counts, sites=["chrM:310"]) table = loaded.counts.merge(loaded.depth, on=["cell", "site"]) table["frequency"] = table["count"] / table["depth"] "sites" takes site labels, or positions on an artifact of one contig, and "cells" a cell selection within the counts’ cells or cell IDs, as a list, a NumPy array or a pandas index; "cells=[]" reads only the sites with their sums. The tables are checked against the memory budget, at about 120 bytes per row, with the two stored chunks that a read’s first and last rows decompress whole, the cells’ IDs as text, whose read decodes a stored chunk of them, and 40 bytes for each cell of the counts, whose offsets and per-cell arrays are read whole, so loading a few cells of many counts needs room for every cell. Cells named by ID are found among the IDs of every cell of the counts, which a cell selection does not need. The variant_counts artifact =========================== The artifact holds a row per selected cell in its canonical "mean_depth", the cell’s depth summed over the sites and divided by their number, NaN for a cell that no file names, "covered_sites" and "listed", so "ds.select_cells(counts, low=10)" keeps cells with a mean depth above 10. mgatk scores the cells whose mean depth exceeds 10 too, but averages over every position of its reference, the "N" at 3107 of rCRS included, so a cell within a few thousandths of 10 can fall on the other side. Read per-cell values lists them, and Artifacts, lineage, and summaries API reference describes the sites, alleles and long tables, which are sorted by cell in row bands. Plot depth ========== "ds.plots.variant_depth(counts)" draws two panels: each site’s mean depth over the cells that the files name, on a log scale, along the contig for positions and as bars for named sites; and how many cells have each mean depth over the sites. "group_by" names a cell metadata column, such as the sample, and "groups" a cell-label artifact; with either, the second panel draws one histogram per group. Its tables are "sites" ("site", "contig", "position", "mean_depth", "covered_fraction") and "cells" ("cell", "group", "mean_depth", "covered_sites"). Find informative variants ========================= Most bases that differ from the reference in a few reads are sequencing errors, and a variant marks clones only when some cells carry it and others do not. "informative" scores every variant that a chosen cell shows, an alternative allele at a site such as "310T>C", with the statistics of mgatk’s "variant_stats" table, computed as mgatk 0.9.0 computes them. A cell’s allele frequency is its count of the allele over its depth at the site, and is undefined in a cell without depth there. +----------------------------------------------------+----------------------------------------------------+ | Column | What it holds | |====================================================|====================================================| | "mean" | the allele’s reads over the depth at its site, | | | both summed over the cells | +----------------------------------------------------+----------------------------------------------------+ | "variance" | the sample variance of the covered cells’ allele | | | frequencies | +----------------------------------------------------+----------------------------------------------------+ | "vmr" | "variance / (mean + 1e-11)", the variance-mean | | | ratio | +----------------------------------------------------+----------------------------------------------------+ | "n_cells_conf_detected" | cells with at least 2 reads of the allele on each | | | strand | +----------------------------------------------------+----------------------------------------------------+ | "strand_correlation" | the Pearson correlation of the forward and reverse | | | reads of the cells that show the allele | +----------------------------------------------------+----------------------------------------------------+ | "n_cells_over_5", "n_cells_over_10", | cells whose allele frequency reaches 5, 10, 20 and | | "n_cells_over_20", "n_cells_over_95" | 95% | +----------------------------------------------------+----------------------------------------------------+ | "max_heteroplasmy" | the highest allele frequency of a cell | +----------------------------------------------------+----------------------------------------------------+ | "mean_coverage" | the depth at the site summed over the cells, over | | | their number | +----------------------------------------------------+----------------------------------------------------+ | "n_cells_covered" | cells with depth at the site, which the variance | | | covers | +----------------------------------------------------+----------------------------------------------------+ | "mean_quality" | the mean base quality of the allele, as MAESTER | | | computes it: the mean over the cells and strands | | | with reads of the allele of their mean quality | +----------------------------------------------------+----------------------------------------------------+ A variant is informative when it passes every threshold that is set and no mask lists it. "preset" names a published rule, and a threshold given by name replaces the preset’s: +-----------------------------------+-----------------------------------+-----------------------------------+ | "preset" | Thresholds | Source | |===================================|===================================|===================================| | ""atac"" | "min_cells=5", | the thresholds that Lareau et al. | | | "min_strand_correlation=0.65", | (2020) set for mtscATAC-seq, as | | | "min_vmr=0.01" | Signac’s vignette applies them | +-----------------------------------+-----------------------------------+-----------------------------------+ | ""rna"" | "min_depth=20", "min_quality=30", | the depth and quality that Miller | | | "min_vmr=0.01" | et al. (2022) required of MAESTER | | | | variants, with Signac’s VMR | +-----------------------------------+-----------------------------------+-----------------------------------+ VMR, depth and quality must exceed their thresholds, and the confidently detected cells and the strand correlation must reach theirs, as the published rules compare them; a statistic that is undefined, such as the strand correlation of a variant that one cell shows, fails. Errors appear on one strand, so the strand tests separate them from variants in DNA libraries such as mtscATAC-seq. In RNA libraries, strand follows the transcript, so ""rna"" tests no strand. MAESTER then chose variants by percentiles of their allele frequencies, which mgatk’s columns do not hold; the VMR stands in for that test of variability and keeps many more variants: on 50 cells of MiTo’s MAESTER data, ""rna"" keeps 470 of the 5,664 variants they show, where MiTo’s own filter keeps 3, and a larger "min_vmr" keeps fewer. A threshold that the counts cannot test, such as a strand correlation of cellsnp-lite counts or a base quality of four-column mgatk files, raises "ValueError" before any row is read. "cells" restricts the statistics to a cell selection within the counts’ cells. mgatk scores cells whose mean depth exceeds 10, and "ds.select_cells(counts, low=10)" selects nearly the same cells, for the reason that the variant_counts artifact above gives. Signac scores the cells of its object, often filtered the same way. Cells that no file names are left out either way. The statistics read the stored rows in bands of cells within the memory budget, which counts each band’s rows, the two stored chunks that a band’s reads decompress beyond its rows, and 48 bytes for each cell of the counts, chosen or not, whose offsets are read whole. Masks ----- "mask" names lists of variants that artifacts mimic, in the coordinates of the revised Cambridge Reference Sequence, which hg38’s "chrM" and GRCh38’s "MT" share. A masked variant is scored but never informative. A listed variant matches by position and alternative allele, and the masks refuse counts whose contig is not 16,569 bases long, such as hg19’s "chrM". +-----------------------------------+-----------------------------------+-----------------------------------+ | "mask" | Variants | Source | |===================================|===================================|===================================| | ""numt"" | the 25 false positives that two | Supplemental Table S1 of | | | common polymorphic NUMTs, nuclear | Laricchia et al. (2022) | | | copies of mitochondrial DNA, | | | | cause when their reads align to | | | | the mitochondrial genome, such as | | | | "16293A>C" | | +-----------------------------------+-----------------------------------+-----------------------------------+ | ""rna_editing"" | the RNA-DNA differences of human | Bar-Yaacov et al. (2013) | | | mitochondrial transcripts, | | | | "295C>T", "2617A>G", "2617A>T", | | | | "13710A>G" and "13710A>T", which | | | | RNA libraries read as variants | | +-----------------------------------+-----------------------------------+-----------------------------------+ | ""artifact_prone"" | every variant at 301, 302, 310, | Laricchia et al. (2022) | | | 316, 3107 and 16182, homopolymers | | | | that sequencers misread and the | | | | reference’s "N", which gnomAD | | | | filters | | +-----------------------------------+-----------------------------------+-----------------------------------+ How the statistics differ from Signac’s --------------------------------------- Signac’s "IdentifyVariants" computes the same columns but gives a cell without depth an allele frequency of 0, and by default replaces the frequency of a cell with depth below 10 by the pooled mean before taking the variance. CyteArc keeps the frequency of an uncovered cell undefined, as mgatk does: the variance and the cell counts cover the covered cells. Where every cell is covered, the statistics equal Signac’s without its variance stabilisation; where coverage is sparse, a variant that covered cells rarely carry has a larger variance and VMR in CyteArc than in Signac, whose zeros for uncovered cells lower both. Signac also scores every alternative base at every position, while "informative" lists those that some chosen cell shows, and mgatk those seen on both strands. "load_informative" returns a "VariantSelection" record: its "variants" table holds one row per scored variant, indexed by a label such as "chrM:310T>C", with "site", "position", "reference", "alternative", the statistics, "masked" and "informative", and the record names the counts, the cells, the preset, the thresholds as applied, the masks and the number of cells, and says whether the counts give strands and base qualities. "ds.plots.variant_statistics(informative)" draws each variant’s VMR against its strand correlation, as mgatk draws it, with the informative and masked variants marked and the thresholds drawn as lines; for counts without strands, the x axis is the mean depth instead. Call clones =========== "clones(counts, variants=informative)" calls clones from the allele frequencies of the informative variants, as Signac’s "FindClonotypes" does. Each cell’s neighbourhood holds the cell itself and the "n_neighbors - 1" cells most similar to it by the cosine similarity of the square roots of their allele frequencies; two cells that share a neighbour are joined with the Jaccard index of their neighbourhoods, their shared cells over the cells in either, and joins below 1/15 are dropped, as Seurat’s "FindNeighbors" builds its shared-neighbour graph; CyteArc’s Leiden algorithm then partitions the graph at "resolution", where Signac uses Seurat’s "FindClusters". "n_neighbors=10" and "resolution=1" are Signac’s defaults. Every cell is compared with every other cell, so the time grows with the square of the cells, and "within" compares the cells of each group apart. Cells with collinear profiles, such as cells that each show one informative variant, tie at similarity 1 and take the same few cells as neighbours, so every pair of them shares those and the graph grows with the square of their number too; a graph that does not fit the memory budget raises "MemoryError" before anything is saved. Signac turns an uncovered cell’s allele frequency into 0, which makes a cell that was not read at a variant look like one that lacks it. CyteArc compares two cells over the variants that both cover, and a pair that shares no covered variant, or no alternative allele at the variants it shares, has similarity 0 and is never joined. Cells covered at no informative variant, cells that show none of them and cells similar to no other cell carry no evidence of a clone and get none, where Signac assigns every cell to a cluster; the artifact counts them in "n_cells_without_coverage", "n_cells_without_alternative" and "n_cells_without_neighbors". The cosine similarity ignores how frequent an allele is, so a cell that shows a clone’s variant in one read is joined to that clone, as in Signac. "within" names a cell metadata column, such as the donor: clones never join cells of two groups, each clone’s label starts with its group, as in "donor1:3", and cells without a group get no clone. The groups are those that the clone summaries group cells into by a column, read and coded within the memory budget, and the artifact records them by each cell’s group and the groups’ labels. The clones cover the cells that "variants" scored, over the counts’ cells, and are saved as a "clones" artifact, which "ds.clones.load_clones", every clone summary and "ds.plots.clone_abundance" read. The bands of frequencies count what the statistics’ bands count and 24 bytes more for each cell of the counts, its rows and group, and saving the clones is checked at 34 bytes for each of the counts’ cells, each clone’s label as a Python string and what saving clones holds, which stores every cell’s label at the width of the longest, so counts over many cells of which the files name few need room for every cell. "ds.plots.allele_frequencies(informative, clones=clones)" draws each cell’s allele frequency of each informative variant, cells ordered by clone under a bar of clone colours, with uncovered entries grey and apart from a frequency of 0. Without "clones", it draws the scored cells in their order; "max_cells" caps the cells drawn, chosen evenly along the order. Call genotypes ============== Targeted genotyping asks what each cell carries at a few known mutation sites. "genotypes" calls each cell’s genotype at the sites of a "variant_counts" artifact in one vocabulary and saves a "variant_genotypes" artifact, "ds.variants.load_genotypes" reads it as a "VariantGenotypes" record, "genotype_clones" calls clones from the combinations of mutant calls that cells carry, and "ds.plots.genotypes" draws the calls. counts = ds.variants.import_counts( ["Patient1_CALR_summTable.txt", "Patient1_SF3B1_summTable.txt"], cell_selection=cells, ) calls = ds.variants.genotypes(counts, method="threshold", min_mutant_umis=1) ds.variants.load_genotypes(calls).calls.head() # each cell's call at each site clones = ds.variants.genotype_clones(calls, within="donor") ds.plots.genotypes(calls, group_by="cell_type", layout_key="UMAP", site="Patient1_CALR") Calls ----- A call says what a cell carries at a site, and keeps the source’s own label, such as GoT’s "MUT": +----------------------------------------------------+----------------------------------------------------+ | Call | What the cell carries | |====================================================|====================================================| | "wild_type" | the wild type only | +----------------------------------------------------+----------------------------------------------------+ | "heterozygous" | the alternative allele on one copy | +----------------------------------------------------+----------------------------------------------------+ | "homozygous_alternative" | the alternative allele on every copy | +----------------------------------------------------+----------------------------------------------------+ | "alternative_present" | the alternative allele, without a claim of how | | | many copies carry it, as GoT’s UMIs of one | | | transcript give | +----------------------------------------------------+----------------------------------------------------+ | "missing" | no evidence at the site | +----------------------------------------------------+----------------------------------------------------+ | "filtered" | evidence that the rule did not call, such as | | | ambiguous UMIs alone | +----------------------------------------------------+----------------------------------------------------+ A cell that the calls do not list at a site is "missing" there. "cytearc.variants.CALLS" lists the calls in the order of their stored codes, and "cytearc.variants.SOURCE_CALLS" maps the labels of each tool’s call column to them. Methods ------- "method" chooses how the calls are made: * ""source"" adopts the calls that the tool wrote into the imported files: GoT’s "Genotype" on GEO, GoT-Multi’s "genotype_merged" and its refined "genotype_final", whose "FalsePositive", a mutant call that its model rejected, is "filtered", Gotcha’s "genotype_pred", whose "HET" is "heterozygous", and Tapestri’s "NGT" and DAb-seq’s "GT", whose 0 to 3 are "wild_type", "heterozygous", "homozygous_alternative" and no call. GoT’s "MUT" is "alternative_present". Tapestri and DAb-seq call each variant of a position apart, so a site with coordinates is genotyped once per alternative that holds calls, labelled by the site and its alleles, such as "chr1:115256513:G/A", with the site’s reference and the alternative as its "reference" and "alternative"; the label is Tapestri’s name of the variant unless its reference was extended, as "chr1:2000:T/G" beside "TC/T" reads as "chr1:2000:TC/GC". A named site, such as a GoT target or DAb-seq’s FLT3-ITD, is genotyped once. Where the files wrote several call columns for a site, each cell takes its call from the first that holds one for it: the refined call, then the merged call, then a single call column, then the channels’ calls, so the cells that GoT-Multi’s model did not refine keep their merged call. A label of no call, such as GoT-Multi’s "Unprofiled", is "filtered" where the cell holds UMIs at the site or at its channels and "missing" where it holds none, and a cell with UMIs but no call is "filtered". "min_quality" and "min_depth" filter the adopted calls, as Tapestri’s pipeline filters genotypes: a call whose quality, such as Tapestri’s "GQ", is below "min_quality" or missing, or whose cell’s depth at the site, or at the deepest of its channels, is below "min_depth", is "filtered", and keeps its label. Only Tapestri’s and DAb-seq’s calls carry a quality, so counts of other formats refuse "min_quality" with "ValueError" rather than filter every call, and so does "min_depth" at a genotyped site where no cell holds depth, as in a GoT table of calls alone. "min_quality=30, min_depth=10" applies the thresholds that Mission Bio’s example records as "gq_cutoff" and "dp_cutoff". Tapestri’s "FILTER_MASK" layer, which marks the calls that fail those thresholds or the allele frequency’s, "af_cutoff", is not read, since the frequency is the counts’ alternative reads over the depth. * ""threshold"" calls each cell from its UMIs of the site’s wild type and alternative, by the rule given, since GoT studies apply different ones, and refuses Tapestri’s counts, which hold no reference reads. "min_mutant_umis=1" calls a cell "alternative_present" with one mutant UMI or more and "wild_type" with none, the rule that GEO GSE204845 states for its patient samples and that Single-cell_BPDCN’s XV-seq analysis applies. "min_fraction=0.9" calls a cell mutant or wild type when that allele holds more than 90% of its UMIs, the rule that Nam et al. (2019) applied to cell lines, whose many transcripts make PCR recombination more likely; the two may be combined. "min_umis", 1 by default, is the fewest UMIs of the two alleles that a call needs, and a cell that no rule calls is "filtered". On GoT-Multi’s counts each channel of a target is a site of its own, so a cell’s "n_mutant_sites" counts a target once per channel. * ""gift"" applies the model of GIFT (Blattman et al. 2026) to giftwrap’s counts, as the next section describes. ""threshold"" and ""gift"" genotype the sites that have one wild type and one alternative, and ""source"" those with calls, unless "sites" names them. "alleles" names the wild type and alternative of a site whose files do not, such as a GIFT probe whose manifest lacks its gapfills: "alleles={"JAK2 c.1849G>T": ("CAC", "AAC")}", in the probe’s orientation, as the files give the gapfills. A wild type that no cell shows at the site raises "ValueError", and an alternative that none shows gives a warning. ""gift"" reads the gapfills of giftwrap’s counts and refuses counts of other tools. Saved genotypes are reused when the counts, the method and its rules are the same. GIFT’s model ------------ GIFT reads a mutation from the gapfill that a probe pair captures, and a UMI may come from template switching during PCR, a swap, which carries no genotype information. Each UMI scores each genotype with the larger of its swap probability and its gapfill’s probability under the genotype. A gapfill’s probability of the wild type is its likelihood under the wild type over the sum of its likelihoods under the wild type and the mutant, where each edit between the gapfill and an allele multiplies the likelihood by 0.01. GIFT gives any gapfill that does not match an allele a likelihood of 0.01 in cell lines, which equals this for a gapfill one edit away; the heterozygote scores the larger of half of each allele’s probability. A cell’s likelihood of a genotype is the product of its UMIs’ scores, summed as logarithms so that many UMIs do not underflow, and the posterior divides it by the sum over the three genotypes, a uniform prior. The most probable genotype is called when its posterior exceeds "min_posterior", 0.6 by default, and the cell holds at least "min_umis" UMIs, 1 by default, whose swap probability lies below 0.1; otherwise the cell is "filtered". Each call keeps the three posteriors. "swap_probabilities" gives the probability of a swap by a UMI’s reads, starting at 1 read, such as "{1: 0.5, 2: 0.2, 4: 0.05}", and a UMI takes the value of the most reads given that it reaches. GIFT’s reproducibility code stores each library’s likelihood that a UMI of so many reads is no swap, so give one minus it, and a probability of 1 for UMIs of few reads drops them, as GIFT drops them in cell lines. giftwrap stores each cell’s UMIs and reads per probe and gapfill rather than per UMI, so each UMI takes its entry’s reads per UMI, rounded down. Without "swap_probabilities" every UMI counts, as GIFT treats cell lines once their UMIs of few reads are dropped. GIFT learns the swap probabilities and the gapfill error rates of patient samples from controls; here the swap probabilities are given, and a gapfill that matches neither allele is weighed by its edits. The variant_genotypes artifact ------------------------------ The artifact holds a row per cell of the counts in its canonical "n_mutant_sites", the sites called "heterozygous", "homozygous_alternative" or "alternative_present", and "n_called_sites", those and the sites called "wild_type", so "ds.select_cells(calls, low=0)" keeps the cells with a mutant call. Artifacts, lineage, and summaries API reference describes its arrays. "load_genotypes" returns a "VariantGenotypes" record. Its "sites" table holds each genotyped site’s "reference" and "alternative", "source", the call columns that ""source"" adopted in the order a cell’s call is taken from them, and the cells of each call, those that the calls do not list counted "missing"; "cells" holds each cell’s "n_mutant_sites" and "n_called_sites" at the sites shown; and "calls" holds one row per cell and site with evidence: "cell", "site", "call", "label", "reference_umis" and "alternative_umis", -1 for adopted calls, and GIFT’s "posterior_wild_type", "posterior_heterozygous" and "posterior_homozygous_alternative", NaN for the other methods. "sites" reads the calls of some sites, which the "cells" table then counts. Calling reads the counts in bands of cells within the memory budget, which counts each band’s rows, about 144 bytes per row for ""threshold"", 168 for ""source"" and 360 for ""gift"", whose scores and posteriors hold the most; 72 bytes for each cell of the counts; the chunks that a band’s reads decompress and that appending its calls rewrites; and 1 MiB that the store’s objects hold whatever the counts. A band holds at least 1,000 rows, or every row, and every row of the cell that holds the most. Loading checks about 220 bytes per call row, 200 per cell and the read of the cells’ IDs. Clones from genotypes --------------------- Cells that descend from one mutated ancestor share its mutations. "genotype_clones(calls)" takes each cell’s mutant sites as a set and the sites at which it is "missing" or "filtered" as unknown rather than wild type. Two cells are linked when they agree at every site that both were called at, a Jaccard index of 1 over the sites that both know. Cells whose calls agree with two cells that disagree at a site they were not called at cannot be placed: cells that carry a founder mutation and were not read where a subclone differs from the rest of the founder’s cells could belong to either. Such cells are left out as ambiguous, however many there are, so they never join the founder and its subclone into one clone. A cell read where its clone’s partly read cells were not stays with them, since cells that share no read site do not disagree. The linked cells are the clones, numbered by size as "clone_1" and on, and a cell without a mutant call gets no clone. A doublet of two clones read at the sites where they differ carries a combination of its own and forms a clone of one cell. "sites" chooses the sites whose calls define the clones, and "within" names cell metadata columns, such as the donor, whose groups clones never cross. Exact agreement makes the choice of sites matter. Single-cell DNA amplification loses one allele of a cell at some sites, allelic dropout, which turns a heterozygous call into a wild-type or homozygous one, and a cell whose call turns disagrees with the rest of its clone. Every site adds such chances, while only the sites that some but not most cells carry tell clones apart. In Mission Bio’s example of four cell lines, the 40 variants that pass give 587 clones of at most 51 cells, and the 16 variants that 15% to 40% of the called cells carry give four clones of 194 to 351 cells, beside clones of 43 cells or fewer: genotypes = ds.variants.genotypes(counts, method="source") sites = ds.variants.load_genotypes(genotypes).sites mutant = sites["heterozygous"] + sites["homozygous_alternative"] share = mutant / (mutant + sites["wild_type"]) clones = ds.variants.genotype_clones( genotypes, sites=share[share.between(0.15, 0.4)].index.tolist() ) Filtering calls by "min_quality" and "min_depth" leaves more sites unknown, and so more cells ambiguous: with "min_quality=30, min_depth=10", the 17 variants that 15% to 40% of the called cells then carry leave 809 of the 1,313 cells ambiguous. The clones are saved as a "clones" artifact, which "ds.clones.load_clones", every clone summary and "ds.plots.clone_abundance" read; its "ambiguous" array marks the cells left out as ambiguous, and its attributes count the clones, the sites, the mutant cells and the ambiguous cells. The calls are read whole, at about 48 bytes per call row and 80 per cell, beside a mask of each mutant cell’s known sites and the 1 MiB that the store’s objects hold. Plot genotypes -------------- "ds.plots.genotypes(calls)" draws, for each site, one bar per group whose segments are the shares of the group’s cells with each call; cells that the calls do not list at a site count as "missing". "group_by" names a cell metadata column, such as the cell type, and "groups" a cell-label artifact; without either, each site has one bar of all cells, and cells without a group are left out. With "layout_key", such as ""UMAP"", a second panel draws the cells at the cell metadata columns "UMAP1" and "UMAP2", coloured by their call at "site", which may be left out when one site is drawn. Its tables are "calls" ("site", "group", "call", "cells", "share") and, with a layout, "cells" ("cell", "x", "y", "call"). VariantCounts ============= class cytearc.variants.VariantCounts(sites, alleles, cells, depth, counts, calls, contigs, format, count_unit, depth_definition) Per-cell allele evidence at variant sites, read from a "variant_counts" artifact. Long tables hold one row per stored entry, so a cell that a table does not list at a site was not covered there ("depth"), did not show the allele ("counts") or has no call from the source ("calls"). A frequency is a count over the depth at its cell and site, and is undefined where the depth is 0. sites One row per site, indexed by label (such as "chrM:310" or "IDH2_Mut1"): "contig", "position" (-1 for a named site), "name", "channel", "reference" (its reference allele, or None), "depth" (summed over all of the artifact’s cells) and "cells" (cells covered there). alleles One row per allele: "site", "allele" and "role". cells One row per cell, indexed by cell ID: "listed" (whether a file names the cell), "mean_depth" (over the sites, NaN for a cell no file names) and "covered_sites". depth "cell", "site" and "depth", one row per covered cell and site. counts "cell", "site", "allele" and "count", then the fields the source has: "forward" and "reverse" (strand counts), "quality_forward" and "quality_reverse" (mean base quality, NaN where the strand count is 0), "reads" and "support". calls "cell", "site", "allele", "source" (the column that held the call), "call", as the source wrote it, and "quality", the call’s quality where the source gives one, such as Tapestri’s genotype quality "GQ", and NaN elsewhere. contigs One row per contig: "name", "length" and "sequence_sha256". format The format of the imported files, such as ""mgatk"". count_unit ""reads"", ""umis"" or ""molecules"". depth_definition ""all"" or ""reference_and_alternative"". VariantSelection ================ class cytearc.variants.VariantSelection(variants, counts, cells, preset, thresholds, masks, n_cells, stranded, qualities) Per-variant statistics and the informative variants they select. variants One row per variant that some cell shows, indexed by its label (such as "chrM:310T>C"): "site", "position" (-1 for a named site), "reference" and "alternative"; mgatk’s "variant_stats" columns "mean", "variance", "vmr", "n_cells_conf_detected" and "strand_correlation" (NA without strand counts), "n_cells_over_5" to "n_cells_over_95", "max_heteroplasmy" and "mean_coverage"; "n_cells_covered", "mean_quality" (NaN without base qualities), "masked" and "informative". counts The "variant_counts" artifact that the statistics read. cells The cell selection that the statistics were restricted to, or None for every cell that the counts’ files name. preset The preset that set the thresholds, or None. thresholds Each threshold as the selection applied it, None for a test not made. masks The masks whose variants are never informative. n_cells The cells that the statistics cover. stranded Whether the counts give forward and reverse reads, and so "n_cells_conf_detected" and "strand_correlation". qualities Whether the counts give base qualities, and so "mean_quality". VariantGenotypes ================ class cytearc.variants.VariantGenotypes(sites, cells, calls, counts, method, rules) Genotype calls of the cells at targeted sites, read from a "variant_genotypes" artifact. A cell that the "calls" table does not list at a site holds no evidence there, which is the call "missing". sites One row per site shown, indexed by its label: "reference" and "alternative", the alleles the calls compare (empty for adopted source calls), "source" (the call columns that "method="source"" adopted, in the order a cell’s call is taken from them, else empty) and the cells of each call, those that the calls do not list counted "missing". cells One row per cell, indexed by cell ID: "n_mutant_sites" (sites shown called heterozygous, homozygous alternative or alternative present) and "n_called_sites" (those and the sites called wild type). calls One row per cell and site with evidence: "cell", "site", "call" (one of "cytearc.variants.CALLS", ordered), "label" (the source’s own label, empty for the other methods), "reference_umis" and "alternative_umis" (-1 where the counts hold none), and the posteriors "posterior_wild_type", "posterior_heterozygous" and "posterior_homozygous_alternative" (NaN unless "method="gift""). counts The "variant_counts" artifact the calls came from. method ""source"", ""threshold"" or ""gift"". rules The rules as the method applied them.