Merge pull request #164 from sbslee/0.26.0-dev

0.26.0 dev

Author: sbslee

.readthedocs.yml

@@ -9,7 +9,7 @@ version: 2
 build:
   os: ubuntu-22.04
   tools:
-    python: "3.7"
+    python: "3.10"
 
 # Build documentation in the docs/ directory with Sphinx
 sphinx:

CHANGELOG.rst

@@ -1,6 +1,13 @@
 Changelog
 *********
 
+0.26.0 (2026-01-04)
+-------------------
+* :issue:`135`: Fix GRCh38 coordinate of a CYP17A1 variant (thanks `@NeiH4207 <https://github.com/NeiH4207>`__).
+* Add the ``fileformat`` field to the VCF header when writing VCF files.
+* :issue:`150`: Fix bug in :meth:`api.utils.estimate_phase_beagle` method when overlapping samples exist between the input VCF and reference panel, but the window contains only a single position (thanks `@toddknutson <https://github.com/toddknutson>`__).
+* Add new gene MT-RNR1.
+
 0.25.0 (2024-06-16)
 -------------------
 

README.rst

@@ -52,7 +52,7 @@ glossary.html#structural-variation-sv>`__ using a machine learning-based
 approach. Finally, note that PyPGx is compatible with both of the Genome
 Reference Consortium Human (GRCh) builds, GRCh37 (hg19) and GRCh38 (hg38).
 
-There are currently 87 pharmacogenes in PyPGx:
+There are currently 88 pharmacogenes in PyPGx:
 
 .. list-table::
 
@@ -114,38 +114,38 @@ There are currently 87 pharmacogenes in PyPGx:
    * - IFNL3
      - ITGB3
      - ITPA
+     - MT-RNR1
      - MTHFR
-     - NAT1
-   * - NAT2
+   * - NAT1
+     - NAT2
      - NUDT15
      - OPRK1
      - OPRM1
-     - POR
-   * - PTGIS
+   * - POR
+     - PTGIS
      - RARG
      - RYR1
      - SLC6A4
-     - SLC15A2
-   * - SLC22A2
+   * - SLC15A2
+     - SLC22A2
      - SLC28A3
      - SLC47A2
      - SLCO1B1
-     - SLCO1B3
-   * - SLCO2B1
+   * - SLCO1B3
+     - SLCO2B1
      - SULT1A1
      - TBXAS1
      - TPMT
-     - UGT1A1
-   * - UGT1A4
+   * - UGT1A1
+     - UGT1A4
      - UGT1A6
      - UGT2B7
      - UGT2B15
-     - UGT2B17
-   * - VKORC1
+   * - UGT2B17
+     - VKORC1
      - XPC
      -
      -
-     -
 
 Your contributions (e.g. feature ideas, pull requests) are most welcome.
 

docs/create.py

@@ -79,7 +79,7 @@
 approach. Finally, note that PyPGx is compatible with both of the Genome
 Reference Consortium Human (GRCh) builds, GRCh37 (hg19) and GRCh38 (hg38).
 
-There are currently 87 pharmacogenes in PyPGx:
+There are currently 88 pharmacogenes in PyPGx:
 
 .. list-table::
 
@@ -141,38 +141,38 @@
    * - IFNL3
      - ITGB3
      - ITPA
+     - MT-RNR1
      - MTHFR
-     - NAT1
-   * - NAT2
+   * - NAT1
+     - NAT2
      - NUDT15
      - OPRK1
      - OPRM1
-     - POR
-   * - PTGIS
+   * - POR
+     - PTGIS
      - RARG
      - RYR1
      - SLC6A4
-     - SLC15A2
-   * - SLC22A2
+   * - SLC15A2
+     - SLC22A2
      - SLC28A3
      - SLC47A2
      - SLCO1B1
-     - SLCO1B3
-   * - SLCO2B1
+   * - SLCO1B3
+     - SLCO2B1
      - SULT1A1
      - TBXAS1
      - TPMT
-     - UGT1A1
-   * - UGT1A4
+   * - UGT1A1
+     - UGT1A4
      - UGT1A6
      - UGT2B7
      - UGT2B15
-     - UGT2B17
-   * - VKORC1
+   * - UGT2B17
+     - VKORC1
      - XPC
      -
      -
-     -
 
 Your contributions (e.g. feature ideas, pull requests) are most welcome.
 
@@ -577,18 +577,18 @@
 
 - Method 1. Simple diplotype-phenotype mapping: This method directly uses the
   diplotype-phenotype mapping as defined by CPIC or PharmGKB. Using the
-  CYP2B6 gene as an example, the diplotypes \*6/\*6, \*1/\*29, \*1/\*2,
-  \*1/\*4, and \*4/\*4 correspond to Poor Metabolizer, Intermediate
+  CYP2B6 gene as an example, the diplotypes \\*6/\\*6, \\*1/\\*29, \\*1/\\*2,
+  \\*1/\\*4, and \\*4/\\*4 correspond to Poor Metabolizer, Intermediate
   Metabolizer, Normal Metabolizer, Rapid Metabolizer, and Ultrarapid
   Metabolizer.
 - Method 2. Summation of haplotype activity scores: This method uses a
   standard unit of enzyme activity known as an activity score. Using the
-  CYP2D6 gene as an example, the fully functional reference \*1 allele is
-  assigned a value of 1, decreased-function alleles such as \*9 and \*17
-  receive a value of 0.5, and nonfunctional alleles including \*4 and \*5
+  CYP2D6 gene as an example, the fully functional reference \\*1 allele is
+  assigned a value of 1, decreased-function alleles such as \\*9 and \\*17
+  receive a value of 0.5, and nonfunctional alleles including \\*4 and \\*5
   have a value of 0. The sum of values assigned to both alleles constitutes
-  the activity score of a diplotype. Consequently, subjects with \*1/\*1,
-  \*1/\*4, and \*4/\*5 diplotypes have an activity score of 2 (Normal
+  the activity score of a diplotype. Consequently, subjects with \\*1/\\*1,
+  \\*1/\\*4, and \\*4/\\*5 diplotypes have an activity score of 2 (Normal
   Metabolizer), 1 (Intermediate Metabolizer), and 0 (Poor Metabolizer),
   respectively.
 
@@ -718,45 +718,45 @@
      - VariantData
      - CNV
    * - NA11839
-     - \*1/\*2
+     - \\*1/\\*2
      - Normal Metabolizer
-     - \*1;
-     - \*2;
+     - \\*1;
+     - \\*2;
      - ;
-     - \*1:22-42522613-G-C,22-42523943-A-G:0.5,0.488;\*2:default
+     - \\*1:22-42522613-G-C,22-42523943-A-G:0.5,0.488;\\*2:default
      - Normal
    * - NA12006
-     - \*4/\*41
+     - \\*4/\\*41
      - Intermediate Metabolizer
-     - \*41;\*2;
-     - \*4;\*10;\*2;
-     - \*69;
-     - \*69:22-42526694-G-A,22-42523805-C-T:0.5,0.551;\*4:22-42524947-C-T:0.444;\*10:22-42523943-A-G,22-42526694-G-A:0.55,0.5;\*41:22-42523805-C-T:0.551;\*2:default;
+     - \\*41;\\*2;
+     - \\*4;\\*10;\\*2;
+     - \\*69;
+     - \\*69:22-42526694-G-A,22-42523805-C-T:0.5,0.551;\\*4:22-42524947-C-T:0.444;\\*10:22-42523943-A-G,22-42526694-G-A:0.55,0.5;\\*41:22-42523805-C-T:0.551;\\*2:default;
      - Normal
    * - HG00276
-     - \*4/\*5
+     - \\*4/\\*5
      - Poor Metabolizer
-     - \*4;\*10;\*74;\*2;
-     - \*10;\*74;\*2;
+     - \\*4;\\*10;\\*74;\\*2;
+     - \\*10;\\*74;\\*2;
      - ;
-     - \*4:22-42524947-C-T:0.913;\*10:22-42523943-A-G,22-42526694-G-A:1.0,1.0;\*74:22-42525821-G-T:1.0;\*2:default;
+     - \\*4:22-42524947-C-T:0.913;\\*10:22-42523943-A-G,22-42526694-G-A:1.0,1.0;\\*74:22-42525821-G-T:1.0;\\*2:default;
      - WholeDel1
    * - NA19207
-     - \*2x2/\*10
+     - \\*2x2/\\*10
      - Normal Metabolizer
-     - \*10;\*2;
-     - \*2;
+     - \\*10;\\*2;
+     - \\*2;
      - ;
-     - \*10:22-42523943-A-G,22-42526694-G-A:0.361,0.25;\*2:default;
+     - \\*10:22-42523943-A-G,22-42526694-G-A:0.361,0.25;\\*2:default;
      - WholeDup1
 
 This list explains each of the columns in the example results.
 
 - **Genotype**: Diplotype call. When there is no SV this simply combines the two top-ranked star alleles from **Haplotype1** and **Haplotype2** with the delimiter '/'. In the presence of SV the final diplotype is determined using one of the genotypers in the ``pypgx.api.genotype`` module (e.g. `CYP2D6Genotyper <https://pypgx.readthedocs.io/en/latest/api.html#pypgx.api.genotype.CYP2D6Genotyper>`__).
 - **Phenotype**: Phenotype call.
-- **Haplotype1**, **Haplotype2**: List of candidate star alleles for each haplotype. For example, if a given haplotype contains three variants ``22-42523943-A-G``, ``22-42524947-C-T``, and ``22-42526694-G-A``, then it will get assigned ``*4;*10;`` because the haplotype pattern can fit both \*4 (``22-42524947-C-T``) and \*10 (``22-42523943-A-G`` and ``22-42526694-G-A``). Note that \*4 comes first before \*10 because it has higher priority for reporting purposes (see the ``pypgx.sort_alleles`` `method <https://pypgx.readthedocs.io/en/latest/api.html#pypgx.api.core.sort_alleles>`__ for detailed implementation).
-- **AlternativePhase**: List of star alleles that could be missed due to potentially incorrect statistical phasing. For example, let's assume that statistical phasing has put ``22-42526694-G-A`` for **Haplotype1** and ``22-42523805-C-T`` for **Haplotype2**. Even though the two variants are in trans orientation, PyPGx will also consider alternative phase in case the two variants are actually in cis orientation, resulting in ``*69;`` as **AlternativePhase** because \*69 is defined by ``22-42526694-G-A`` and ``22-42523805-C-T``.
-- **VariantData**: Information for SNVs/indels used to define observed star alleles, including allele fraction which is important for allelic decomposition after identifying CNV (e.g. the sample NA19207). In some situations, there will not be any variants for a given star allele because the allele itself is "default" allele for the selected reference assembly (e.g. GRCh37 has \*2 as default while GRCh38 has \*1).
+- **Haplotype1**, **Haplotype2**: List of candidate star alleles for each haplotype. For example, if a given haplotype contains three variants ``22-42523943-A-G``, ``22-42524947-C-T``, and ``22-42526694-G-A``, then it will get assigned ``*4;*10;`` because the haplotype pattern can fit both \\*4 (``22-42524947-C-T``) and \\*10 (``22-42523943-A-G`` and ``22-42526694-G-A``). Note that \\*4 comes first before \\*10 because it has higher priority for reporting purposes (see the ``pypgx.sort_alleles`` `method <https://pypgx.readthedocs.io/en/latest/api.html#pypgx.api.core.sort_alleles>`__ for detailed implementation).
+- **AlternativePhase**: List of star alleles that could be missed due to potentially incorrect statistical phasing. For example, let's assume that statistical phasing has put ``22-42526694-G-A`` for **Haplotype1** and ``22-42523805-C-T`` for **Haplotype2**. Even though the two variants are in trans orientation, PyPGx will also consider alternative phase in case the two variants are actually in cis orientation, resulting in ``*69;`` as **AlternativePhase** because \\*69 is defined by ``22-42526694-G-A`` and ``22-42523805-C-T``.
+- **VariantData**: Information for SNVs/indels used to define observed star alleles, including allele fraction which is important for allelic decomposition after identifying CNV (e.g. the sample NA19207). In some situations, there will not be any variants for a given star allele because the allele itself is "default" allele for the selected reference assembly (e.g. GRCh37 has \\*2 as default while GRCh38 has \\*1).
 - **CNV**: Structural variation call. See the `Structural variation detection <https://pypgx.readthedocs.io/en/latest/readme.html#structural-variation-detection>`__ section for more details.
 
 Getting help

docs/genes.rst

@@ -2,7 +2,7 @@ Genes
 *****
 
 This page describes gene-specific information. PyPGx currently supports
-genotyping of a total of 87 pharmacogenes.
+genotyping of a total of 88 pharmacogenes.
 
 In order to provide the most accurate information, this page borrows heavily
 from the works of the :ref:`glossary:Clinical Pharmacogenetics Implementation
@@ -616,6 +616,16 @@ Below is a summary table:
      - `chr20:3187143-3207506 <https://genome.ucsc.edu/cgi-bin/hgTracks?db=hg19&lastVirtModeType=default&lastVirtModeExtraState=&virtModeType=default&virtMode=0&nonVirtPosition=&position=chr20%3A3187143%2D3207506&hgsid=2144291390_4bmNOa6wq7Mk8crO6QNojSla7rfr>`__
      - `chr20:3206497-3226860 <https://genome.ucsc.edu/cgi-bin/hgTracks?db=hg38&lastVirtModeType=default&lastVirtModeExtraState=&virtModeType=default&virtMode=0&nonVirtPosition=&position=chr20%3A3206497%2D3226860&hgsid=2144291444_9qcI5Mpt1r7ap1R1ceKlp0aFJjS3>`__
      - NM_033453.4 was used as the main transcript.
+   * - :ref:`genes:MT-RNR1`
+     - ✅
+     -
+     - ✅
+     -
+     -
+     - Disease
+     - `chrM:410-1840 <https://genome.ucsc.edu/cgi-bin/hgTracks?db=hg19&lastVirtModeType=default&lastVirtModeExtraState=&virtModeType=default&virtMode=0&nonVirtPosition=&position=chrMT%3A410%2D1840&hgsid=3332880867_aVXAZCwJTtLSGNhbga0xleh3qxVY>`__
+     - `chrM:410-1840 <https://genome.ucsc.edu/cgi-bin/hgTracks?db=hg38&lastVirtModeType=default&lastVirtModeExtraState=&virtModeType=default&virtMode=0&nonVirtPosition=&position=chrM%3A410%2D1840&hgsid=3332880867_aVXAZCwJTtLSGNhbga0xleh3qxVY>`__
+     -
    * - MTHFR
      - ✅
      -
@@ -3296,6 +3306,30 @@ Resources for IFNL3
 - `PharmGKB: Annotation of CPIC Guideline for peginterferon alfa-2a,peginterferon alfa-2b,ribavirin and IFNL3 <https://www.pharmgkb.org/guidelineAnnotation/PA166110235>`__
 - `CPIC® Guideline for PEG Interferon-Alpha-Based Regimens and IFNL3 <https://cpicpgx.org/guidelines/guideline-for-peg-interferon-alpha-based-regimens-and-ifnl3/>`__
 
+MT-RNR1
+=======
+
+Phenotype summary for MT-RNR1
+-----------------------------
+
+Diplotype-phenotype mapping is used for phenotype prediction.
+
+.. list-table::
+   :header-rows: 1
+
+   * - Phenotype
+     - Example
+     - Priority
+   * - Normal Risk of Aminoglycoside-Induced Hearing Loss
+     - Reference/Reference
+     - None
+   * - Increased Risk of Aminoglycoside-Induced Hearing Loss
+     - m.1494C>T/m.1494C>T
+     - None
+   * - Uncertain Risk of Aminoglycoside-Induced Hearing Loss
+     - m.663A>G/m.663A>G
+     - None
+
 NUDT15
 ======
 

docs/requirements.txt

@@ -1,6 +1,6 @@
 autodocsumm
 fuc
 scikit-learn
-sphinx==4.1.2
-sphinx_rtd_theme==0.5.2
+sphinx
+sphinx_rtd_theme
 sphinx_issues

docs/tutorials.rst

@@ -347,6 +347,66 @@ Now you have successfully genotyped CYP3A5 as well!
     gene requires SV analysis. In other words, users can provide the same
     input files for all target genes.
 
+Genotyping with CRAM files
+--------------------------
+
+PyPGx also supports CRAM files as input.
+
+.. code-block:: text
+
+    $ mkdir grch37-cram
+    $ wget -P grch37-cram https://storage.googleapis.com/sbslee-bucket/pypgx/getrm-wgs-tutorial/grch37-cram/HG00276_PyPGx.sorted.markdup.recal.cram.crai
+    $ wget -P grch37-cram https://storage.googleapis.com/sbslee-bucket/pypgx/getrm-wgs-tutorial/grch37-cram/HG00276_PyPGx.sorted.markdup.recal.cram
+    $ wget -P grch37-cram https://storage.googleapis.com/sbslee-bucket/pypgx/getrm-wgs-tutorial/grch37-cram/HG00436_PyPGx.sorted.markdup.recal.cram.crai
+    $ wget -P grch37-cram https://storage.googleapis.com/sbslee-bucket/pypgx/getrm-wgs-tutorial/grch37-cram/HG00436_PyPGx.sorted.markdup.recal.cram
+
+Similar to before, we will need the reference FASTA file used to create the CRAM files:
+
+.. code-block:: text
+
+    $ wget https://storage.googleapis.com/sbslee-bucket/ref/grch37/genome.fa
+    $ wget https://storage.googleapis.com/sbslee-bucket/ref/grch37/genome.fa.fai
+
+We can create input files as usual:
+
+.. code-block:: text
+
+    $ pypgx create-input-vcf \
+    cram-variants.vcf.gz \
+    genome.fa \
+    grch37-cram/*.cram
+
+.. code-block:: text
+
+    $ pypgx prepare-depth-of-coverage \
+    cram-depth-of-coverage.zip \
+    grch37-cram/*.cram
+
+.. code-block:: text
+
+    $ pypgx compute-control-statistics \
+    VDR \
+    cram-control-statistics-VDR.zip \
+    grch37-cram/*.cram
+
+Finally, we can run the NGS pipeline:
+
+.. code-block:: text
+
+    $ pypgx run-ngs-pipeline \
+    CYP2D6 \
+    cram-CYP2D6-pipeline \
+    --variants cram-variants.vcf.gz \
+    --depth-of-coverage cram-depth-of-coverage.zip \
+    --control-statistics cram-control-statistics-VDR.zip
+
+.. code-block:: text
+
+    $ pypgx print-data cram-CYP2D6/results.zip
+        Genotype	Phenotype	Haplotype1	Haplotype2	AlternativePhase	VariantData	CNV
+    HG00276_PyPGx	*4/*5	Poor Metabolizer	*4;*10;*74;*2;	*4;*10;*74;*2;	;	*4:22-42524947-C-T:0.87;*10:22-42523943-A-G,22-42526694-G-A:1.0,1.0;*74:22-42525821-G-T:1.0;*2:default;	WholeDel1
+    HG00436_PyPGx	*2x2/*71	Indeterminate	*71;*1;	*2;	;	*71:22-42526669-C-T:0.433;*1:22-42523943-A-G,22-42522613-G-C:0.353,0.462;*2:default;	WholeDup1
+
 Genotyping with GRCh38 data
 ---------------------------
 

pypgx/api/core.py

@@ -1498,17 +1498,17 @@ def sort_alleles(
     impact protein coding (e.g. two misssense variants > one missense variant
     plus one intron variant), and 4. reference allele status (e.g.
     non-reference allele with two SNVs > reference allele with two SNVs such
-    that CYP2D6\*46 > CYP2D6\*1 in GRCh37). Note that the priority of allele
+    that CYP2D6\\*46 > CYP2D6\\*1 in GRCh37). Note that the priority of allele
     function decreases in the following order: 'No Function', 'Decreased
     Function', 'Possible Decreased Function', 'Increased Function', 'Possible
     Increased Function', 'Uncertain Function', 'Unknown Function', 'Normal
     Function'.
 
     When ``by='name'`` the method will report alleles with a smaller
-    number first. This means, for example, '\*4' will come before '\*10'
+    number first. This means, for example, '\\*4' will come before '\\*10'
     whereas lexicographic sorting would produce the opposite result. This is
-    particularly useful when forming a diplotype (e.g. '\*4/\*10' vs.
-    '\*10/\*4').
+    particularly useful when forming a diplotype (e.g. '\\*4/\\*10' vs.
+    '\\*10/\\*4').
 
     Parameters
     ----------