feat: a new GenotypesPLINKTR class for reading TRs from PGEN files

docs/api/data.rst

@@ -168,7 +168,7 @@ All of the other methods in the :class:`Genotypes` class are inherited, but the
 
 GenotypesTR
 ++++++++++++
-The :class:`GenotypesTR` class *extends* :class:`Genotypes` class. The :class:`GenotypesTR` class follows the same structure of :class:`GenotypesVCF`, but can now load repeat count of tandem repeats as the alleles.
+The :class:`GenotypesTR` class *extends* the :class:`Genotypes` class. The :class:`GenotypesTR` class follows the same structure of :class:`GenotypesVCF`, but can now load repeat counts of tandem repeats as the alleles.
 
 All of the other methods in the :class:`Genotypes` class are inherited, but the :class:`GenotypesTR` class' ``load()`` function is unique to loading tandem repeat variants.
 
@@ -178,6 +178,11 @@ All of the other methods in the :class:`Genotypes` class are inherited, but the
 	# make the first sample have 4 and 7 repeats for the alleles of the fourth variant
 	genotypes.data[0, 3] = (4, 7)
 
+The following methods from the :class:`Genotypes` class are disabled, however.
+
+1. ``check_biallelic``
+2. ``check_maf``
+
 .. _api-data-genotypestr:
 
 GenotypesPLINK
@@ -239,6 +244,24 @@ A large ``chunk_size`` is more likely to result in memory over-use while a small
 	genotypes = data.GenotypesPLINK('tests/data/simple.pgen', chunk_size=500)
 	genotypes.read()
 
+GenotypesPLINKTR
+++++++++++++++++
+The :class:`GenotypesPLINKTR`` class extends the :class:`GenotypesPLINK` class to support loading tandem repeat variants.
+The :class:`GenotypesPLINKTR` class works similarly to :class:`GenotypesTR` by filling the ``data`` property with repeat counts for each allele.
+
+The following methods from the :class:`GenotypesPLINK` class are disabled, however.
+
+1. ``write``
+2. ``check_maf``
+3. ``write_variants``
+4. ``check_biallelic``
+
+The :class:`GenotypesPLINKTR` uses INFO fields from the PVAR file to determine the repeat unit and the number of repeats for each allele. To ensure your PVAR file contains the necessary information, use the following command when converting from VCF.
+
+.. code-block:: bash
+
+	plink2 --vcf-half-call m --make-pgen 'pvar-cols=vcfheader,qual,filter,info' --vcf input.vcf --make-pgen --out output
+
 haplotypes.py
 ~~~~~~~~~~~~~
 Overview

docs/formats/genotypes.rst

@@ -9,14 +9,14 @@ Genotypes
 	The time required to load various genotype file formats.
 
 VCF/BCF
--------
+~~~~~~~
 
 Genotype files must be specified as VCF or BCF files. They can be bgzip-compressed.
 
 .. _formats-genotypesplink:
 
 PLINK2 PGEN
------------
+~~~~~~~~~~~
 
 There is also experimental support for `PLINK2 PGEN <https://github.com/chrchang/plink-ng/blob/master/pgen_spec/pgen_spec.pdf>`_ files in some commands. These files can be loaded and created much more quickly than VCFs, so we highly recommend using them if you're working with large datasets. See the documentation for the :class:`GenotypesPLINK` class in :ref:`the API docs <api-data-genotypesplink>` for more information.
 
@@ -29,9 +29,16 @@ If you run out memory when using PGEN files, consider reading/writing variants f
 
 		pip install haptools[files]
 
-.. warning::
-	At the moment, only biallelic SNPs can be encoded in PGEN files because of limitations in the ``Pgenlib`` python library. It doesn't properly support multiallelic variants yet (`source <https://github.com/chrchang/plink-ng/blob/c4b8d4361de74c58f0cc11361062eca4f34210d3/2.0/Python/python_api.txt#L88-L89>`_). To ensure your PGEN files only contain SNPs, we recommend use the following command to convert from VCF to PGEN.
+Converting from VCF to PGEN
+---------------------------
+To convert a VCF containing only biallelic SNPs to PGEN, use the following command.
 
-	.. code-block:: bash
+.. code-block:: bash
+
+	plink2 --snps-only 'just-acgt' --max-alleles 2 --vcf input.vcf --make-pgen --out output
+
+To convert a VCF containing tandem repeats to PGEN, use this command, instead.
+
+.. code-block:: bash
 
-		plink2 --snps-only 'just-acgt' --max-alleles 2 --vcf input.vcf --make-pgen --out output
+	plink2 --vcf-half-call m --make-pgen 'pvar-cols=vcfheader,qual,filter,info' --vcf input.vcf --make-pgen --out output

docs/project_info/contributing.rst

@@ -70,7 +70,7 @@ Follow these steps to set up a development environment.
 
     .. code-block:: bash
 
-        conda create -n haptools-dev -c conda-forge 'poetry>=1.1.15' 'python=3.7' 'gxx_linux-64'
+        conda create -n haptools-dev -c conda-forge 'poetry>=1.1.15' 'python=3.7'
 2. Activate the environment
 
     .. code-block:: bash

haptools/clump.py

@@ -1,13 +1,12 @@
 #!/usr/bin/env python
 
 # To test: ./clumpSTR.py --summstats-snps tests/eur_gwas_pvalue_chr19.LDL.glm.linear --clump-snp-field ID --clump-field p-value --clump-chrom-field CHROM --clump-pos-field position --clump-p1 0.2 --out test.clump
-import sys
 import math
 from logging import Logger, getLogger
 
 import numpy as np
 
-from .data import Genotypes, GenotypesVCF, GenotypesTR
+from .data import Genotypes, GenotypesVCF, GenotypesTR, GenotypesPLINKTR
 
 
 class Variant:
@@ -557,15 +556,17 @@ def clumpstr(
     strgts = None
     gts = None
     if gts_snps:
+        log.debug("Loading SNP Genotypes.")
         if str(gts_snps).endswith("pgen"):
-            log.debug("Loading SNP Genotypes.")
             snpgts = GenotypesPLINK.load(gts_snps)
         else:
-            log.debug("Loading SNP Genotypes.")
             snpgts = GenotypesVCF.load(gts_snps)
     if gts_strs:
         log.debug("Loading STR Genotypes.")
-        strgts = GenotypesTR.load(gts_strs)
+        if str(gts_strs).endswith("pgen"):
+            strgts = GenotypesPLINKTR.load(gts_strs)
+        else:
+            strgts = GenotypesTR.load(gts_strs)
 
     if gts_snps and gts_strs:
         log.debug("Calculating set of overlapping samples between STRs and SNPs.")

haptools/data/__init__.py

@@ -3,4 +3,10 @@
 from .covariates import Covariates
 from .breakpoints import Breakpoints, HapBlock
 from .haplotypes import Extra, Repeat, Variant, Haplotype, Haplotypes
-from .genotypes import Genotypes, GenotypesVCF, GenotypesTR, GenotypesPLINK
+from .genotypes import (
+    Genotypes,
+    GenotypesVCF,
+    GenotypesTR,
+    GenotypesPLINK,
+    GenotypesPLINKTR,
+)