0
$\begingroup$

I have two bed files where each column has a numeric score. I want to sum these numbers where there is overlap between the two bed files. How do I do that? Also, I am only interested in the score in specific genes.

Here is my setup:

import numpy as np
np.random.seed(0)
import pyranges as pr
import pandas as pd

entries = int(1e6)
cpg1 = pr.random(entries)
cpg2 = pr.random(entries)

cpg1.CPG = np.random.randint(100, size=entries)
cpg2.CPG = np.random.randint(100, size=entries)

regulatory_elements = pr.random(int(1e5), length=int(1e4))
regulatory_elements.Gene = np.arange(len(regulatory_elements), dtype=int)

So my data looks like:

cpg1
# +--------------+-----------+-----------+--------------+-----------+
# | Chromosome   | Start     | End       | Strand       | CPG       |
# | (category)   | (int32)   | (int32)   | (category)   | (int64)   |
# |--------------+-----------+-----------+--------------+-----------|
# | chr1         | 8830650   | 8830750   | +            | 72        |
# | chr1         | 9564361   | 9564461   | +            | 57        |
# | chr1         | 44977425  | 44977525  | +            | 99        |
# | chr1         | 239741543 | 239741643 | +            | 54        |
# | ...          | ...       | ...       | ...          | ...       |
# | chrY         | 29437476  | 29437576  | -            | 44        |
# | chrY         | 49995298  | 49995398  | -            | 43        |
# | chrY         | 50840129  | 50840229  | -            | 10        |
# | chrY         | 38069647  | 38069747  | -            | 89        |
# +--------------+-----------+-----------+--------------+-----------+
# Stranded PyRanges object has 1,000,000 rows and 5 columns from 25 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.


regulatory_elements
# +--------------+-----------+-----------+--------------+-----------+
# | Chromosome   | Start     | End       | Strand       | Gene      |
# | (category)   | (int32)   | (int32)   | (category)   | (int64)   |
# |--------------+-----------+-----------+--------------+-----------|
# | chr1         | 217921634 | 217931634 | +            | 0         |
# | chr1         | 166226804 | 166236804 | +            | 1         |
# | chr1         | 170688210 | 170698210 | +            | 2         |
# | chr1         | 57958563  | 57968563  | +            | 3         |
# | ...          | ...       | ...       | ...          | ...       |
# | chrY         | 35870043  | 35880043  | -            | 99996     |
# | chrY         | 35057634  | 35067634  | -            | 99997     |
# | chrY         | 52718975  | 52728975  | -            | 99998     |
# | chrY         | 19196507  | 19206507  | -            | 99999     |
# +--------------+-----------+-----------+--------------+-----------+
# Stranded PyRanges object has 100,000 rows and 5 columns from 25 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.

(This is a q from privateish correspondence).

$\endgroup$
1
  • $\begingroup$ Mods: would love tags for pandas, pyranges and run-length-encodings. All three have tags on SO you can copy :) $\endgroup$ Oct 29 '19 at 10:58
0
$\begingroup$

If you want to sum two columns based on overlaps, converting them to RLEs is a good idea.

# get the cpg scores only where there is overlap with the genes you are interested in 
cr1 = cpg1.join(regulatory_elements).new_position("intersection").drop(like="_(a|b)")
cr2 = cpg2.join(regulatory_elements).new_position("intersection").drop(like="_(a|b)")

# convert these to run length encodings, using the score from the CPG column
rle1 = cr1.to_rle("CPG", rpm=False)
rle2 = cr2.to_rle("CPG", rpm=False)

rle_sum = rle1 + rle2

# convert the rles back to ranges
# (I will give the to_ranges function an option dtype eventually)
range_sum = rle_sum.to_ranges().apply(lambda df: df.astype({"Start": np.int32, "End": np.int32}))

# add the gene names to the scores
gr = range_sum.join(regulatory_elements) # .drop(like="_b")

Now the result is:

print(gr)
# +--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------+
# | Chromosome   | Start     | End       | Score       | Strand     | Start_b   | End_b     | Strand_b     | Gene      |
# | (object)     | (int32)   | (int32)   | (float64)   | (object)   | (int32)   | (int32)   | (category)   | (int64)   |
# |--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------|
# | chr1         | 15623     | 15723     | 77.0        | +          | 9886      | 19886     | +            | 3714      |
# | chr1         | 16034     | 16134     | 84.0        | +          | 9886      | 19886     | +            | 3714      |
# | chr1         | 36588     | 36688     | 2.0         | +          | 36479     | 46479     | +            | 2871      |
# | chr1         | 37271     | 37371     | 30.0        | +          | 36479     | 46479     | +            | 2871      |
# | ...          | ...       | ...       | ...         | ...        | ...       | ...       | ...          | ...       |
# | chrY         | 59312851  | 59312951  | 28.0        | -          | 59308107  | 59318107  | +            | 98599     |
# | chrY         | 59332843  | 59332943  | 82.0        | -          | 59332664  | 59342664  | -            | 99215     |
# | chrY         | 59337592  | 59337692  | 41.0        | -          | 59332664  | 59342664  | -            | 99215     |
# | chrY         | 59341719  | 59341819  | 15.0        | -          | 59332664  | 59342664  | -            | 99215     |
# +--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------+
# Stranded PyRanges object has 668,317 rows and 9 columns from 25 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.

print(gr[gr.Score > 100])
# +--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------+
# | Chromosome   | Start     | End       | Score       | Strand     | Start_b   | End_b     | Strand_b     | Gene      |
# | (object)     | (int32)   | (int32)   | (float64)   | (object)   | (int32)   | (int32)   | (category)   | (int64)   |
# |--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------|
# | chr1         | 96874     | 96974     | 116.0       | +          | 87668     | 97668     | -            | 7391      |
# | chr1         | 96874     | 96974     | 116.0       | +          | 89794     | 99794     | -            | 6357      |
# | chr1         | 749686    | 749697    | 180.0       | +          | 744041    | 754041    | +            | 1150      |
# | chr1         | 749686    | 749697    | 180.0       | +          | 749277    | 759277    | -            | 6395      |
# | ...          | ...       | ...       | ...         | ...        | ...       | ...       | ...          | ...       |
# | chrY         | 59060667  | 59060767  | 196.0       | -          | 59052942  | 59062942  | -            | 99957     |
# | chrY         | 59060667  | 59060767  | 196.0       | -          | 59059904  | 59069904  | +            | 98401     |
# | chrY         | 59063570  | 59063576  | 105.0       | -          | 59059904  | 59069904  | +            | 98401     |
# | chrY         | 59269010  | 59269085  | 152.0       | -          | 59264704  | 59274704  | -            | 99499     |
# +--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------+
# Stranded PyRanges object has 106,643 rows and 9 columns from 25 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.