This is an implementation of various haplotype frequency spectrum statistics useful for detecting hard and soft selective sweeps in genomes. This program implements the following statistics:
- saltiLASSI: DeGriorgio and Szpiech (2022) PLoS Genetics 18: e1010134.
- LASSI: Harris and DeGiorgio (2020) MBE doi.org/10.1093/molbev/msaa115.
- H12: Garud et al. (2015) PLoS Genetics 11:e1005004.
- H2/H1: Garud et al. (2015) PLoS Genetics 11:e1005004.
- G123: Harris et al. (2018) Genetics 210:1429-1452.
- G2/G1: Harris et al. (2018) Genetics 210:1429-1452.
- Number of Unique Haplotypes at Locus
lassip accepts VCF files, either phased (default, "hap" output files) or unphased (set --unphased, "mlg" output files), with or without missing data. ***SEE CHANGELOG 19JAN2024 for update on how missing data is handeled. lassip reads one or more vcf files, each of which may contain one or more contigs; --vcf takes a list, a file's records must be grouped by contig, and no contig may appear twice across the input. You must provide a population file that specifies population IDs for each individual ID you wish to analyse. Only IDs listed in the population file will be analyzed, and if multiple populations are present, all statistics will be computed on a per-population basis.
Use --hapstats to compute H/G stats in sliding windows along the genome, whether or not data are phased determines whether H or G statistics are used.
Use --calc-spec to compute the top K haplotype frequency specra in sliding windows along a contig. Pass multiple spectra files (e.g. from multiple contigs) with --spectra and --lassi to run the CLR computation for detecting sweeps (Harris and DeGiorgio 2020) or with --salti to run CLR computation from (DeGiorgio and Szpiech 2021).
If only --hapstats is given, files are named <basename>.lassip.[hap|mlg].stats.gz with format:
<chr> <start> <end> <nSNPs> <ppos> <nHaps> <uniqHaps> <h12|g123> <h2h1|g2g1>
The final 4 columns are repeated for each population in the analysis, with the population code prepended on the appropriate header label.
If only --calc-spec is given with --vcf, files are named <basename>.lassip.[hap|mlg].spectra.gz with format:
<header line for use with lassip when reading with --spectra>
<chr> <start> <end> <nSNPs> <ppos> <nHaps> <uniqHaps> <hfs_1> ... <hfs_K>
The final K+2 columns are repeated for each population in the analysis, with the population code prepended on the appropriate header label. hfs_N gives the frequency of the Nth most common haplotype in the given window. ppos is the window's physical midpoint, or the window index when --dist-type nw was used; stage 2 names the same column pos.
If both --calc-spec and --hapstats are given with --vcf, files are named <basename>.lassip.[hap|mlg].spectra.gz with format:
<header line for use with lassip when reading with --spectra>
<chr> <start> <end> <nSNPs> <ppos> <nHaps> <uniqHaps> <h12|g123> <h2h1|g2g1> <hfs_1> ... <hfs_K>
The final K+4 columns are repeated for each population in the analysis, with the population code prepended on the appropriate header label. hfs_N gives the frequency of the Nth most common haplotype in the given window. ppos is the window's physical midpoint, or the window index when --dist-type nw was used; stage 2 names the same column pos.
The header line referred to above records what the file was written with, and is what lets stage 2 read it:
#phased <0|1> hapstats <0|1> wins <N> K <K> npop <P> <pop1> ... <popP> [contigs <C> <name> <rows> ...]
wins is the number of window rows in the file. The trailing contigs field is written only when a file holds more than one contig; it lists each contig with the number of rows it contributes, and stage 2 checks it against the rows themselves. Contigs are delimited by the chr column rather than by the file, so a file may hold any number of them provided its rows are grouped by contig and ascending within one, and no contig may appear twice across the files passed to --spectra. A file written by an earlier version has no contigs field and is read exactly as before.
Passing *.spectra.gz files to lassip with --spectra <file1> ... <fileN> will compute the LASSI CLR if --lassi is set and the saltiLASSI CLR if --salti is set. Output is a single file concatenating the results from all contigs, named <basename>.lassip.[hap|mlg].out.gz. For LASSI computations, one of two formats is output:
<chr> <start> <end> <nSNPs> <pos> <nHaps> <uniqHaps> <h12|g123> <h2h1|g2g1> <m> <T>
or
<chr> <start> <end> <nSNPs> <pos> <nHaps> <uniqHaps> <m> <T>
For saltiLASSI computations, one of two formats is output:
<chr> <start> <end> <nSNPs> <pos> <nHaps> <uniqHaps> <h12|g123> <h2h1|g2g1> <m> <A> <L>
or
<chr> <start> <end> <nSNPs> <pos> <nHaps> <uniqHaps> <m> <A> <L>
Depending on whether --hapstats was set when the *.spectra.gz files were generated. The final columns are repeated for each population in the analysis, with the population code prepended on the appropriate header label. m gives the inferred number of sweeping haplotypes, A gives the sweep “width” and T/L gives the CLR test statistic.
If --filter-level 2 is given (default), loci are filtered on a population basis and all monomorphic snps in a given population are filtered. As a consequence, total number of windows and the window coordinates will not necessarily match between pops, and so each population has results output in a separate file.
RELEASES
The Binaries workflow builds the distributable binaries. It runs when a release is published, when a version tag is pushed, and on demand from the Actions tab. It compiles on a clean Linux and a clean macOS runner, runs the regression suite against each binary it has just built, refuses to continue if the tag does not name the version the binary reports, and produces one archive per platform:
lassip-v<version>-linux-x86_64.tar.gz
lassip-v<version>-macos-arm64.tar.gz
each with a .sha256 beside it, and each containing the binary, this README, the LICENSE and the manual. An archive rather than a bare binary because a file downloaded through a browser arrives without its executable bit. The Linux build links the C++ runtime statically, so it runs on distributions older than the runner; libz and libc stay dynamic.
Where the archives end up depends on how the run started. A published release gets them attached to it. A pushed tag gets them attached to the release for that tag if one already exists, and otherwise leaves them as run artifacts with a note -- pushing a tag from a terminal is not the same as announcing a release, so the workflow does not create one. Publishing the release afterwards runs the workflow again and attaches them then, so either order works. A manual run always leaves them as run artifacts, which is the way to rehearse a release before tagging one.
Nothing is committed to bin/ by any workflow.
BUILDING AND TESTING
cd src && make # portable -O3 build
make ARCHFLAGS="-m64 -mmmx -msse -msse2" # the flags lassip <= 1.2.2 used, x86-64 only
make check # build, then run the regression suite in tests/
tests/run_tests.sh compares the current build against recorded outputs for both stages. Its only data input is example/YRI.chr22.vcf.gz; every fixture is derived from that at run time, so the suite runs from a clean checkout and no test data is committed. The cases cover phased, unphased, two populations, each --filter-level, missing genotypes, --avg-spec, and stage 2 under --lassi and --salti with each --dist-type. Stage-1 outputs are compared by hash; stage-2 outputs column by column, so a failure names the statistic that moved. Run it after any change to src/.
Usage: lassip --vcf <file> --pop <file> --calc-spec [--hapstats] --winsize <int> --winstep <int> --out <prefix>
lassip --spectra <file> [<file> ...] (--lassi | --salti | --avg-spec) --out <prefix>
lassip computes haplotype frequency spectrum statistics in two stages.
Stage 1 reads genotypes and writes a haplotype frequency spectrum per window:
lassip --vcf YRI.chr22.vcf.gz --pop YRI.ids.pop.txt --calc-spec --hapstats \
--winsize 117 --winstep 12 --k 10 --out YRI.chr22
Stage 2 reads those spectra (all contigs of a population at once) and writes
the likelihood ratio statistics:
lassip --spectra YRI.chr*.lassip.hap.spectra.gz --salti --out YRI
Flags belonging to the other stage are ignored, so check which stage you are in
if an option seems to have no effect.
Methods: saltiLASSI (DeGiorgio and Szpiech 2022, PLoS Genetics 18:e1010134),
LASSI (Harris and DeGiorgio 2020, MBE doi.org/10.1093/molbev/msaa115),
H12 and H2/H1 (Garud et al. 2015, PLoS Genetics 11:e1005004),
G123 and G2/G1, reported in place of H12 and H2/H1 under --unphased
(Harris et al. 2018, Genetics 210:1429-1452).
----------General----------
--threads <int>: The number of threads to spawn during computations.
Default: 1
--seed <int>: Seed for the shuffle used when clustering haplotypes that
carry missing data. Runs with the same seed, inputs and flags are reproducible
regardless of --threads, operating system or compiler. Set 0 to seed from the
system clock instead, which reproduces the non-reproducible behaviour of
lassip <= 1.2.2.
Default: 1
----------Input and output----------
--out <string>: The basename for all output files.
Default: outfile
--map <string>: A map file formatted <chr#> <locusID> <genetic pos> <physical pos>.
Sites in VCF not in map file will be interpolated.
Default:
--vcf <string1> ... <stringN>: One or more VCF files containing haplotype data.
Variants should be coded 0/1 and every file must contain all of the samples
named by --pop. Several contigs may be analysed in one run, given either as
one file per contig or as files holding several; the run then writes a single
spectra file per population covering all of them. A contig may not appear
twice, and a file's records must be grouped by contig.
Default:
--pop <string>: A file containing <ind ID> <pop ID>.
Default:
--spectra <string1> ... <stringN>: A list of spectra files for finalization.
Contigs are delimited by the chr column rather than by the file, so a file may
hold any number of them and the results are the same either way. A file's rows
must be grouped by contig and ascending within one, and no contig may appear
twice across the files given here.
Default:
--unphased <bool>: Set this flag to indicate data are unphased.
Default: false
----------Windows----------
--winsize <int>: The window size within which to calculate statistics.
Default: 0
--winstep <int>: The sliding window step size.
Default: 0
----------Statistics----------
--calc-spec <bool>: Set this flag to compute K-truncated haplotype
frequency spectra.
Default: false
--avg-spec <bool>: Set this flag to compute and output the average
K-truncated haplotype frequency spectrum from a set of .spectra files.
Default: false
--null-spec <string>: A file containing a null K-truncated
haplotype spectrum for use computing LASSI or saltiLASSI.
Default:
--lassi <bool>: Set this flag to use the LASSI method.
Default: false
--lassi-choice <int>: Set this flag to change the way LASSI
distributes mass across sweeping haplotype classes. Takes an integer in {1..5}.
Default: 4
--hapstats <bool>: Set this flag to calculate haplotype statistics.
Default: false
--salti <bool>: Set this flag to use the saltiLASSI method.
Default: false
--k <int>: Top K haplotypes for LASSI computations.
Default: 10
----------Filtering----------
--filter-level <int>: Filter monomorphic sites and sites
with missing data: 0-no filtering, 1-compute freq for all samples, 2-compute freq per pop.
Default: 2
--max-lmiss <double>: Filter loci with > this proportion of missing data.
Default: 0.10
--max-hmiss <double>: Drop haplotypes with > this proportion of missing data when computing the HFS.
Default: 0.20
--hap-cluster <string>: How haplotypes carrying missing genotypes are grouped into
classes. Missing sites act as wildcards, so such a haplotype can be compatible with
several distinct haplotypes and something has to choose. One of:
best-comp (default) each haplotype joins the most frequent class it is
compatible with, taking the best-observed haplotypes first.
Deterministic; --seed has no effect on it.
garud-shuffle the pre-1.3 behaviour: shuffle the haplotypes, then let each in
turn absorb every compatible one not yet claimed. The result
depends on the shuffle, so runs differ unless --seed is fixed.
soft-em EXPERIMENTAL. Divide an ambiguous haplotype's count across the
classes it could belong to in proportion to their frequencies,
iterated to convergence. Class sizes become fractional.
Without missing genotypes and at --match-tol 0 all three give the same classes.
Default: best-comp
--match-tol <int>: Group haplotypes with missing data into
the same class as a haplotype with no missing data if they have <= this many pairwise differences.
Default: 0
--keep-monomorphic <bool>: Set this flag to retain any monomorphic sites in the data.
Default: false
----------saltiLASSI----------
--dist-type <string>: Distance measure for saltiLASSI: bp, cm, nw.
Default: bp
--max-gap <double>: Widest gap in basepairs between two genetic map positions that
--dist-type cm will interpolate across. A window falling in a wider gap cannot be
placed and lassip stops rather than guess its genetic position. Real maps have
gaps at centromeres and other low-recombination regions, so raise this if your
map is sparse where your data are dense. Set 0 to interpolate across any gap.
Default: 3000000.00
--max-extend-bp <double>: Maximum distance in basepairs from core window to consider for saltiLASSI.
Default: 100000.00
--max-extend-cm <double>: Maximum distance in centimorgans from core window to consider for saltiLASSI.
Default: 0.05
--max-extend-nw <double>: Maximum distance in number of windows from core window to consider for saltiLASSI.
Default: 5.00
CHANGE LOG
17SEP2026 - v1.3.1.
Exit codes now distinguish a file that cannot be opened from a file whose contents are wrong. A missing or unopenable input, and an output path that cannot be created, exit 74; a file that opens and then turns out to hold the wrong thing still exits 65. The two are distinguishable at the same flag: --map with a path that does not exist gives 74, --map with a map that does not cover the contig gives 65. Bad command lines remain 64. One message was wrong as well: a --null-spec file that could not be opened for reading was reported as having failed to open "for writing".
The manual is rewritten for this version. LASSI-Plus-Manual.pdf described v1.1.2 and predated --seed, --hap-cluster, --max-gap, --version, the grouped help, the exit codes, the reproducibility work and multi-contig input.
It now has a source in the repository: doc/MANUAL.tex, built with
cd doc && make # writes ../LASSI-Plus-Manual.pdf
which needs a LaTeX installation and python3. The option reference is not
written by hand -- doc/gen-options.py generates doc/options.tex from
lassip --help, so it describes the binary it was built against, and
make check in doc/ fails if the two have diverged. The statistic
definitions and the structure are carried over from the v1.1.2 manual.
The old manual was written in Word and the .docx was never tracked here, which is why the PDF drifted three releases behind the program.
15SEP2026 - v1.3.0. Performance, reproducibility and CLI work on top of v1.2.2.
Reproducibility. Runs on data containing missing genotypes were not reproducible: the haplotype shuffle described in the 06FEB2025 entry was seeded from the system clock once per window, so repeated identical runs disagreed even at the default --match-tol 0, and even single-threaded. The seed is now --seed (default 1), mixed with each window's SNP boundaries, so a run depends only on the seed, the inputs and the flags -- not on --threads or on window order. --seed 0 restores clock seeding if you want to investigate how clustering changes across runs, which is what that feature was for. Note that which complete haplotype an incomplete one merges into is still order-dependent; the seed makes the order reproducible, it does not remove the dependence.
Speed and memory, measured on example/YRI.chr22.vcf.gz (278,607 SNPs, 216 haplotypes, --winsize 117 --winstep 12, 23,208 windows, K = 10):
stage 1 (--calc-spec --hapstats) 30.2 s -> 1.9 s at 1 thread
12.6 s -> 1.4 s at 8 threads
229 MB -> 87 MB peak memory
stage 2 (--salti) ~50 min -> 17 s at 8 threads
~6 h -> 114 s of CPU
1,009 MB -> 48 MB peak memory
Genotypes are now stored two bits per locus rather than one char, the VCF is read once rather than twice, and each window's haplotype frequency spectrum is tallied from the packed window rather than from a freshly built character string per haplotype. Spectra are byte-identical to v1.2.2. For saltiLASSI, m and the likelihood ratio are unchanged; the reported A moves to an adjacent point on the 101-point grid in about 1% of windows, because the likelihood is flat in A there and which grid point wins is decided by rounding -- v1.2.2 has the same instability, not a stable value that has been perturbed.
--match-tol now matches its documentation. The help says haplotypes are grouped when they have "<= this many pairwise differences", but the test was d < MATCH_TOL, so every setting behaved as the one below it and the default --match-tol 0 pooled nothing at all. It is now d <= MATCH_TOL. On data with missing genotypes the new --match-tol t reproduces the old --match-tol t+1 exactly (mean unique haplotypes per window on a test fixture: old 81.0 / 77.1 / 48.8 at t = 0 / 1 / 2, new 77.1 / 48.8 / 30.7). Output is unchanged for --match-tol 0 on data without missing genotypes, which is the default path. Note that the 06FEB2025 entry below describes the old rule ("< MATCH_TOL"); the flag help was the correct description of the intent.
Threading. The workers ran through a cast to void ()(void *), which is undefined behaviour, their work orders were never freed, and the null-window counter was incremented from every thread without synchronisation -- so the window count in the spectra header, which the finalize stage reads, was wrong and irreproducible whenever a run had null windows and --threads > 1. Workers now run as std::thread with their real signatures, work orders are owned by the caller, and each thread counts null windows into its own slot. Windows are claimed from an atomic cursor rather than by a fixed stride, and per-window scratch buffers are per-thread and reused. The build now passes -std=c++11 explicitly (override with make STDFLAG=...); the code has used C++11 library features since 1.2.0 but compiled as C++98 by default.
Structure. main() was 440 lines holding flag registration, flag unpacking,
165 lines of validation and both pipelines; it is now 23 lines reading
registerFlags -> parseCommandLine -> readConfig -> validate -> run, with the
flags unpacked once into a Config struct. Stage-1 results were seven
map<string,T*> keyed by population name, looked up by name inside the window
loops; they are now fields of one struct per population, held in a vector.
Haplotype clustering is selectable, and the default changed. When a haplotype has missing genotypes it can be compatible with several distinct haplotypes -- missing sites act as wildcards -- and something has to decide which class it joins. That decision used to be made by shuffling the haplotypes and letting each in turn absorb whatever it reached first, so the spectrum depended on the shuffle: on one window of the test data with 5% of genotypes missing, 120 runs differing only in --seed gave H12 between 0.066 and 0.100 at --match-tol 0, and between 0.132 and 0.390 at --match-tol 1.
--hap-cluster now selects the rule:
best-comp (default) each haplotype joins the most frequent class it is
compatible with, taking the best-observed haplotypes first.
Deterministic, and in simulation the most accurate of the
hard-assignment rules at every missing rate above 1%.
garud-shuffle the pre-1.3 behaviour, for reproducing earlier output.
soft-em EXPERIMENTAL. Splits an ambiguous haplotype's count across
the classes it could belong to in proportion to their
frequencies. Class sizes become fractional.
Without missing genotypes and at --match-tol 0 all three produce identical output, so runs on complete data are unaffected. On chr22 of the example with 2% of genotypes blanked, best-comp is also 9x faster than garud-shuffle (11.4 s vs 107.8 s at one thread), because it compares each haplotype against the existing classes rather than against every other haplotype.
Effect on --unphased data. Multilocus genotype strings are drawn from {0,1,2,-} rather than {0,1,-}, so two of them agreeing at every jointly observed site is much less likely. On one 50-locus window with 5% of genotypes missing, the phased reading had 149 distinct haplotypes of which 72 were compatible with more than one other, and the three rules returned 83, 82 and 82 classes; the unphased reading of the same data had 81 distinct genotype strings of which 3 were ambiguous, and all three rules returned 74. So the default change has much less effect on --unphased runs, at least at --match-tol 0.
Multiple contigs in one spectra file. Stage 2 delimits contigs by the chr
column rather than by the input file, so --spectra accepts one file holding
several contigs as readily as one file per contig, and the two give identical
results. Rows must be grouped by contig and ascending within one; a file that
interleaves contigs, repeats one, or disagrees with its own window count is
refused. The header may carry an optional 'contigs <n> <name> <nwins> ...'
field after the population names, validated against the rows when present.
Stage 1 takes several contigs in one run. --vcf is a list flag, so
lassip --vcf chr1.vcf.gz chr2.vcf.gz --pop pops.txt --calc-spec \
--winsize 117 --winstep 12 --out scan
writes scan.<pop>.lassip.hap.spectra.gz covering both contigs, which stage 2
reads as one file. A single VCF holding several contigs works too, and the
two forms mix; what is not allowed is the same contig twice, or a contig
whose records are not all together in one file -- lassip reads a VCF in one
pass and will not buffer contigs to reassemble them, so sort by position
first. Contigs are processed one at a time and each one's genotypes are
released before the next is read, so peak memory grows sub-linearly and
flattens. On one to five copies of the chr22 example the median of three runs
each is 86, 141, 186, 211, 211 MB; how the contigs are packaged does not
matter, three contigs giving 156 MB as three files and 172 MB as one, a
difference smaller than the spread between repeats of either. The rows are
identical to what separate per-contig runs produce, concatenated.
This is also a bug fix. Previously a file holding more than one contig was analysed as a single block, and saltiLASSI's flanking-window scan crossed the contig boundary: with --dist-type nw the reported positions of every window after the first contig were wrong, and with any --dist-type the windows within --max-extend of a boundary took neighbouring-contig windows into the likelihood. Anyone who concatenated spectra files by hand should rerun. Stage 1 writes multi-contig spectra directly, so the hand concatenation that motivated the fix is no longer necessary.
Two things to know when moving between versions. A spectra file written by this version and read by lassip <= 1.2.2 is not refused: the contigs field is past the population names, where the older positional header parse stops, so an older binary reads the file and analyses every contig in it as one block -- the same silent error described above. If you share spectra files with someone on an older version, split them by contig first. In the other direction, files written by older versions are read exactly as before: they carry no contigs field, and one file per contig is still a supported input.
Genetic map handling. On the --dist-type cm path, getMapInfo could interpolate between two uninitialised indices when its forward scan found no bracketing pair, and it dereferenced a null pointer -- a segfault -- when the spectra named a contig the map does not. It now brackets by binary search and returns false when it cannot place a window. fillCMDistance had been ignoring that return value and writing the previous window's genetic position instead: with a map holding a 6 Mb hole, 3,951 windows of the YRI example came out sharing two genetic positions. Unplaced windows are now counted and reported, and the run stops. The gap threshold itself, previously hardcoded at 3 Mb, is now --max-gap: raise it for a map that is sparse where the data are dense, or set 0 to interpolate across any gap.
Bug fixes. The null-spectrum check at the --spectra stage used the K from the command line rather than from the spectra file, so its threshold was wrong unless you repeated the same --k. --calc-spec --dist-type cm wrote a column with no name in the header. A GT field shorter than three characters (haploid or truncated) was read past its end. A VCF containing more than one chromosome was silently labelled with the first record's contig; it is now rejected. An error raised while reading input escaped main and aborted the process instead of exiting with a code. A record beginning with whitespace was misread: the rewritten reader took the CHROM token from the first character of the line, so a leading space or tab gave an empty contig name and shifted every remaining field by one, putting the FORMAT column where the first sample's genotype belonged and failing with "Alleles must be coded 0/1/. only" -- an error naming the genotypes, which were fine. Leading whitespace is now skipped, a record with no CHROM field is reported as such, and the allele-coding error names the genotype, sample and position it read. lassip <= 1.2.2 took the fixed columns by stream extraction, which skips whitespace, so the fields were not shifted there; that concerns field extraction only, and whether any particular file completed a run depended on the usual filtering.
CLI. --help goes to stdout and exits 0 (-h also works), --version prints the version, a bare lassip prints a usage synopsis, and the help is grouped by topic with worked examples for both stages. Flags belonging to the other stage now warn instead of being silently ignored. Exit codes follow sysexits: 64 for a command line error, 65 for bad input data.
Housekeeping. 1,357 lines of unreachable code removed. The Makefile tracks header dependencies, keeps the x86-only architecture flags in an overridable ARCHFLAGS (empty by default, so builds work on arm64), and has a check target. tests/ holds a regression suite.
The example outputs are regenerated. example/YRI.chr22.lassip.hap.spectra.gz
dated from v1.0.0 and predated the ppos column; it is replaced by
example/YRI.chr22.YRI.lassip.hap.spectra.gz, which is what the script
actually writes -- at the default --filter-level 2 each population is
filtered separately and gets its own file, <out>.<pop>.lassip.hap.spectra.gz.
The script had two faults: it passed --lassi to stage 1, where it does
nothing, and without --calc-spec wrote .stats.gz; and its second command read
the population-free filename, so it silently consumed the stale 2020 file
rather than the one it had just produced. Both are fixed, and running
example/do_lassip_YRI.bash now reproduces the committed outputs.
08AUG2025 - v1.2.2. Fixed more bugs under unusual use cases related to filtering monomorphic sites and missing data. Missing data filtering and monomorphic filtering is now done separately, and you will receive a message with the number of sites filtered for each reason. There is a new command line argument --keep-monomorphic, if you wish to skip the monomorphic filtering step.
06FEB2025 - v1.2.1. Fixed a bug related to filtering monomorphic sites when using --unphased. I have also simplified the haplotype clustering method in the context of missing data. Haplotypes are merged if their hamming distance is < MATCH_TOL. Locations where at least 1 haplotype has a missing genotype call are ignored. This follows the approach from https://github.com/ngarud/SelectionHapStats/blob/master/scripts/H12_H2H1.py. NOTE: This approach can result in different hapotype cllusters depending on the order the haplotypes are processed. Generally, this seems to have relatively small impact on the number of unique haplotypes inferred, but I have implemented a random shuffle such that on subsequent runs the haplotype order is different. You can use this feature to investigate how clustering might change across runs.
29FEB2024 - Fixed more bugs that appear in the 1.1.2a version. DO NOT USE 1.1.2a.
18JAN2024 - v1.2.0. Changing the way missing data is handled. Three new command line arguments:
--max-hmiss <double>: Drop haplotypes with > this proportion of missing data when computing the HFS.
Default: 0.20
Previous behavior is that as soon as a missing genotype is encountered this haplotype is excluded from HFS calculation. You can now modify this behavior with --max-hmiss. A haplotype will only be dropped if there is more than MAX_HMISS proportion of alleles missing in that haplotype. e.g. 10-00--011 is a 10-snp haplotype with missing data represented by '-'. The default MAX_HMISS is 0.2, and since 3/10 > 0.2 this haplotype would be excluded. 10-0011011 has only 1/10 missing data, so it is kept.
--max-lmiss <double>: Filter loci with > this proportion of missing data.
Default: 0.10
Previous behavior did not filter loci based on missing data. Here you can set the proportion of missing data beyond which a locus will be excluded. --filter-level controls how this is calculated, 0 fir no filtering at all, 1 to calculate the denominator from all provided samples, 2 to calculate denominator based on population groups.
--match-tol <int>: Group haplotypes with missing data into
the same class as a haplotype with no missing data if they have <= this many pairwise differences.
Default: 0
The current implementation now tries to cluster any haplotypes that have missing data (that pass the MAX_HMISS check) with haplotypes in the data that are complete (i.e. no missing data). To do this, I calculate # of pairwise differences, ignoring any locations with missing data, and cluster the missing-data haplotype with the complete haplotype(s) to which it has the fewest differences (enforcing that the number of differences must at least be <= MATCH_TOL). If there are ties, it is counted fractionally with equal contribution to each matching haplotype. Let me illustrate:
Complete haps
A 101
B 111
C 110
Haps w/missing data
D 1-1
E 0-1
Let d(X,Y) be the # pairwise differences and MATCH_TOL == 0, then d(D,A) = 0, d(D,B) = 0, d(D,C) = 1. In this case D would be counted 0.5 as A and 0.5 as B.
d(E,A) = 1, d(E,B) = 1, d(E,C) = 2. E would be kept as its own haplotype class.
Let MATCH_TOL == 1, then D would still be counted as 0.5 towards both A and B. The clustering always prefers the fewest number of differences even if MATCH_TOL would allow more. E would also cluster 0.5 with A and 0.5 with B.
Even if MATCH_TOL > 0, complete haplotypes are never merged (this could be changed, thinking about it). Haps w/ missing data are never merged among themselves (some non-trivial issues here, but thinking about solutions).
03DEC2023 - v1.1.2a Reverting the changes from v1.1.2 regarding grid search. This introduced more bugs than it helped to fix. Please do not use v1.1.2.
12SEPT2023 - v1.1.2 Grid search across epsilon now starts at 0 and not 1/(100K). This fixes a bug that tends to manifest when there is lots of missing data. When the average/null truncated HFS has its smallest frequency class p_K < 1/(100K), this created a condition where the likelihoods were not computed and only placeholder values were output for all windows, since we set p_K = U but also expected p_K = U >= 1/(100*K).
03DEC2021 - Bug fixes relating to genetic map distance interpolation.
02SEPT2021 - v1.1.1 Introduced saltiLASSI options for measuring distance in bp (basepairs), nw (number of windows), or cm (centimorgans). New command line flags:
--max-extend-nw <double>: Maximum distance in number of windows from core window to consider for saltiLASSI.
--max-extend-cm <double>: Maximum distance in centimorgans from core window to consider for saltiLASSI.
--map <string>: A map file formatted <chr#> <locusID> <genetic pos> <physical pos>. Sites in VCF not in map file will be interpolated.
--dist-type <string>: Distance measure for saltiLASSI: bp, cm, nw.
These flags are used with --salti and --spectra flags when finalizing the computations.
07MAY2021 - v1.1.0 Introduced saltiLASSI framework. Must specify --lassi or --salti when passing spectra files with --spectra in order to choose a method. --filter-level default changes to 2. --null-spec introduced. Can be used to provide a null HFS for LASSI or saltiLASSI methods. --avg-spec introduced. When used with --spectra will compute and output an averaged HFS which can be passed with --null-spec to use as background null HFS. --max-extend-bp introduced. Maximum distance in basepairs from core window to consider for saltiLASSI.
10JUN2020 - v1.0.2a Fixed a bug that caused chromosome name to be blank with filter level > 0.
07JUN2020 - v1.0.2 Memory improvements for filter levels < 2.
06JUN2020 - v1.0.1 Options for filtering monomorphic loci with --filter-level 0 - no filtering 1 - filter all loci that are monomorphic in the samples selected for analysis 2 - filter monomorphic loci on a population by population basis
Also fixed error in computation of H2/H1 and G2/G1.
26MAY2020 - Initial release v1.0.0