Answers.Body: <p>Setup:</p>
  
  <pre><code>import pandas as pd
  import pyranges as pr
  from pyranges import PyRanges
  from scipy.stats import fisher_exact
  import numpy as np
  
  # ! zcat dataset1.tsv.gz | head -2
  # chromosome    start   end num_motifs_in_group called_sites    called_sites_methylated methylated_frequency    group_sequence
  # chr21 5010053 5010053 1   3   0   0.000   CACCACGTCCA
  
  # avoid using more columns than we need, especially seq which contains object
  # is heavy to lug around
  colnames = ["Chromosome", "Start", "End", "calls", "methylated"]
  usecols = [0, 1, 2, 4, 5]
  
  h1 = PyRanges(pd.read_csv("dataset1.tsv.gz", sep="\t", names=colnames, header=0, usecols=usecols))
  h2 = PyRanges(pd.read_csv("dataset2.tsv.gz", sep="\t", names=colnames, header=0, usecols=usecols))
  
  regions = pr.read_bed("chr21.bed.gz")
  regions = regions[~regions.Chromosome.str.contains("_")].merge()
  
  regions.ID = list(range(len(regions)))
  
  regions
  # +--------------+-----------+-----------+-----------+
  # | Chromosome   | Start     | End       | ID        |
  # | (category)   | (int32)   | (int32)   | (int64)   |
  # |--------------+-----------+-----------+-----------|
  # | chr21        | 5065660   | 5065810   | 0         |
  # | chr21        | 5068140   | 5068290   | 1         |
  # | chr21        | 5128400   | 5128550   | 2         |
  # | chr21        | 5154780   | 5154930   | 3         |
  # | ...          | ...       | ...       | ...       |
  # | chr21        | 46637040  | 46637190  | 2377      |
  # | chr21        | 46638480  | 46638630  | 2378      |
  # | chr21        | 46639820  | 46639970  | 2379      |
  # | chr21        | 46692680  | 46692830  | 2380      |
  # +--------------+-----------+-----------+-----------+
  # Unstranded PyRanges object has 2,381 rows and 4 columns from 1 chromosomes.
  # For printing, the PyRanges was sorted on Chromosome.
  
  h1
  # +--------------+-----------+-----------+-----------+--------------+
  # | Chromosome   | Start     | End       | calls     | methylated   |
  # | (category)   | (int32)   | (int32)   | (int64)   | (int64)      |
  # |--------------+-----------+-----------+-----------+--------------|
  # | chr21        | 5010053   | 5010053   | 3         | 0            |
  # | chr21        | 5010215   | 5010215   | 5         | 1            |
  # | chr21        | 5010331   | 5010335   | 6         | 0            |
  # | chr21        | 5010525   | 5010525   | 2         | 2            |
  # | ...          | ...       | ...       | ...       | ...          |
  # | chr21        | 46699713  | 46699718  | 274       | 254          |
  # | chr21        | 46699739  | 46699739  | 133       | 113          |
  # | chr21        | 46699756  | 46699756  | 37        | 32           |
  # | chr21        | 46699770  | 46699773  | 156       | 130          |
  # +--------------+-----------+-----------+-----------+--------------+
  # Unstranded PyRanges object has 331,373 rows and 5 columns from 1 chromosomes.
  # For printing, the PyRanges was sorted on Chromosome.
  </code></pre>
  
  <p>Get total number of methylated per region:</p>
  
  <pre><code>def count_total_and_methylated(df):
      """For each region in bed, get sum of calls and methylated"""
  
      grpby = df.groupby("ID")
      total = grpby.calls.sum()
      methylated = grpby.methylated.sum()
  
      chromosome = df.Chromosome.iloc[0]
      start = grpby.Start.first()
      end = grpby.End.first()
      region_id = grpby.ID.first()
  
      return pd.DataFrame({"Chromosome": chromosome, "Start": start, "End": end, "ID": region_id,
                           "calls": total, "methylated": methylated})
  
  rh1 = regions.join(h1).drop(like="_b").apply(count_total_and_methylated)
  rh2 = regions.join(h2).drop(like="_b").apply(count_total_and_methylated)
  
  rh1.File = 1
  rh2.File = 2
  
  def add_missing(df, df2):
  
      """Since join is inner rh1 and rh2 might not have the exact same intervals. Fix this."""
  
      missing_in_first = np.setdiff1d(df2.ID, df.ID)
      first_missing = df2[df2.ID.isin(missing_in_first)].copy()
  
      f = df.File.head().iloc[0]
  
      first_missing.loc[:, "File"] = f
      first_missing.loc[:, "calls"] = 0
      first_missing.loc[:, "methylated"] = 0
  
      return pd.concat([df, first_missing])
  
  rh1 = rh1.apply_pair(rh2, add_missing)
  rh2 = rh2.apply_pair(rh1, add_missing)
  
  
  j = rh1.join(rh2, suffix=2).drop(like="File*|ID*|(Start|End)2")
  j
  # +--------------+-----------+-----------+-----------+--------------+-----------+---------------+
  # | Chromosome   | Start     | End       | calls     | methylated   | calls2    | methylated2   |
  # | (object)     | (int32)   | (int32)   | (int64)   | (int64)      | (int64)   | (int64)       |
  # |--------------+-----------+-----------+-----------+--------------+-----------+---------------|
  # | chr21        | 5065660   | 5065810   | 188       | 8            | 70        | 0             |
  # | chr21        | 5068140   | 5068290   | 168       | 166          | 74        | 72            |
  # | chr21        | 5128400   | 5128550   | 901       | 14           | 143       | 1             |
  # | chr21        | 5154780   | 5154930   | 653       | 0            | 593       | 0             |
  # | ...          | ...       | ...       | ...       | ...          | ...       | ...           |
  # | chr21        | 46637040  | 46637190  | 66        | 63           | 69        | 62            |
  # | chr21        | 46638480  | 46638630  | 28        | 0            | 30        | 3             |
  # | chr21        | 46639820  | 46639970  | 22        | 0            | 18        | 14            |
  # | chr21        | 46692680  | 46692830  | 71        | 25           | 57        | 23            |
  # +--------------+-----------+-----------+-----------+--------------+-----------+---------------+
  # Unstranded PyRanges object has 1,922 rows and 7 columns from 1 chromosomes.
  # For printing, the PyRanges was sorted on Chromosome.
  </code></pre>
  
  <p>Compute fisher exact using <a href="https://github.com/brentp/fishers_exact_test" rel="nofollow noreferrer">fisher</a> instead of scipy:</p>
  
  <pre><code>def fe(df):
      from fisher import pvalue_npy
  
      cols = "calls methylated calls2 methylated2".split()
      df[cols] = df[cols].astype(np.uint)
      _, _, two_sided = pvalue_npy(df.calls.values, df.methylated.values, df.calls2.values, df.methylated2.values)
  
      t1, t2, m1, m2 = df.calls, df.calls2, df.methylated, df.methylated2
  
      _or =   (t1 * m2) / (t2 * m1)
  
      df.insert(df.shape[1], "PValue", two_sided)
      df.insert(df.shape[1], "OR", _or)
  
      return df
  
  result = j.apply(fe, nb_cpu=1)
  result.drop(like="methylated*|calls*")
  # +--------------+-----------+-----------+------------------------+--------------------+
  # | Chromosome   | Start     | End       | PValue                 | OR                 |
  # | (object)     | (int32)   | (int32)   | (float64)              | (float64)          |
  # |--------------+-----------+-----------+------------------------+--------------------|
  # | chr21        | 5065660   | 5065810   | 0.1152954430992567     | 0.0                |
  # | chr21        | 5068140   | 5068290   | 0.9999999999995065     | 0.9846955389124064 |
  # | chr21        | 5128400   | 5128550   | 0.7075919513120225     | 0.4500499500499501 |
  # | chr21        | 5154780   | 5154930   | 1.0                    | nan                |
  # | ...          | ...       | ...       | ...                    | ...                |
  # | chr21        | 46638480  | 46638630  | 0.24262295081967644    | inf                |
  # | chr21        | 46639820  | 46639970  | 0.00024909594620266284 | inf                |
  # | chr21        | 46692680  | 46692830  | 0.7354565511951735     | 1.1459649122807019 |
  # | chr21        | 36784100  | 36784250  | 1.0                    | nan                |
  # +--------------+-----------+-----------+------------------------+--------------------+
  # Unstranded PyRanges object has 1,932 rows and 5 columns from 1 chromosomes.
  # For printing, the PyRanges was sorted on Chromosome.
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-11-12T14:09:57.567Z
Answers.Text: Thanks a lot for your reply. It seems like a decent solution but I need some more time to test it, as I'm focussing on getting a grant written these days...
Answers.UserDisplayName: Wouter De Coster
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-11-12T14:18:11.833Z
Answers.Text: Yeah, I'm writing a PR now. But with better code than here :)
Answers.UserDisplayName: The Unfun Cat
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-11-12T14:18:35.027Z
Answers.Text: I've tested that it gives the same results as your code. And it is super fast. Just wait for my PR plz
Answers.UserDisplayName: The Unfun Cat
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-11-01T10:01:40.383Z
Answers.LastActivityDate: 2019-11-01T14:13:46.460Z
Answers.LastEditDate: 2019-11-01T14:13:46.460Z
Answers.UserDisplayName: The Unfun Cat
Answers.Body: <p>Unless you need to use Python and can't use <code>subprocess</code>, here's a quick CLI one-liner which sums signal from sorted BED5 files, over the genomic space where they overlap:</p>
  
  <pre><code>$ bedmap --echo --sum --delim '\t' &lt;(bedops --merge A.bed B.bed ... N.bed) &lt;(bedops --everything A.bed B.bed ... N.bed) &gt; answer.bed
  </code></pre>
  
  <p>(Requires <code>bash</code> for process substitutions.)</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-10-31T15:48:44.357Z
Answers.LastActivityDate: 2019-10-31T15:48:44.357Z
Answers.UserDisplayName: Alex Reynolds
Body: <p>I am using <a href="https://github.com/biocore-ntnu/pyranges" rel="nofollow noreferrer">pyranges</a> and I have two datasets that I want to compare in several genome-wide intervals based on a bed file.</p>
  
  <p>For every interval in the bed file, I want to get all overlapping positions from the two datasets and sum two columns. That gives me 4 numbers for a fisher exact test, but that's not the most important part.</p>
  
  <p>My current <a href="https://github.com/wdecoster/methplotlib/blob/master/scripts/allele_specific_modification.py" rel="nofollow noreferrer">code</a> looks like below, and just iterates over the itertuples() of the bed file as pandas dataframe, and uses that to slice the pyranges object:</p>
  
  <pre><code>
  import pandas as pd
  import pyranges as pr
  from pyranges import PyRanges
  from scipy.stats import fisher_exact
  import numpy as np
  
  
  def main():
      colnames = ["Chromosome", "Start", "End", "num_motifs", "calls", "methylated", "freq", "seq"]
      h1 = PyRanges(pd.read_csv("dataset1.tsv.gz", sep="\t", names=colnames, header=0))
      h2 = PyRanges(pd.read_csv("dataset2.tsv.gz", sep="\t", names=colnames, header=0))
      for chr, begin, end in pr.read_bed("chr21.bed.gz").merge().as_df().itertuples(index=False, name=None):
          odr, pv = fisher_exact([meth_counts(h1[chr, begin:end]),
                                  meth_counts(h2[chr, begin:end])])
  
  
  def meth_counts(interval):
      if not interval.empty:
          all, methylated = interval.as_df()[['calls', "methylated"]].sum()
          return [all - methylated, methylated]
      else:
          return [0, 0]
  </code></pre>
  
  <p>This code works well, but I'm not sure if I use the "full power" of what pyranges could offer me for iterating and subsetting. I uploaded some example data on <a href="https://drive.google.com/drive/folders/1N5SIcmeFFGaQgrQSuinOyrmCaLw1f-cX?usp=sharing" rel="nofollow noreferrer">google drive</a> (250GB).</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-11-01T10:06:24.747Z
Comments.Text: Note that `gr[chr, begin:end]` will create a new NCLS to do lookups for every single iteration. This is sloooow. It is better to use functions like overlap, intersect, join etc.
Comments.UserDisplayName: The Unfun Cat
Tags: <python><bed><pyranges>
Title: Efficiently iterate over bed and subset pyranges objects
ThreadText: # Efficiently iterate over bed and subset pyranges objects
  
  (Asked by: username_0 on 2019-10-31)
  
  I am using [pyranges](https://github.com/biocore-ntnu/pyranges) and I have two datasets that I want to compare in several genome-wide intervals based on a bed file.
  
  For every interval in the bed file, I want to get all overlapping positions from the two datasets and sum two columns. That gives me 4 numbers for a fisher exact test, but that's not the most important part.
  
  My current [code](https://github.com/wdecoster/methplotlib/blob/master/scripts/allele_specific_modification.py) looks like below, and just iterates over the itertuples() of the bed file as pandas dataframe, and uses that to slice the pyranges object:
  
  ```
  
  import pandas as pd
  import pyranges as pr
  from pyranges import PyRanges
  from scipy.stats import fisher_exact
  import numpy as np
  
  def main():
      colnames = ["Chromosome", "Start", "End", "num_motifs", "calls", "methylated", "freq", "seq"]
      h1 = PyRanges(pd.read_csv("dataset1.tsv.gz", sep="\t", names=colnames, header=0))
      h2 = PyRanges(pd.read_csv("dataset2.tsv.gz", sep="\t", names=colnames, header=0))
      for chr, begin, end in pr.read_bed("chr21.bed.gz").merge().as_df().itertuples(index=False, name=None):
          odr, pv = fisher_exact([meth_counts(h1[chr, begin:end]),
                                  meth_counts(h2[chr, begin:end])])
  
  def meth_counts(interval):
      if not interval.empty:
          all, methylated = interval.as_df()[['calls', "methylated"]].sum()
          return [all - methylated, methylated]
      else:
          return [0, 0]
  
  ```
  
  This code works well, but I'm not sure if I use the "full power" of what pyranges could offer me for iterating and subsetting. I uploaded some example data on [google drive](https://drive.google.com/drive/folders/1N5SIcmeFFGaQgrQSuinOyrmCaLw1f-cX?usp=sharing) (250GB).
  
  **Tags:** python, bed, pyranges
  
  > **Comments:**
  > **username_1**: Note that `gr[chr, begin:end]` will create a new NCLS to do lookups for every single iteration. This is sloooow. It is better to use functions like overlap, intersect, join etc.  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_1 on 2019-11-01. Score: 2 | Accepted answer
  
  Setup:
  
  ```
  import pandas as pd
  import pyranges as pr
  from pyranges import PyRanges
  from scipy.stats import fisher_exact
  import numpy as np
  
  # ! zcat dataset1.tsv.gz | head -2
  # chromosome    start   end num_motifs_in_group called_sites    called_sites_methylated methylated_frequency    group_sequence
  # chr21 5010053 5010053 1   3   0   0.000   CACCACGTCCA
  
  # avoid using more columns than we need, especially seq which contains object
  # is heavy to lug around
  colnames = ["Chromosome", "Start", "End", "calls", "methylated"]
  usecols = [0, 1, 2, 4, 5]
  
  h1 = PyRanges(pd.read_csv("dataset1.tsv.gz", sep="\t", names=colnames, header=0, usecols=usecols))
  h2 = PyRanges(pd.read_csv("dataset2.tsv.gz", sep="\t", names=colnames, header=0, usecols=usecols))
  
  regions = pr.read_bed("chr21.bed.gz")
  regions = regions[~regions.Chromosome.str.contains("_")].merge()
  
  regions.ID = list(range(len(regions)))
  
  regions
  # +--------------+-----------+-----------+-----------+
  # | Chromosome   | Start     | End       | ID        |
  # | (category)   | (int32)   | (int32)   | (int64)   |
  # |--------------+-----------+-----------+-----------|
  # | chr21        | 5065660   | 5065810   | 0         |
  # | chr21        | 5068140   | 5068290   | 1         |
  # | chr21        | 5128400   | 5128550   | 2         |
  # | chr21        | 5154780   | 5154930   | 3         |
  # | ...          | ...       | ...       | ...       |
  # | chr21        | 46637040  | 46637190  | 2377      |
  # | chr21        | 46638480  | 46638630  | 2378      |
  # | chr21        | 46639820  | 46639970  | 2379      |
  # | chr21        | 46692680  | 46692830  | 2380      |
  # +--------------+-----------+-----------+-----------+
  # Unstranded PyRanges object has 2,381 rows and 4 columns from 1 chromosomes.
  # For printing, the PyRanges was sorted on Chromosome.
  
  h1
  # +--------------+-----------+-----------+-----------+--------------+
  # | Chromosome   | Start     | End       | calls     | methylated   |
  # | (category)   | (int32)   | (int32)   | (int64)   | (int64)      |
  # |--------------+-----------+-----------+-----------+--------------|
  # | chr21        | 5010053   | 5010053   | 3         | 0            |
  # | chr21        | 5010215   | 5010215   | 5         | 1            |
  # | chr21        | 5010331   | 5010335   | 6         | 0            |
  # | chr21        | 5010525   | 5010525   | 2         | 2            |
  # | ...          | ...       | ...       | ...       | ...          |
  # | chr21        | 46699713  | 46699718  | 274       | 254          |
  # | chr21        | 46699739  | 46699739  | 133       | 113          |
  # | chr21        | 46699756  | 46699756  | 37        | 32           |
  # | chr21        | 46699770  | 46699773  | 156       | 130          |
  # +--------------+-----------+-----------+-----------+--------------+
  # Unstranded PyRanges object has 331,373 rows and 5 columns from 1 chromosomes.
  # For printing, the PyRanges was sorted on Chromosome.
  
  ```
  
  Get total number of methylated per region:
  
  ```
  def count_total_and_methylated(df):
      """For each region in bed, get sum of calls and methylated"""
  
      grpby = df.groupby("ID")
      total = grpby.calls.sum()
      methylated = grpby.methylated.sum()
  
      chromosome = df.Chromosome.iloc[0]
      start = grpby.Start.first()
      end = grpby.End.first()
      region_id = grpby.ID.first()
  
      return pd.DataFrame({"Chromosome": chromosome, "Start": start, "End": end, "ID": region_id,
                           "calls": total, "methylated": methylated})
  
  rh1 = regions.join(h1).drop(like="_b").apply(count_total_and_methylated)
  rh2 = regions.join(h2).drop(like="_b").apply(count_total_and_methylated)
  
  rh1.File = 1
  rh2.File = 2
  
  def add_missing(df, df2):
  
      """Since join is inner rh1 and rh2 might not have the exact same intervals. Fix this."""
  
      missing_in_first = np.setdiff1d(df2.ID, df.ID)
      first_missing = df2[df2.ID.isin(missing_in_first)].copy()
  
      f = df.File.head().iloc[0]
  
      first_missing.loc[:, "File"] = f
      first_missing.loc[:, "calls"] = 0
      first_missing.loc[:, "methylated"] = 0
  
      return pd.concat([df, first_missing])
  
  rh1 = rh1.apply_pair(rh2, add_missing)
  rh2 = rh2.apply_pair(rh1, add_missing)
  
  j = rh1.join(rh2, suffix=2).drop(like="File*|ID*|(Start|End)2")
  j
  # +--------------+-----------+-----------+-----------+--------------+-----------+---------------+
  # | Chromosome   | Start     | End       | calls     | methylated   | calls2    | methylated2   |
  # | (object)     | (int32)   | (int32)   | (int64)   | (int64)      | (int64)   | (int64)       |
  # |--------------+-----------+-----------+-----------+--------------+-----------+---------------|
  # | chr21        | 5065660   | 5065810   | 188       | 8            | 70        | 0             |
  # | chr21        | 5068140   | 5068290   | 168       | 166          | 74        | 72            |
  # | chr21        | 5128400   | 5128550   | 901       | 14           | 143       | 1             |
  # | chr21        | 5154780   | 5154930   | 653       | 0            | 593       | 0             |
  # | ...          | ...       | ...       | ...       | ...          | ...       | ...           |
  # | chr21        | 46637040  | 46637190  | 66        | 63           | 69        | 62            |
  # | chr21        | 46638480  | 46638630  | 28        | 0            | 30        | 3             |
  # | chr21        | 46639820  | 46639970  | 22        | 0            | 18        | 14            |
  # | chr21        | 46692680  | 46692830  | 71        | 25           | 57        | 23            |
  # +--------------+-----------+-----------+-----------+--------------+-----------+---------------+
  # Unstranded PyRanges object has 1,922 rows and 7 columns from 1 chromosomes.
  # For printing, the PyRanges was sorted on Chromosome.
  
  ```
  
  Compute fisher exact using [fisher](https://github.com/brentp/fishers_exact_test) instead of scipy:
  
  ```
  def fe(df):
      from fisher import pvalue_npy
  
      cols = "calls methylated calls2 methylated2".split()
      df[cols] = df[cols].astype(np.uint)
      _, _, two_sided = pvalue_npy(df.calls.values, df.methylated.values, df.calls2.values, df.methylated2.values)
  
      t1, t2, m1, m2 = df.calls, df.calls2, df.methylated, df.methylated2
  
      _or =   (t1 * m2) / (t2 * m1)
  
      df.insert(df.shape[1], "PValue", two_sided)
      df.insert(df.shape[1], "OR", _or)
  
      return df
  
  result = j.apply(fe, nb_cpu=1)
  result.drop(like="methylated*|calls*")
  # +--------------+-----------+-----------+------------------------+--------------------+
  # | Chromosome   | Start     | End       | PValue                 | OR                 |
  # | (object)     | (int32)   | (int32)   | (float64)              | (float64)          |
  # |--------------+-----------+-----------+------------------------+--------------------|
  # | chr21        | 5065660   | 5065810   | 0.1152954430992567     | 0.0                |
  # | chr21        | 5068140   | 5068290   | 0.9999999999995065     | 0.9846955389124064 |
  # | chr21        | 5128400   | 5128550   | 0.7075919513120225     | 0.4500499500499501 |
  # | chr21        | 5154780   | 5154930   | 1.0                    | nan                |
  # | ...          | ...       | ...       | ...                    | ...                |
  # | chr21        | 46638480  | 46638630  | 0.24262295081967644    | inf                |
  # | chr21        | 46639820  | 46639970  | 0.00024909594620266284 | inf                |
  # | chr21        | 46692680  | 46692830  | 0.7354565511951735     | 1.1459649122807019 |
  # | chr21        | 36784100  | 36784250  | 1.0                    | nan                |
  # +--------------+-----------+-----------+------------------------+--------------------+
  # Unstranded PyRanges object has 1,932 rows and 5 columns from 1 chromosomes.
  # For printing, the PyRanges was sorted on Chromosome.
  
  ```
  
  > **Comments:**
  > **username_0**: Thanks a lot for your reply. It seems like a decent solution but I need some more time to test it, as I'm focussing on getting a grant written these days...  
  > **username_1**: Yeah, I'm writing a PR now. But with better code than here :)  
  > **username_1**: I've tested that it gives the same results as your code. And it is super fast. Just wait for my PR plz  
  
  ---
  
  ### Answer by username_2 on 2019-10-31. Score: 1
  
  Unless you need to use Python and can't use `subprocess`, here's a quick CLI one-liner which sums signal from sorted BED5 files, over the genomic space where they overlap:
  
  ```
  $ bedmap --echo --sum --delim '\t' <(bedops --merge A.bed B.bed ... N.bed) <(bedops --everything A.bed B.bed ... N.bed) > answer.bed
  
  ```
  
  (Requires `bash` for process substitutions.)
  
  
  ---
  
  