CAST-genomics · aryarm · Jan 22, 2022 · Jan 24, 2022 · Jan 26, 2022 · Feb 8, 2022
diff --git a/.gitignore b/.gitignore
@@ -9,4 +9,13 @@ __pycache__
 # pytest cache
 .pytest_cache
 # poetry
-dist/
+dist/
+
+# OSX
+*.DS_Store*
+
+# Test output files
+test.par
+test.phen
+example_simgenotype.bp
+example_simgenotype.vcf
diff --git a/README.md b/README.md
@@ -2,7 +2,31 @@
 Please wait until we have published our first tagged release before using our code.
 
 # haptools
-Simulate phenotypes for fine-mapping. Use real variants to simulate real, biological LD patterns.
-The Snakemake pipeline in the `snakemake/` directory uses the results of the simulation to test several fine-mapping methods, including FINEMAP and SuSiE.
+
+Haptools is a collection of tools for simulating and analyzing genotypes and phenotypes while taking into account haplotype information. It is particularly designed for analysis of individuals with admixed ancestries, although the tools can also be used for non-admixed individuals.
 
 Homepage: https://haptools.readthedocs.io/
+
+## Installation
+
+UNDER CONSTRUCTION
+
+## Haptools utilities
+
+Haptools consists of multiple utilities listed below. Click on a utility to see more detailed usage information.
+
+* [`haptools simgenotype`](docs/commands/simgenotype.md): Simulate genotypes for admixed individuals under user-specified demographic histories. 
+
+* [`haptools simphenotype`](docs/commands/simphenotype.md): Simulate a complex trait, taking into account local ancestry- or haplotype- specific effects. `haptools simphenotype` takes as input a VCF file and outputs simulated phenotypes for each sample.
+
+* [`haptools karyogram`](docs/commands/karyogram.md): Visualize a "chromosome painting" of local ancestry labels based on breakpoints output by `haptools simgenome`.
+
+Outputs produced by these utilities are compatible with each other. For example
+`haptools simgenome` outputs a VCF file with local ancestry information annotated for each variant. The output VCF file can be used as input to `haptools simphenotype` to simulate phenotype information. `haptools simgenome` also outputs a list of local ancestry breakpoints which can be visualized using `haptools karyogram`. 
+
+## Contributing
+
+If you are interested in contributing to `haptools`, please get in touch by submitting a Github issue or contacting us at mlamkin@ucsd.edu.
+
+
+
diff --git a/docs/api/haptools.rst b/docs/api/haptools.rst
@@ -36,3 +36,26 @@ haptools.data.phenotypes module
    :undoc-members:
    :show-inheritance:
 
+haptools.data.covariates module
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. automodule:: haptools.data.covariates
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+haptools.data.haplotypes module
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. automodule:: haptools.data.haplotypes
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+haptools.sim_genotype module
+~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. automodule:: haptools.sim_genotype
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/docs/commands/karyogram.md b/docs/commands/karyogram.md
@@ -0,0 +1,36 @@
+# Haptools karyogram
+
+`haptools karyogram` takes as input a breakpoints file (e.g. as output by `haptools simgenotype`) and a sample name, and plots a karyogram depicting local ancestry tracks.
+
+
+## Basic usage
+
+```
+haptools karyogram \
+  --bp tests/data/5gen.bp \
+  --sample Sample_1 \
+  --out karyogram.png
+```
+
+See details of the breakpoints file [here](simgenotype.md). If you specify `--sample $SAMPLE`, the breakpoints file must have breakpoints for `$SAMPLE_1` and `$SAMPLE_2` (the two haplotypes of `$SAMPLE`).
+
+## Additional options
+
+You may also specify the following options:
+
+* `--centromeres <FILE>`: path to a file describing the locations of chromosome ends and centromeres. An example file is given here: `tests/data/centromeres_hg19.txt`. The columns are: chromosome, chrom_start, centromere, chrom_end. For acrocentric chromosomes, the centromere field is ommitted. This file format was taken from [here](https://github.com/armartin/ancestry_pipeline).
+* `--colors "pop1:color1,pop2:color2..."`: You can optionally specify which colors should be used for each population. If colors are not given, the script chooses reasonable defaults.
+
+## Example command
+
+The following example can be run with files in this repository:
+
+```
+haptools karyogram --bp tests/data/5gen.bp --sample Sample_1 \
+       --out test_karyogram.png --centromeres tests/data/centromeres_hg19.txt \
+       --colors 'CEU:blue,YRI:red'
+```
+
+This will output a file `test_karyogram.png`. The example is shown below.
+
+![Example karyogram](../images/test_karyogram.png)
diff --git a/docs/commands/karyogram.rst b/docs/commands/karyogram.rst
@@ -0,0 +1,4 @@
+.. _subcommands-karyogram:
+
+.. include:: karyogram.md
+   :parser: myst_parser.sphinx_
diff --git a/docs/commands/simgenotype.md b/docs/commands/simgenotype.md
@@ -0,0 +1,87 @@
+# simgenotype
+
+`haptools simgenotype` takes as input a reference set of haplotypes in VCF format and a user-specified admixture model. It outputs a VCF file with simulated genotype information for admixed genotypes, as well as a breakpoints file that can be used for visualization.
+
+## Basic usage
+
+```
+haptools simgenotype \
+  --invcf REFVCF \
+  --sample_info SAMPLEINFOFILE \
+  --model MODELFILE \
+  --map GENETICMAP \
+  --out OUTPREFIX
+```
+
+Detailed information about each option, and example commands using publicly available files, are shown below.
+
+## Detailed usage
+
+`--invcf` - Input VCF file used to simulate specifiic haplotypes for resulting samples  
+`--sample_info` - File used to map samples in `REFVCF` to populations found in `MODELFILE`   
+`--model` - Parameters for simulating admixture across generations  
+`--map` - .map file used to determine recombination events during the simulation  
+`--out` - Output prefix of the structure `/path/to/output` which results in the vcf file `output.vcf.gz` and breakpoints file `output.bp`   
+
+## File formats
+
+Model Format
+
+Structure of model.dat file
+
+`num_samples` - Total number of samples to be output by the simulator (`num_samples*2` haplotypes)  
+`num_generations` - Number of generations to simulate admixture, must be > 0  
+`*_freq` - Frequency of populations to be present in the simulated samples
+
+```
+{num_samples} Admixed Pop_label1 Pop_label2 ... Pop_labeln
+{num_generations} {admixed_freq} {pop_label1_freq} {pop_label2_freq} ... {pop_labeln_freq}
+```
+
+Example model.dat file
+
+```
+40   Admixed   CEU    YRI
+6    0         0.2    0.8
+```
+Simulating 6 generations in this case implies the first generation has population freqs `Admixed=0, CEU=0.2, YRI=0.8` and the remaining 2-6 generations have population frequency `Admixed=1, CEU=0, YRI=0`  
+
+Map Format
+
+`chr` - chromosome of coordinate (1-22, X)  
+`var` - variant identifier   
+`pos cM` - Position in centimorgans   
+`pos bp` - Base pair coordinate  
+
+```
+{chr}\t{var}\t{pos cM}\t{pos bp}
+```
+Beagle Genetic Maps used in simulation (GRCh38): http://bochet.gcc.biostat.washington.edu/beagle/genetic_maps/
+
+
+Outfile Format
+
+`Sample Header` - Name of sample following the structure `Sample_{number}_{hap}` eg. `Sample_10_1` for sample number 10 haplotype 1  
+`pop` - Population label corresponding to the index of the population in the dat file so in the example above CEU = 1, YRI = 2  
+`chr` - chromosome (1-22, X)  
+
+```
+Sample Header
+{pop}\t{chr}\t{pos bp}
+...
+Sample Header 2
+...
+```
+
+## Examples
+
+Example Command
+```
+haptools simgenotype \
+  --model ./tests/data/outvcf_gen.dat \
+  --mapdir ./tests/data/map/ \
+  --chroms 1,2 \
+  --invcf ./tests/data/outvcf_test.vcf \
+  --sample_info ./tests/data/outvcf_info.tab \
+  --out ./tests/data/example_simgenotype
+```
diff --git a/docs/commands/simgenotype.rst b/docs/commands/simgenotype.rst
@@ -0,0 +1,4 @@
+.. _commands-simgenotype:
+
+.. include:: simgenotype.md
+   :parser: myst_parser.sphinx_
diff --git a/docs/commands/simphenotype.md b/docs/commands/simphenotype.md
@@ -0,0 +1,76 @@
+# simphenotype
+
+Haptools simphenotype simulates a complex trait, taking into account haplotype- or local-ancestry- specific effects as well as traditional variant-level effects. It takes causal effects and genotypes as input and outputs simulated phenotypes.
+
+Usage is modeled based on the [GCTA GWAS Simulation](https://yanglab.westlake.edu.cn/software/gcta/#GWASSimulation) utility.
+
+## Usage
+
+Below is a basic `haptools simphenotype` command:
+
+```
+haptools simphenotype \
+  --vcf <gzipped vcf file> \
+  --hap <gzipped hap file> \
+  --out <outprefix> \
+  < --simu-qt | --simu-cc > \
+  [simulation options]
+```
+
+Required parameters:
+
+* `--vcf <string>`: A bgzipped, tabix-indexed, phased VCF file. If you are simulating local-ancestry effects, the VCF file must contain the `FORMAT/LA` tag included in output of `haptools simgenotype`. See [haptools file formats](../../docs/formats/inputs.rst) for more details.
+
+* `--hap <string>`: A bgzipped, tabix-indexed HAP file, which specifies causal effects. This is a custom format described in more detail in [haptools file formats](../../docs/formats/haplotypes.rst). The HAP format enables flexible specification of a range of effect types including traditional variant-level effects, haplotype-level effects, associations with repeat lengths at short tandem repeats, and interaction of these effects with local ancestry labels. See [Examples](#examples) below for detailed examples of how to specify effects.
+
+* `--out <string>`: Prefix to name output files.
+
+* `--simu-qt` or `simu-cc` indicate whether to simulate a quantitative or case control trait. One of these two options must be specified.
+
+Additional parameters:
+
+* `--simu-rep <int>`: Number of phenotypes to simulate. Default: 1.
+
+* `--simu-hsq <float>`: Trait heritability. Default: 0.1.
+
+* `--simu-k <float>`: Disease prevalence (for case-control). Default: 0.1
+
+## Output files
+
+`haptools simphenotypes` outputs the following files:
+
+* `<outprefix>.phen`: Based on the phenotype files accepted by [plink](https://www.cog-genomics.org/plink/1.9/input#pheno). Tab-delimited file with one row per sample. The first and second columns give the sample ID. The third column gives the simulated phenotype. If `--simu-rep` was set to greater than 1, additional columns are output for each simulated trait. Example file:
+
+```
+HG00096	HG00096	0.058008375906919506
+HG00097	HG00097	0.03472768471423458
+HG00099	HG00099	-0.20850127859705808
+HG00100	HG00100	-0.21206803352471154
+HG00101	HG00101	0.3157913763428323
+```
+
+* `<outprefix>.par`: summarizes the frequencies and effects of simulated haplotypes. The columns are: haplotype ID (from the HAP file), haplotype frequency, and effect. Example file:
+
+```
+Haplotype	Frequency	Effect
+H-001	0.6	-0.2
+```
+
+<a name="examples"></a>
+## Examples
+
+#### Simulate a single haplotype-effect based on a 2 SNP haplotype:
+
+```
+haptools simphenotype \
+  --vcf tests/data/simple.vcf.gz \
+  --hap tests/data/simple.hap.gz \
+  --out test \
+  --simu-qt --simu-hsq 0.8 --simu-rep 1
+```
+
+based on this HAP file (available in `tests/data`)
+
+```
+H-001	1	10114	10116	1:10114:T:C,1:10116:A:G	T,G	*	-0.2
+```
diff --git a/docs/commands/simphenotype.rst b/docs/commands/simphenotype.rst
@@ -0,0 +1,4 @@
+.. _commands-simphenotype:
+
+.. include:: simphenotype.md
+   :parser: myst_parser.sphinx_
diff --git a/docs/commands/transform.rst b/docs/commands/transform.rst
@@ -0,0 +1,39 @@
+.. _commands-transform:
+
+
+transform
+=========
+
+Transform a set of genotypes via a list of haplotypes. Create a new VCF containing haplotypes instead of variants.
+
+The ``transform`` command takes as input a set of genotypes in VCF and a list of haplotypes (specified as a :doc:`.hap file </formats/haplotypes>`) and outputs a set of haplotype "genotypes" in VCF.
+
+Usage
+~~~~~
+.. code-block:: bash
+
+	haptools transform \
+	--region TEXT \
+	--sample SAMPLE \
+	--samples-file FILENAME \
+	--output PATH \
+	--verbosity [CRITICAL|ERROR|WARNING|INFO|DEBUG|NOTSET] \
+	GENOTYPES HAPLOTYPES
+
+Examples
+~~~~~~~~
+.. code-block:: bash
+
+	haptools transform tests/data/example.vcf.gz tests/data/example.hap.gz | less
+
+.. code-block:: bash
+
+	haptools transform -o output.vcf.gz -s NA12878 tests/data/apoe.vcf.gz tests/data/apoe4.hap
+
+Detailed Usage
+~~~~~~~~~~~~~~
+
+.. click:: haptools.__main__:main
+   :prog: haptools
+   :show-nested:
+   :commands: transform
diff --git a/docs/conf.py b/docs/conf.py
@@ -19,8 +19,8 @@
 # -- Project information -----------------------------------------------------
 
 project = "haptools"
-copyright = "2021, Arya Massarat, Michael Lamkin"
-author = "Arya Massarat, Michael Lamkin"
+copyright = "2021, Michael Lamkin, Arya Massarat"
+author = "Michael Lamkin, Arya Massarat"
 
 
 # -- General configuration ---------------------------------------------------
@@ -34,6 +34,7 @@
     "sphinx_rtd_theme",
     "numpydoc",
     "sphinx_click",
+    "myst_parser",
 ]
 
 # Add any paths that contain templates here, relative to this directory.
@@ -47,6 +48,8 @@
 # -- Extension configuration -------------------------------------------------
 autosummary_generate = True
 numpydoc_show_class_members = False
+# allow for both rst and md syntax
+source_suffix = ['.rst']
 
 # -- Options for HTML output -------------------------------------------------
 

diff --git a/docs/executing/inputs.rst b/docs/executing/inputs.rst