Merge pull request #43 from sbslee/0.11.0-dev

0.11.0 dev

Author: sbslee

CHANGELOG.rst

@@ -1,6 +1,23 @@
 Changelog
 *********
 
+0.11.0 (2022-01-01)
+-------------------
+
+* Add CNV caller for CYP4F2 and SULT1A1.
+* Fix minor bug in :command:`compute-copy-number` command.
+* Update :command:`plot-cn-af` command to check input files more rigorously.
+* Improve CNV caller for CYP2A6, CYP2D6, SLC22A2, and SULT1A1.
+* Add new method :meth:`sdk.utils.add_cn_samples` method.
+* Update :command:`compare-genotypes` command to output CNV comparison results as well.
+* Update :command:`estimate-phase-beagle` command. From now on, the 'chr' prefix in contig names (e.g. 'chr1' vs. '1') will be automatically added or removed as necessary to match the reference VCF’s contig names.
+* Add index files for 1KGP reference haplotype panels.
+* Add new argument ``--panel`` to :command:`run-chip-pipeline` command.
+* Remove 1KGP reference haplotype panels for GSTT1 and UGT2B17 because these genes only have star alleles defined with SV.
+* Change 1KGP reference haplotype panels for GRCh38. Previously, PyPGx was using the panels from `Lowy-Gallego et al., 2019 <https://wellcomeopenresearch.org/articles/4-50>`__ where the authors had aligned sequence reads against the full GRCh38 reference, including ALT contigs, decoy, and EBV/IMGT/HLA sequences. This resulted in poor phasing/imputation performance for highly polymorphic PGx genes (e.g. CYP2D6) presumably because the panels were missing haplotype information for lots of SNVs/indels as sequence reads with those variants were mapped to ALT contigs; however, the panels were still the best option at the time (definitely better than lifting over GRCh37 panels). Fortunately, `Byrska-Bishop et al., 2021 <https://www.biorxiv.org/content/10.1101/2021.02.06.430068v2>`__ from New York Genome Center has recently published a new set of GRCh38 panels which apparently has less of this problem despite still having sequence reads aligned in the presence of ALT contigs, etc. When empirically tested, these panels showed a significant increase in phasing/imputation performance. Therefore, from now on, PyPGx will use these panels for GRCh38 data.
+* Update GRCh38 variant information for following alleles: CYP2D6\*35, CYP2D6\*45, CYP2D6\*46.
+* Update gene region for SLC22A2 to match GRCh37 and GRCh38.
+
 0.10.1 (2021-12-20)
 -------------------
 
@@ -72,7 +89,7 @@ Changelog
 * Fix minor bug in :command:`estimate-phase-beagle` command on not properly exiting the program even though there was an error raised by Beagle.
 * Update :meth:`api.utils.create_consolidated_vcf` method to check synonymous variants as well when performing phase-extension algorithm.
 * Update :command:`run-ngs-pipeline` command to output a warning when user provides CovFrame[DepthOfCoverage] even though target gene does not have any star alleles defined by SVs.
-* Add new argument ``--fontsize`` argument to :command:`plot-bam-copy-number` command.
+* Add new argument ``--fontsize`` to :command:`plot-bam-copy-number` command.
 * Remove ``--ymin`` and ``--ymax`` arguments from :command:`plot-vcf-allele-fraction` command.
 * Update ``--ymin`` and ``--ymax`` arguments of :command:`plot-bam-copy-number` command to have a default value.
 * Add new command :command:`plot-cn-af`.

README.rst

@@ -158,6 +158,46 @@ you can access a development branch with the ``git checkout`` command. When
 you do this, please make sure your environment already has all the
 dependencies installed.
 
+Structural variation detection
+==============================
+
+Many pharmacogenes are known to have `structural variation (SV)
+<https://pypgx.readthedocs.io/en/latest/glossary.html#structural-variation-
+sv>`__ such as gene deletions, duplications, and hybrids. You can visit the
+`Genes <https://pypgx.readthedocs.io/en/latest/genes.html>`__ page to see the
+list of genes with SV.
+
+Some of the SV events can be quite challenging to detect accurately with
+next-generation sequencing (NGS) data due to misalignment of sequence reads
+caused by sequence homology with other gene family members (e.g. CYP2D6 and
+CYP2D7). PyPGx attempts to address this issue by training a `support vector
+machine (SVM) <https://scikit-learn.org/stable/modules/generated/sk
+learn.svm.SVC.html>`__-based multiclass classifier using the `one-vs-rest
+strategy <https://scikit-learn.org/stable/modules/generated/sklearn.multi
+class.OneVsRestClassifier.html>`__ for each gene for each GRCh build. Each
+classifier is trained using copy number profiles of real NGS samples as well
+as simulated ones.
+
+You can plot copy number profile and allele fraction profile with PyPGX to
+visually inspect SV calls. Below are CYP2D6 examples:
+
+.. list-table::
+   :header-rows: 1
+   :widths: 20 80
+
+   * - SV Name
+     - Profile
+   * - Normal
+     - .. image:: https://raw.githubusercontent.com/sbslee/pypgx-data/main/dpsv/GRCh37-CYP2D6-8.png
+   * - DeletionHet
+     - .. image:: https://raw.githubusercontent.com/sbslee/pypgx-data/main/dpsv/GRCh37-CYP2D6-1.png
+   * - Duplication
+     - .. image:: https://raw.githubusercontent.com/sbslee/pypgx-data/main/dpsv/GRCh37-CYP2D6-2.png
+   * - Tandem3
+     - .. image:: https://raw.githubusercontent.com/sbslee/pypgx-data/main/dpsv/GRCh37-CYP2D6-9.png
+   * - Tandem2C
+     - .. image:: https://raw.githubusercontent.com/sbslee/pypgx-data/main/dpsv/GRCh37-CYP2D6-7.png
+
 GRCh37 vs. GRCh38
 =================
 
@@ -170,6 +210,15 @@ do this). The good news is, PyPGx supports both of the builds!
 In many of the PyPGx actions, you can simply indicate which human genome
 build to use. For example, you can use ``assembly`` for the API and
 ``--assembly`` for the CLI. **Note that GRCh37 will always be the default.**
+Below is an example of using the API:
+
+.. code:: python3
+
+    >>> import pypgx
+    >>> pypgx.list_variants('CYP2D6', alleles=['*4'], assembly='GRCh37')
+    ['22-42524947-C-T']
+    >>> pypgx.list_variants('CYP2D6', alleles=['*4'], assembly='GRCh38')
+    ['22-42128945-C-T']
 
 However, there is one important caveat to consider if your sequencing data is
 GRCh38. That is, sequence reads must be aligned only to the main contigs
@@ -314,6 +363,29 @@ Please visit the `Genes <https://pypgx.readthedocs.io/en/latest/
 genes.html>`__ page to see the list of genes with a genotype-phenotype
 table and each of their prediction method.
 
+To perform phenotype prediction with the API, you can use the
+``pypgx.predict_phenotype`` method:
+
+.. code:: python3
+
+    >>> import pypgx
+    >>> pypgx.predict_phenotype('CYP2D6', '*4', '*5')   # Both alleles have no function
+    'Poor Metabolizer'
+    >>> pypgx.predict_phenotype('CYP2D6', '*5', '*4')   # The order of alleles does not matter
+    'Poor Metabolizer'
+    >>> pypgx.predict_phenotype('CYP2D6', '*1', '*22')  # *22 has uncertain function
+    'Indeterminate'
+    >>> pypgx.predict_phenotype('CYP2D6', '*1', '*1x2') # Gene duplication
+    'Ultrarapid Metabolizer'
+
+To perform phenotype prediction with the CLI, you can use the
+``call-phenotypes`` command. It takes a ``SampleTable[Genotypes]`` file as
+input and outputs a ``SampleTable[Phenotypes]`` file:
+
+.. code-block:: text
+
+   $ pypgx call-phenotypes genotypes.zip phenotypes.zip
+
 Getting help
 ============
 
@@ -333,7 +405,7 @@ For getting help on the CLI:
        call-genotypes      Call genotypes for the target gene.
        call-phenotypes     Call phenotypes for the target gene.
        combine-results     Combine various results for the target gene.
-       compare-genotypes   Calculate concordance rate between two genotype results.
+       compare-genotypes   Calculate concordance between two genotype results.
        compute-control-statistics
                            Compute summary statistics for the control gene from BAM files.
        compute-copy-number
@@ -385,20 +457,23 @@ Below is the list of submodules available in the API:
 - **plot** : The plot submodule is used to plot various kinds of profiles such as read depth, copy number, and allele fraction.
 - **utils** : The utils submodule contains main actions of PyPGx.
 
-For getting help on a specific submodule (e.g. utils):
+For getting help on a specific submodule (e.g. ``utils``):
 
 .. code:: python3
 
    >>> from pypgx.api import utils
    >>> help(utils)
 
-For getting help on a specific method (e.g. predict_phenotype):
+For getting help on a specific method (e.g. ``predict_phenotype``):
 
 .. code:: python3
 
    >>> import pypgx
    >>> help(pypgx.predict_phenotype)
 
+In Jupyter Notebook and Lab, you can see the documentation for a python
+function by hitting ``SHIFT + TAB``. Hit it twice to expand the view.
+
 CLI examples
 ============
 

docs/cli.rst

@@ -23,7 +23,7 @@ For getting help on the CLI:
        call-genotypes      Call genotypes for the target gene.
        call-phenotypes     Call phenotypes for the target gene.
        combine-results     Combine various results for the target gene.
-       compare-genotypes   Calculate concordance rate between two genotype results.
+       compare-genotypes   Calculate concordance between two genotype results.
        compute-control-statistics
                            Compute summary statistics for the control gene from BAM files.
        compute-copy-number
@@ -137,9 +137,10 @@ compare-genotypes
    $ pypgx compare-genotypes -h
    usage: pypgx compare-genotypes [-h] [--verbose] first second
    
-   Calculate concordance rate between two genotype results.
+   Calculate concordance between two genotype results.
    
-   The command will only use samples that appear in both genotype results.
+   Only samples that appear in both genotype results will be used to calculate
+   concordance for genotype calls as well as CNV calls.
    
    Positional arguments:
      first       First archive file with the semantic type 
@@ -149,7 +150,8 @@ compare-genotypes
    
    Optional arguments:
      -h, --help  Show this help message and exit.
-     --verbose   Whether to print the verbose version of output.
+     --verbose   Whether to print the verbose version of output, including 
+                 discordant calls.
 
 compute-control-statistics
 ==========================
@@ -327,6 +329,10 @@ estimate-phase-beagle
    
    Estimate haplotype phase of observed variants with the Beagle program.
    
+   Note that the 'chr' prefix in contig names (e.g. 'chr1' vs. '1') will be
+   automatically added or removed as necessary to match the reference VCF’s
+   contig names.
+   
    Positional arguments:
      imported-variants  Archive file with the semantic type VcfFrame[Imported].
      phased-variants    Archive file with the semantic type VcfFrame[Phased].
@@ -674,7 +680,8 @@ run-chip-pipeline
 .. code-block:: text
 
    $ pypgx run-chip-pipeline -h
-   usage: pypgx run-chip-pipeline [-h] [--assembly TEXT] [--impute] [--force]
+   usage: pypgx run-chip-pipeline [-h] [--assembly TEXT] [--panel PATH]
+                                  [--impute] [--force]
                                   [--samples TEXT [TEXT ...]] [--exclude]
                                   gene output variants
    
@@ -692,6 +699,9 @@ run-chip-pipeline
      -h, --help            Show this help message and exit.
      --assembly TEXT       Reference genome assembly (default: 'GRCh37') (choices: 
                            'GRCh37', 'GRCh38').
+     --panel PATH          VCF file corresponding to a reference haplotype panel 
+                           (compressed or uncompressed). By default, the 1KGP 
+                           panel is used.
      --impute              Perform imputation of missing genotypes.
      --force               Overwrite output directory if it already exists.
      --samples TEXT [TEXT ...]

docs/create.py

@@ -185,6 +185,46 @@
 you do this, please make sure your environment already has all the
 dependencies installed.
 
+Structural variation detection
+==============================
+
+Many pharmacogenes are known to have `structural variation (SV)
+<https://pypgx.readthedocs.io/en/latest/glossary.html#structural-variation-
+sv>`__ such as gene deletions, duplications, and hybrids. You can visit the
+`Genes <https://pypgx.readthedocs.io/en/latest/genes.html>`__ page to see the
+list of genes with SV.
+
+Some of the SV events can be quite challenging to detect accurately with
+next-generation sequencing (NGS) data due to misalignment of sequence reads
+caused by sequence homology with other gene family members (e.g. CYP2D6 and
+CYP2D7). PyPGx attempts to address this issue by training a `support vector
+machine (SVM) <https://scikit-learn.org/stable/modules/generated/sk
+learn.svm.SVC.html>`__-based multiclass classifier using the `one-vs-rest
+strategy <https://scikit-learn.org/stable/modules/generated/sklearn.multi
+class.OneVsRestClassifier.html>`__ for each gene for each GRCh build. Each
+classifier is trained using copy number profiles of real NGS samples as well
+as simulated ones.
+
+You can plot copy number profile and allele fraction profile with PyPGX to
+visually inspect SV calls. Below are CYP2D6 examples:
+
+.. list-table::
+   :header-rows: 1
+   :widths: 20 80
+
+   * - SV Name
+     - Profile
+   * - Normal
+     - .. image:: https://raw.githubusercontent.com/sbslee/pypgx-data/main/dpsv/GRCh37-CYP2D6-8.png
+   * - DeletionHet
+     - .. image:: https://raw.githubusercontent.com/sbslee/pypgx-data/main/dpsv/GRCh37-CYP2D6-1.png
+   * - Duplication
+     - .. image:: https://raw.githubusercontent.com/sbslee/pypgx-data/main/dpsv/GRCh37-CYP2D6-2.png
+   * - Tandem3
+     - .. image:: https://raw.githubusercontent.com/sbslee/pypgx-data/main/dpsv/GRCh37-CYP2D6-9.png
+   * - Tandem2C
+     - .. image:: https://raw.githubusercontent.com/sbslee/pypgx-data/main/dpsv/GRCh37-CYP2D6-7.png
+
 GRCh37 vs. GRCh38
 =================
 
@@ -197,6 +237,15 @@
 In many of the PyPGx actions, you can simply indicate which human genome
 build to use. For example, you can use ``assembly`` for the API and
 ``--assembly`` for the CLI. **Note that GRCh37 will always be the default.**
+Below is an example of using the API:
+
+.. code:: python3
+
+    >>> import pypgx
+    >>> pypgx.list_variants('CYP2D6', alleles=['*4'], assembly='GRCh37')
+    ['22-42524947-C-T']
+    >>> pypgx.list_variants('CYP2D6', alleles=['*4'], assembly='GRCh38')
+    ['22-42128945-C-T']
 
 However, there is one important caveat to consider if your sequencing data is
 GRCh38. That is, sequence reads must be aligned only to the main contigs
@@ -341,6 +390,29 @@
 genes.html>`__ page to see the list of genes with a genotype-phenotype
 table and each of their prediction method.
 
+To perform phenotype prediction with the API, you can use the
+``pypgx.predict_phenotype`` method:
+
+.. code:: python3
+
+    >>> import pypgx
+    >>> pypgx.predict_phenotype('CYP2D6', '*4', '*5')   # Both alleles have no function
+    'Poor Metabolizer'
+    >>> pypgx.predict_phenotype('CYP2D6', '*5', '*4')   # The order of alleles does not matter
+    'Poor Metabolizer'
+    >>> pypgx.predict_phenotype('CYP2D6', '*1', '*22')  # *22 has uncertain function
+    'Indeterminate'
+    >>> pypgx.predict_phenotype('CYP2D6', '*1', '*1x2') # Gene duplication
+    'Ultrarapid Metabolizer'
+
+To perform phenotype prediction with the CLI, you can use the
+``call-phenotypes`` command. It takes a ``SampleTable[Genotypes]`` file as
+input and outputs a ``SampleTable[Phenotypes]`` file:
+
+.. code-block:: text
+
+   $ pypgx call-phenotypes genotypes.zip phenotypes.zip
+
 Getting help
 ============
 
@@ -364,20 +436,23 @@
 Below is the list of submodules available in the API:
 
 {submodule_help}
-For getting help on a specific submodule (e.g. utils):
+For getting help on a specific submodule (e.g. ``utils``):
 
 .. code:: python3
 
    >>> from pypgx.api import utils
    >>> help(utils)
 
-For getting help on a specific method (e.g. predict_phenotype):
+For getting help on a specific method (e.g. ``predict_phenotype``):
 
 .. code:: python3
 
    >>> import pypgx
    >>> help(pypgx.predict_phenotype)
 
+In Jupyter Notebook and Lab, you can see the documentation for a python
+function by hitting ``SHIFT + TAB``. Hit it twice to expand the view.
+
 CLI examples
 ============