Merge branch 'master' of github.com:biocore-ntnu/pyranges

Author: endre bakken stovner

CHANGELOG.txt

@@ -1,3 +1,8 @@
+# 0.0.109 (16.09.21)
+- fix overlap invert behavior
+- add intersect invert flag
+- fix subtract in cases where more than 1024 intervals overlapped a single interval
+
 # 0.0.106/107/108(hotfixes) (07/8.09.21)
 - fix join with slack mutating first arg
 - add flag use_other_strand in join, nearest, k_nearest

pyranges/methods/intersection.py

@@ -61,7 +61,6 @@ def _intersection(scdf, ocdf, **kwargs):
 
 def _overlap(scdf, ocdf, **kwargs):
 
-    invert = kwargs["invert"]
     return_indexes = kwargs.get("return_indexes", False)
 
     if scdf.empty or ocdf.empty:
@@ -83,9 +82,6 @@ def _overlap(scdf, ocdf, **kwargs):
     else:
         _indexes = it.has_overlaps(starts, ends, indexes)
 
-    if invert:
-        _indexes = scdf.index.difference(_indexes)
-
     if return_indexes:
         return _indexes
 

pyranges/methods/subtraction.py

@@ -1,47 +1,87 @@
 import pandas as pd
-from ncls import NCLS
+import numpy as np
 
 
-def _subtraction(scdf, ocdf, **kwargs):
+def add_rows_per_group(df):
+    last_rows = df.groupby("__ix__").last().reset_index()
+    last_rows.loc[:, "__last__"] = True
+    df = pd.concat([df, last_rows], ignore_index=True)
+    df = df.sort_values("__ix__", ascending=True)
+    return df
 
-    if ocdf.empty or scdf.empty:
+def _subtraction(scdf, **kwargs):
+    if scdf.empty:
         return scdf
 
-    strandedness = kwargs["strandedness"]
-    strand = True if strandedness else False
+    falses = np.zeros(len(scdf), dtype=bool)
+    scdf.insert(scdf.shape[1], "__first__", falses)
+    scdf.insert(scdf.shape[1], "__last__", falses)
 
-    chromosome = scdf.Chromosome.head(1).iloc[0]
-    kwargs["chromosome"] = chromosome
+    scdf = add_rows_per_group(scdf)
 
-    if "Strand" in ocdf and strand:
-        strand = scdf.Strand.head(1).iloc[0]
-        kwargs["strand"] = strand
+    scdf.insert(scdf.shape[1], "NewStart", scdf.End__deleteme__.shift(fill_value=-1))
+    scdf.insert(scdf.shape[1], "NewEnd", scdf.Start__deleteme__)
+    scdf.insert(scdf.shape[1], "__ix2__", np.arange(len(scdf)))
 
-    o = NCLS(ocdf.Start.values, ocdf.End.values, ocdf.index.values)
+    first_rows = scdf.groupby(scdf.__ix__, as_index=False).first()
+    scdf.loc[scdf.__ix2__.isin(first_rows.__ix2__), "__first__"] = True
 
-    idx_self, new_starts, new_ends = o.set_difference_helper(
-        scdf.Start.values, scdf.End.values, scdf.index.values)
+    scdf.loc[:, "NewStart"] = np.where(scdf.__first__, scdf.Start, scdf.NewStart)
 
-    missing_idx = pd.Index(scdf.index).difference(idx_self)
+    scdf.loc[scdf.__first__ & ~(scdf.Start__deleteme__ >= scdf.Start), "NewStart"] = -1
+    scdf.loc[:, "NewEnd"] = np.where(scdf.__last__, scdf.End, scdf.NewEnd)
+    scdf.loc[:, "NewStart"] = np.where(scdf.__last__, scdf.End__deleteme__, scdf.NewStart)
+    scdf.loc[scdf.__last__ & ~(scdf.End__deleteme__ <= scdf.End), ["NewEnd", "NewStart"]] = -1
 
-    idx_to_drop = new_starts != -1
 
-    new_starts = new_starts[idx_to_drop]
-    new_ends = new_ends[idx_to_drop]
+    scdf = scdf[~((scdf.NewStart == -1) | (scdf.NewEnd == -1))]
+    scdf = scdf.drop(["Start", "End"], axis=1)
+    scdf.rename(columns={"NewStart": "Start", "NewEnd": "End"}, inplace=True)
 
-    idx_self = idx_self[idx_to_drop]
-    new_starts = pd.Series(new_starts, index=idx_self)
-    new_ends = pd.Series(new_ends, index=idx_self)
+    remove_mask = scdf.Start >= scdf.End
 
-    scdf = scdf.reindex(missing_idx.union(idx_self)).sort_index()
-    new_starts = new_starts.sort_index()
-    new_ends = new_ends.sort_index()
+    return scdf[~remove_mask]
 
-    if len(idx_self):
-        scdf.loc[scdf.index.isin(idx_self), "Start"] = new_starts.values
-        scdf.loc[scdf.index.isin(idx_self), "End"] = new_ends.values
+# def _subtraction(scdf, ocdf, **kwargs):
 
-    if not scdf.empty:
-        return scdf
-    else:
-        return None
+#     if ocdf.empty or scdf.empty:
+#         return scdf
+
+#     strandedness = kwargs["strandedness"]
+#     strand = True if strandedness else False
+
+#     chromosome = scdf.Chromosome.head(1).iloc[0]
+#     kwargs["chromosome"] = chromosome
+
+#     if "Strand" in ocdf and strand:
+#         strand = scdf.Strand.head(1).iloc[0]
+#         kwargs["strand"] = strand
+
+#     o = NCLS(ocdf.Start.values, ocdf.End.values, ocdf.index.values)
+
+#     idx_self, new_starts, new_ends = o.set_difference_helper(
+#         scdf.Start.values, scdf.End.values, scdf.index.values)
+
+#     missing_idx = pd.Index(scdf.index).difference(idx_self)
+
+#     idx_to_drop = new_starts != -1
+
+#     new_starts = new_starts[idx_to_drop]
+#     new_ends = new_ends[idx_to_drop]
+
+#     idx_self = idx_self[idx_to_drop]
+#     new_starts = pd.Series(new_starts, index=idx_self)
+#     new_ends = pd.Series(new_ends, index=idx_self)
+
+#     scdf = scdf.reindex(missing_idx.union(idx_self)).sort_index()
+#     new_starts = new_starts.sort_index()
+#     new_ends = new_ends.sort_index()
+
+#     if len(idx_self):
+#         scdf.loc[scdf.index.isin(idx_self), "Start"] = new_starts.values
+#         scdf.loc[scdf.index.isin(idx_self), "End"] = new_ends.values
+
+#     if not scdf.empty:
+#         return scdf
+#     else:
+#         return None

pyranges/pyranges.py

@@ -1854,7 +1854,7 @@ def insert(self, other, loc=None):
         return self
 
 
-    def intersect(self, other, strandedness=None, how=None, nb_cpu=1):
+    def intersect(self, other, strandedness=None, how=None, invert=False, nb_cpu=1):
 
         """Return overlapping subintervals.
 
@@ -1877,6 +1877,10 @@ def intersect(self, other, strandedness=None, how=None, nb_cpu=1):
             What intervals to report. By default reports all overlapping intervals. "containment"
             reports intervals where the overlapping is contained within it.
 
+        invert : bool, default False
+
+            Whether to return the intervals without overlaps.
+
         nb_cpu: int, default 1
 
             How many cpus to use. Can at most use 1 per chromosome or chromosome/strand tuple.
@@ -1962,9 +1966,21 @@ def intersect(self, other, strandedness=None, how=None, nb_cpu=1):
         kwargs = fill_kwargs(kwargs)
         kwargs["sparse"] = {"self": False, "other": True}
 
+        if len(self) == 0:
+            return self
+
+        if invert:
+            self.__ix__ = np.arange(len(self))
+
         dfs = pyrange_apply(_intersection, self, other, **kwargs)
+        result = pr.PyRanges(dfs)
 
-        return PyRanges(dfs)
+        if invert:
+            found_idxs = getattr(result, "__ix__", [])
+            result = self[~self.__ix__.isin(found_idxs)]
+            result = result.drop("__ix__")
+
+        return result
 
     def items(self):
 
@@ -3195,6 +3211,10 @@ def overlap(self, other, strandedness=None, how="first", invert=False, nb_cpu=1)
             What intervals to report. By default reports every interval in self with overlap once.
             "containment" reports all intervals where the overlapping is contained within it.
 
+        invert : bool, default False
+
+            Whether to return the intervals without overlaps.
+
         nb_cpu: int, default 1
 
             How many cpus to use. Can at most use 1 per chromosome or chromosome/strand tuple.
@@ -3292,9 +3312,22 @@ def overlap(self, other, strandedness=None, how="first", invert=False, nb_cpu=1)
         kwargs["invert"] = invert
         kwargs = fill_kwargs(kwargs)
 
+        if len(self) == 0:
+            return self
+
+        if invert:
+            self = self.copy()
+            self.__ix__ = np.arange(len(self))
+
         dfs = pyrange_apply(_overlap, self, other, **kwargs)
+        result = pr.PyRanges(dfs)
 
-        return pr.PyRanges(dfs)
+        if invert:
+            found_idxs = getattr(result, "__ix__", [])
+            result = self[~self.__ix__.isin(found_idxs)]
+            result = result.drop("__ix__")
+
+        return result
 
     def pc(self, n=8, formatting=None):
 
@@ -4343,16 +4376,27 @@ def subtract(self, other, strandedness=None, nb_cpu=1):
 
         from pyranges.methods.subtraction import _subtraction
 
-        kwargs = {"strandedness": strandedness}
+        strand = True if strandedness in ["same", "opposite"] else False
+        kwargs = {"strandedness": strandedness, "strand": strand}
         kwargs["sparse"] = {"self": False, "other": True}
         kwargs = fill_kwargs(kwargs)
 
         strand = True if strandedness else False
-        other_clusters = other.merge(strand=strand)
-
-        result = pyrange_apply(_subtraction, self, other_clusters, **kwargs)
+        if len(other) != 0:
+            cols_original_order = self.columns.copy()
+            other_clusters = other.merge(strand=strand)
+            self.__ix__ = np.arange(len(self))
+            no_overlap = self.overlap(other_clusters, strandedness=strandedness, invert=True)
+            j = self.sort().join(other_clusters, strandedness=strandedness, suffix="__deleteme__")
+            result = pyrange_apply_single(_subtraction, j, **kwargs)
+            result = {k: v for (k, v) in result.items()}
+            result = PyRanges(result)
+            result = pr.concat([result, no_overlap])
+            result = PyRanges({k: v[cols_original_order] for k, v in result.items()})
+        else:
+            result = self
 
-        return PyRanges(result)
+        return result
 
     def summary(self, to_stdout=True, return_df=False):
 

pyranges/version.py

@@ -1,2 +1 @@
-__version__ = "0.0.108"
-
+__version__ = "0.0.109"

tests/helpers.py

@@ -3,6 +3,12 @@
 
 def assert_df_equal(df1, df2):
 
+    print("-"*100)
+    print("df1")
+    print(df1)
+    print("df2")
+    print(df2)
+
     # df1.loc[:, "Start"] = df1.Start.astype(np.int64)
     # df2.loc[:, "Start"] = df1.Start.astype(np.int64)
     # df1.loc[:, "End"] = df1.End.astype(np.int64)

tests/test_binary.py

@@ -304,9 +304,12 @@ def test_subtraction(gr, gr2, strandedness):
         names="Chromosome Start End Name Score Strand".split(),
         sep="\t")
 
+    print("subtracting" * 50)
     result = gr.subtract(gr2, strandedness=strandedness)
 
+    print("bedtools_result")
     print(bedtools_df)
+    print("PyRanges result:")
     print(result)
 
     compare_results(bedtools_df, result)

tests/test_do_not_error.py

@@ -50,6 +50,8 @@
 # @reproduce_failure('5.5.4', b'AXicY2RAA4xIJDY+AAC2AAY=') # test_three_in_a_row[strandedness_chain24-method_chain24]
 def test_three_in_a_row(gr, gr2, gr3, strandedness_chain, method_chain):
 
+    print(method_chain)
+
     s1, s2 = strandedness_chain
     f1, f2 = method_chain