Merge pull request #64 from gymrek-lab/ref/simphenotype

ref: simphenotype

README.md

@@ -9,7 +9,10 @@ Homepage: https://haptools.readthedocs.io/
 
 ## Installation
 
-UNDER CONSTRUCTION
+We have not officially published `haptools` yet, but in the meantime, you can install it directly from our Github repository.
+```bash
+pip install git+/~https://github.com/gymrek-lab/haptools.git
+```
 
 ## Haptools utilities
 
@@ -21,12 +24,13 @@ Haptools consists of multiple utilities listed below. Click on a utility to see
 
 * [`haptools karyogram`](docs/commands/karyogram.md): Visualize a "chromosome painting" of local ancestry labels based on breakpoints output by `haptools simgenome`.
 
+* [`haptools transform`](https://haptools.readthedocs.io/en/latest/commands/transform.html): Transform a set of genotypes via a list of haplotypes. Create a new VCF containing haplotypes instead of variants.
+
 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.
-
-
+We gladly welcome any contributions to `haptools`!
 
+Please read [our contribution guidelines](https://haptools.readthedocs.io/en/latest/project_info/contributing.html) and then submit a [Github issue](/~https://github.com/gymrek-lab/haptools/issues).

docs/api/haptools.rst

@@ -52,8 +52,16 @@ haptools.data.haplotypes module
    :undoc-members:
    :show-inheritance:
 
+haptools.sim_phenotype module
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. automodule:: haptools.sim_phenotype
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
 haptools.sim_genotype module
-~~~~~~~~~~~~~~~~~~~~~~~~
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 .. automodule:: haptools.sim_genotype
    :members:

docs/commands/karyogram.md

@@ -1,4 +1,4 @@
-# Haptools karyogram
+# 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.
 

docs/commands/simphenotype.rst

@@ -1,4 +1,81 @@
 .. _commands-simphenotype:
 
-.. include:: simphenotype.md
-   :parser: myst_parser.sphinx_
+
+simphenotype
+============
+
+Simulates a complex trait, taking into account haplotype- or local-ancestry- specific effects as well as traditional variant-level effects. The user denotes causal variables to use within the simulation by specifying them in a ``.hap`` file.
+
+The implementation is based on the `GCTA GWAS Simulation <https://yanglab.westlake.edu.cn/software/gcta/#GWASSimulation>`_ utility.
+
+Usage
+~~~~~
+.. code-block:: bash
+
+   haptools simphenotype \
+   --replications INT \
+   --heritability FLOAT \
+   --prevalence FLOAT \
+   --region TEXT \
+   --sample SAMPLE \
+   --samples-file FILENAME \
+   --output PATH \
+   --verbosity [CRITICAL|ERROR|WARNING|INFO|DEBUG|NOTSET] \
+   GENOTYPES HAPLOTYPES
+
+Model
+~~~~~
+Each normalized haplotype :math:`\vec{Z_j}` is encoded as an independent causal variable in a linear model:
+
+.. math::
+
+   \vec{y} = \sum_j \beta_j \vec{Z_j} + \vec \epsilon
+
+where
+
+.. math::
+
+   \epsilon_i \sim N(0, \sigma^2)
+
+.. math::
+
+   \sigma^2 = Var[\sum_j \beta_j \vec{Z_j}] * (\frac 1 {h^2} - 1)
+
+The heritability :math:`h^2` is user-specified, but if it is not provided, then :math:`\sigma^2` will be computed purely from the effect sizes, instead:
+
+.. math::
+
+   \sigma^2 = \Biggl \lbrace {1 - \sum \beta_j^2 \quad \quad {\sum \beta_j^2 \le 1} \atop 0 \quad \quad \quad \quad \quad \text{ otherwise }}
+
+If a prevalence for the disease is specified, the final :math:`\vec{y}` value will be thresholded to produce a binary case/control trait with the desired fraction of diseased individuals.
+
+Output
+~~~~~~
+Phenotypes are output in the PLINK2-style ``.pheno`` file format. If ``--replications`` was set to greater than 1, additional columns are output for each simulated trait.
+
+Note that case/control phenotypes are encoded as 0 (control) + 1 (case) **not** 1 (control) + 2 (case). The latter is used by PLINK2 unless the ``--1`` flag is used (see `the PLIN2 docs <https://www.cog-genomics.org/plink/2.0/input#input_missing_phenotype>`_). Therefore, you must use ``--1`` when providing our ``.pheno`` files to PLINK.
+
+Examples
+~~~~~~~~
+.. code-block:: bash
+
+   haptools simphenotype -o simulated.pheno tests/data/example.vcf.gz tests/data/simphenotype.hap
+
+Simulate two replicates of a case/control trait that occurs in 60% of your samples with a heritability of 0.8. Encode all of the haplotypes in ``tests/data/example.hap.gz`` as independent causal variables.
+
+.. code-block:: bash
+
+   haptools simphenotype \
+   --replications 2 \
+   --heritability 0.8 \
+   --prevalence 0.6 \
+   --output simulated.pheno \
+   tests/data/example.vcf.gz tests/data/example.hap.gz
+
+Detailed Usage
+~~~~~~~~~~~~~~
+
+.. click:: haptools.__main__:main
+   :prog: haptools
+   :show-nested:
+   :commands: simphenotype

docs/commands/transform.rst

@@ -24,7 +24,7 @@ Examples
 ~~~~~~~~
 .. code-block:: bash
 
-	haptools transform tests/data/example.vcf.gz tests/data/example.hap.gz | less
+	haptools transform tests/data/example.vcf.gz tests/data/basic.hap.gz | less
 
 .. code-block:: bash
 

docs/formats/genotypes.rst

@@ -0,0 +1,7 @@
+.. _formats-genotypes:
+
+
+Genotypes
+=========
+
+Genotype files must be specified as VCF or BCF files.

docs/formats/haplotypes.rst

@@ -1,8 +1,8 @@
 .. _formats-haplotypes:
 
 
-.hap
-====
+Haplotypes
+==========
 
 This document describes our custom file format specification for haplotypes: the ``.hap`` file.
 

docs/formats/phenotypes.rst

@@ -0,0 +1,28 @@
+.. _formats-phenotypes:
+
+
+Phenotypes and Covariates
+=========================
+
+Phenotype file format
+---------------------
+Phenotypes are expected to follow `the PLINK2 .pheno file format <https://www.cog-genomics.org/plink/2.0/input#pheno>`_. This is a
+tab-separated format where the first column corresponds to the sample ID, and
+subsequent columns contain each of your phenotypes.
+
+The first line of the file corresponds with the header and must begin with ``#IID``.
+The names of each of your phenotypes belong in the subbsequent columns of the header.
+
+See `tests/data/simple.pheno </~https://github.com/gymrek-lab/haptools/blob/main/tests/data/simple.pheno>`_ for an example of a phenotype file:
+
+.. include:: ../../tests/data/simple.pheno
+   :literal:
+
+Covariate file format
+---------------------
+Covariates follow the same format as phenotypes.
+
+See `tests/data/simple.covar </~https://github.com/gymrek-lab/haptools/blob/main/tests/data/simple.covar>`_ for an example of a covariate file:
+
+.. include:: ../../tests/data/simple.covar
+   :literal:

docs/formats/sample_info.rst

@@ -0,0 +1,17 @@
+.. _formats-sample_info:
+
+
+Sample Info
+===========
+
+1000 Genomes sample_info file format
+------------------------------------
+Within the subcommand ``haptools simgenotype`` we use a file to map samples in the 
+reference to their population listed in the model file. This code is expected to work
+out of the box with 1000 genomes data and we have pre-computed this mapping file as 
+well as given the command to how to compute it if you desire another as well.
+
+``cut -f 1,4 igsr-1000\ genomes\ on\ grch38.tsv | sed '1d' | sed -e 's/ /\t/g' > 1000genomes_sampleinfo.tsv``
+
+See ``example-files/1000genomes_sampleinfo.tsv`` for an example of the 1000genomes 
+GRCh38 samples mapped to their subpopulations.

docs/index.rst

@@ -9,8 +9,10 @@
    :hidden:
    :maxdepth: 1
 
-   formats/inputs.rst
+   formats/genotypes.rst
    formats/haplotypes.rst
+   formats/phenotypes.rst
+   formats/sample_info.rst
 
 .. toctree::
    :caption: Commands
@@ -20,6 +22,7 @@
 
    commands/simgenotype.rst
    commands/simphenotype.rst
+   commands/karyogram.rst
    commands/transform.rst
 
 .. toctree::