Efficiently iterate over bed and subset pyranges objects

Question

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).

Answer 1

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.

Answer 2

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.)