Merge pull request #26 from jrm5100/master

Version 0.11.0 : Add Region

Author: jrm5100

docs/release-history.rst

@@ -2,6 +2,11 @@
 Release History
 ===============
 
+v0.11.0 (2021-07-21)
+--------------------
+
+Add a *Region* scalar with support for loading from BED files and filtering
+
 v0.10.1 (2021-07-16)
 --------------------
 

pandas_genomics/accessors/dataframe_accessor.py

@@ -1,11 +1,12 @@
 from collections import Counter
-from typing import Optional, List
+from typing import Optional, List, Union
 
 import numpy as np
 import pandas as pd
 
 from .utils import calculate_edge_alphas
 from pandas_genomics.arrays import GenotypeDtype
+from ..scalars import Region
 
 
 @pd.api.extensions.register_dataframe_accessor("genomics")
@@ -300,3 +301,53 @@ def filter_variants_hwe(self, cutoff: float = 0.05) -> pd.DataFrame:
             :, (genotype_hwe_pval < cutoff) & ~np.isnan(genotype_hwe_pval)
         ].columns
         return self._obj.drop(columns=removed)
+
+    def in_regions(self, regions: Union[Region, List[Region]]):
+        """
+        Keep variants that exist in one or more of the specified regions
+
+        Parameters
+        ----------
+        regions: Region or List[Region]
+
+        Returns
+        -------
+        pd.DataFrame
+        """
+        if isinstance(regions, Region):
+            regions = [
+                regions,
+            ]
+        genotypes = self._obj.select_dtypes([GenotypeDtype])
+        boolean_array = genotypes.apply(lambda s: False)
+        for r in regions:
+            boolean_array = boolean_array | genotypes.apply(
+                lambda s: s.genomics.contained_by(r)
+            )
+        removed = boolean_array[~boolean_array].index  # Remove where False
+        return self._obj.drop(columns=removed)
+
+    def not_in_regions(self, regions: Union[Region, List[Region]]):
+        """
+        Remove variants that exist in one or more of the specified regions
+
+        Parameters
+        ----------
+        regions: Region or List[Region]
+
+        Returns
+        -------
+        pd.DataFrame
+        """
+        if isinstance(regions, Region):
+            regions = [
+                regions,
+            ]
+        genotypes = self._obj.select_dtypes([GenotypeDtype])
+        boolean_array = genotypes.apply(lambda s: False)
+        for r in regions:
+            boolean_array = boolean_array | genotypes.apply(
+                lambda s: s.genomics.contained_by(r)
+            )
+        removed = boolean_array[boolean_array].index  # Remove where True
+        return self._obj.drop(columns=removed)

pandas_genomics/accessors/series_accessor.py

@@ -1,9 +1,10 @@
-from typing import Optional, List
+from typing import Optional, List, Union
 
 import pandas as pd
 
 from .utils import calculate_edge_alphas
 from pandas_genomics.arrays import GenotypeDtype
+from ..scalars import Region
 
 
 @pd.api.extensions.register_series_accessor("genomics")
@@ -72,9 +73,6 @@ def maf(self):
         See :py:attr:`GenotypeArray.maf`"""
         return self._array.maf
 
-    #########################
-    # Calculated Properties #
-    #########################
     @property
     def hwe_pval(self):
         """Return the probability that the samples are in HWE
@@ -226,3 +224,29 @@ def calculate_edge_encoding_values(
     ##############
     # QC Methods #
     ##############
+
+    # TODO
+
+    #################
+    # Other Methods #
+    #################
+
+    def contained_by(self, regions: Union[Region, List[Region]]):
+        """
+        True if the variant is contained within the specified region(s)
+
+        Parameters
+        ----------
+        regions: Region or List[Region]
+
+        Returns
+        -------
+        bool
+        """
+        if isinstance(regions, Region):
+            return regions.contains_variant(self.variant)
+        else:
+            for r in regions:
+                if r.contains_variant(self.variant):
+                    return True
+            return False

pandas_genomics/io/__init__.py

@@ -1,17 +1,20 @@
 """
 Input/Output
 ------------
-The `io` module provides functions for loading and saving GenotypeArrays to common variant formats
+The `io` module provides functions for loading and saving GenotypeArrays or scalar values to common variant formats
 
 .. autosummary::
      :toctree: io
 
      from_plink
      to_plink
      from_vcf
+
+     from_bed
 """
 
+from .bed import from_bed
 from .plink import from_plink, to_plink
 from .vcf import from_vcf
 
-__all__ = ["from_plink", "to_plink", "from_vcf"]
+__all__ = ["from_bed", "from_plink", "to_plink", "from_vcf"]

pandas_genomics/io/bed.py

@@ -0,0 +1,36 @@
+from pathlib import Path
+from typing import Union
+
+import pandas as pd
+
+from pandas_genomics.scalars import Region
+
+
+def from_bed(filename: Union[str, Path]):
+    """
+    Yields genomic regions from a bed file as Region scalars
+
+    Parameters
+    ----------
+    filename: str or Path
+        bed file
+
+    Returns
+    -------
+    List[Region]
+    """
+    with open(filename, "r") as input:
+        for idx, line in enumerate(input):
+            if (
+                line.startswith("browser")
+                or line.startswith("track")
+                or line.startswith("#")
+            ):
+                continue
+            line = line.split("\t")
+            if len(line) < 3:
+                raise ValueError(
+                    f"Expected at least 3 tab-delimited fields, found {len(line)} in '{line}'."
+                )
+            # Note that positions need to be incremented by 1 to go from 0-based to 1-based.
+            yield Region(line[0], int(line[1]) + 1, int(line[2]) + 1, line[3])

pandas_genomics/scalars.py

@@ -1,18 +1,22 @@
 """
 Scalars
 -------
-This module contains scalar types used in the ExtensionArrays.  They may also be useful on their own.
+This module contains scalar types, some of which are used in the ExtensionArrays.  They may also be useful on their own.
 
 .. autosummary::
      :toctree: scalars
 
      Variant
      Genotype
+     Region
 """
 import uuid
+from dataclasses import dataclass
 from typing import Optional, List, Tuple, Union
 
-MISSING_IDX = 255  # Integer indicating a missing allele or genotype score.  Each variant must have 254 alleles max and the maximum genotype score is 254.
+# Integer indicating a missing allele or genotype score.
+# Each variant must have 254 alleles max and the maximum genotype score is 254.
+MISSING_IDX = 255
 
 
 class Variant:
@@ -509,3 +513,51 @@ def _float_score(self):
             return 255
         else:
             return self.score
+
+
+@dataclass(order=True)
+class Region:
+    chromosome: str
+    start: int
+    end: int
+    name: str = ""
+    """
+    Information associated with a genomic region.
+
+    Parameters
+    ----------
+    chomosome: str
+    start, end: int
+        1-based, open-ended (includes start, not end)
+    name: str
+        name/comment for the region
+
+    Examples
+    --------
+    >>> region1 = Region('chr12', 12345, 12387)
+    >>> region2 = Region('chr12', 5236373, 5246380)
+    >>> print(region1 < region2)
+    True
+    """
+
+    def __post_init__(self):
+        """Run after init, use to validate input"""
+        if not isinstance(self.chromosome, str):
+            raise TypeError(f"chromosome should be a str, not {type(self.chromosome)}")
+        if not isinstance(self.start, int) or not isinstance(self.end, int):
+            raise TypeError("start and end must be integers")
+        if self.start < 1 or self.end < 1:
+            raise ValueError("start and end must both be > 0 (1-based indexing)")
+        if self.start >= self.end:
+            raise ValueError(
+                "start must be <= end, since the start is included and the end is not."
+            )
+
+    def contains_variant(self, var: Variant):
+        if var.chromosome != self.chromosome:
+            return False
+        if var.position < self.start:
+            return False
+        if var.position >= self.end:
+            return False
+        return True

pyproject.toml

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "pandas-genomics"
-version = "0.10.1"
+version = "0.11.0"
 description = "Pandas ExtensionDtypes and ExtensionArray for working with genomics data"
 license = "BSD-3-Clause"
 authors = ["John McGuigan <jrm5100@psu.edu>"]

tests/data/bed/bed_test.bed

@@ -0,0 +1,6 @@
+browser position chr7:127471196-127495720
+browser hide all
+track name="ItemRGBDemo" description="Item RGB demonstration" visibility=2 itemRgb="On"
+chr7	127471196	127472363	Pos1	0	+	127471196	127472363	255,0,0
+chr7	127472363	127473530	Pos2	0	+	127472363	127473530	255,0,0
+chr7	127473530	127474697	Pos3	0	+	127473530	127474697	255,0,0

tests/genotype_array/test_GenotypeArrayAccessors.py

@@ -6,8 +6,12 @@
 import pytest
 from pandas._testing import (
     assert_series_equal,
+    assert_frame_equal,
 )
 
+from pandas_genomics import GenotypeArray
+from pandas_genomics.scalars import Variant, Region
+
 
 def test_variant_score(data, data_for_encoding):
     assert pd.Series(data).genomics.variant.score == float(data.variant.score)
@@ -59,3 +63,52 @@ def test_filter_hwe(ga_inhwe, ga_nothwe, filter_value, num_cols_left):
     else:
         result = data.genomics.filter_variants_hwe(filter_value)
     assert len(result.columns) == num_cols_left + 1
+
+
+def test_region_series():
+    var = Variant("chr1", position=999, ref="A", alt=["a"])
+    s = pd.Series(
+        GenotypeArray(
+            [
+                var.make_genotype_from_str("A/A"),
+            ]
+            * 10
+        )
+    )
+    assert s.genomics.contained_by(Region("chr1", 900, 1000))
+    assert not s.genomics.contained_by(Region("chr2", 900, 1000))
+    assert not s.genomics.contained_by(Region("chr1", 900, 999))
+
+
+def test_region_df():
+    var1 = Variant("chr1", position=999, ref="A", alt=["a"])
+    var2 = Variant("chr1", position=6789, ref="A", alt=["a"])
+    var3 = Variant("chr2", position=999, ref="A", alt=["a"])
+    var4 = Variant("chr3", position=25622, ref="A", alt=["a"])
+    df = pd.DataFrame(
+        {
+            f"var{idx+1}": GenotypeArray(
+                [
+                    var.make_genotype_from_str("A/A"),
+                ]
+                * 10
+            )
+            for idx, var in enumerate([var1, var2, var3, var4])
+        }
+    )
+    assert_frame_equal(
+        df.genomics.in_regions(Region("chr1", 900, 1000)),
+        df[
+            [
+                "var1",
+            ]
+        ],
+    )
+    assert_frame_equal(
+        df.genomics.not_in_regions(Region("chr1", 900, 1000)),
+        df[["var2", "var3", "var4"]],
+    )
+    assert_frame_equal(
+        df.genomics.in_regions([Region("chr1", 900, 1000), Region("chr2", 900, 1000)]),
+        df[["var1", "var3"]],
+    )

tests/io/test_bed.py

@@ -0,0 +1,16 @@
+from pathlib import Path
+
+from pandas_genomics import io
+from pandas_genomics.scalars import Region
+
+DATA_DIR = Path(__file__).parent.parent / "data" / "bed"
+
+
+def test_bed():
+    input = DATA_DIR / "bed_test.bed"
+    result = io.from_bed(input)
+    assert list(result) == [
+        Region("chr7", 127471197, 127472364, "Pos1"),
+        Region("chr7", 127472364, 127473531, "Pos2"),
+        Region("chr7", 127473531, 127474698, "Pos3"),
+    ]