text: user681: I am using 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 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 (250GB).
  
  user776: 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.)
  
  user2688: 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 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.
  
  
original_text: 681: <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>
   776: <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>
   2688: <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>
  