<p>Quick and dirty (not to mention atrociously bad <code>O(n)</code>) solution:</p>
<pre><code>import pyranges as pr
import numpy as np
np.random.seed(42 * 10)

# create large df to test on
gr = pr.random(int(1e5), length=10000, chromsizes={&quot;chr1&quot;: 249250621})
gr.Score = np.random.randint(250, size=len(gr))

def remove_worst_scores_until_no_overlap(gr):
    df = gr.df
    old_length = -1
    new_length = len(df)
    while old_length != new_length:
        df = df.drop([&quot;Cluster&quot;, &quot;Count&quot;], axis=1, errors=&quot;ignore&quot;)
        df = pr.PyRanges(df).cluster(count=True).df.sort_values(&quot;Score&quot;, ascending=False)
        df = df[df.duplicated('Cluster', keep='last') | (df.Count == 1)]
        old_length = new_length
        new_length = len(df)

    return pr.PyRanges(df.drop(&quot;Cluster Count&quot;.split(), axis=1)).sort()

hotspots_no_overlaps = remove_worst_scores_until_no_overlap(gr)
</code></pre>

CC BY-SA 4.0
2021-03-02T15:01:31.793Z
Quick and dirty, but it works!
Whitehot
CC BY-SA 4.0
2021-03-01T19:40:42.117Z
2021-03-01T19:40:42.117Z
The Unfun Cat
<p>Answer to main question:</p>
<pre><code>import pyranges as pr
assert int(pr.__version__.split(&quot;.&quot;)[2]) &gt;= 93, &quot;pip install pyranges==0.0.93&quot;
import numpy as np
np.random.seed(42 * 10)

# create large df to test on
gr = pr.random(int(1e5), length=1, chromsizes={&quot;chr1&quot;: 249250621})
gr.Score = np.random.randint(250, size=len(gr))

result = gr.max_disjoint(slack=5000)

# small example that is easy to validate
gr2 = pr.from_string(&quot;&quot;&quot;Chromosome Start End Score
chr1 1 2 5
chr1 4999 5000 10
chr1 10001 10002 1
chr1 15000 15001 2
chr1 20100 20101 5
&quot;&quot;&quot;)

gr2.max_disjoint(slack=5000)

# +--------------+-----------+-----------+-----------+
# | Chromosome   |     Start |       End |     Score |
# | (category)   |   (int32) |   (int32) |   (int64) |
# |--------------+-----------+-----------+-----------|
# | chr1         |         1 |         2 |         5 |
# | chr1         |     10001 |     10002 |         1 |
# | chr1         |     20100 |     20101 |         5 |
# +--------------+-----------+-----------+-----------+
# Unstranded PyRanges object has 3 rows and 4 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.
</code></pre>
<hr />
<p>Could you please post an example of what data you have and what you want to end up with?</p>
<ol>
<li><p>Inserting row after row is going to be slow, whether it is pandas or pyranges. It is best to collect all the data and then create a new dataframe or PyRanges.</p>
</li>
<li><p><code>pr.itergrs</code> is for iterating over multiple pyranges at the same time. It takes a list of pyranges while you seem to have a collection of tuples. The function iterates over every chromosome/strand pair in the list and for it extracts the data from every pyranges. See the code example at the end.</p>
</li>
</ol>
<p>To iterate over a pyranges you can do <code>for k, df in gr: # ...</code>.</p>
<p>Anyways, to do what you want:</p>
<pre><code>import pyranges as pr
import numpy as np

f1 = pr.data.f1()
f1.Start = ((f1.Start + f1.End)/2).astype(np.int32) # int is 64 bit
f1.End = f1.Start + 1
</code></pre>
<ol start="3">
<li>I understand I have not done the addititonal processing you speak of, but if you can post an example I'm sure I'd be able to help you :)</li>
</ol>
<p>Code example <code>pr.itergrs</code>:</p>
<pre><code>import pyranges as pr
f1 = pr.data.f1()
f2 = pr.data.f2()
for k, dfs in pr.itergrs([f1, f2], keys=True):
    print(k)
     for df in dfs:
         print(&quot;df:&quot;)
         print(df)

('chr1', '+')
df:
  Chromosome  Start  End       Name  Score Strand
0       chr1      3    6  interval1      0      +
2       chr1      8    9  interval3      0      +
df:
  Chromosome  Start  End Name  Score Strand
0       chr1      1    2    a      0      +
('chr1', '-')
df:
  Chromosome  Start  End       Name  Score Strand
1       chr1      5    7  interval2      0      -
df:
  Chromosome  Start  End Name  Score Strand
1       chr1      6    7    b      0      -
</code></pre>

CC BY-SA 4.0
2021-02-24T18:03:43.613Z
Thanks for the pointers, this already helps a lot. I'll post my code in the morning (it's 7PM for me now)
Whitehot
CC BY-SA 4.0
2021-02-25T14:37:08.390Z
Added the information you asked for :)
Whitehot
CC BY-SA 4.0
2021-02-25T15:37:58.293Z
Thanks for the solution, do you think you could comment the function a bit more so I can be sure of what each part of it does? From what I can tell, you sort and cluster the df by score in f, then merge all ranges withing 5kb of each other, is that correct?
Whitehot
CC BY-SA 4.0
2021-02-25T15:48:23.380Z
I collect all ranges within 5k and put them in separate clusters. Then the function f in apply merely picks out the row within each cluster with the highest score and discards the rest. I might have made a mistake in that if you have hotspots at 0, 4000, 8000 and 12000 then you would probably want hotspot 0 and 8000 while I discard all but one. Will think about the problem some more and update the answer tomorrow.
The Unfun Cat
CC BY-SA 4.0
2021-02-25T15:55:01.813Z
Yes I was wondering if it behaved like that. That's why I was thinking of creating "flanking" regions to better define which hotspots to discard, but if you find a more elegant solution, I'm all for it!
Whitehot
CC BY-SA 4.0
2021-02-25T20:53:35.233Z
I've now added a method `max_disjoint` to PyRanges in 0.0.93. Thanks for making me realise there was a fundamental interval operation I was missing.
The Unfun Cat
CC BY-SA 4.0
2021-02-25T22:21:10.023Z
Surely "Dev adds new function to package because of you" should be a medal on bioinformatics SE, right? Thank you for your help in any case! I may be calling upon your services later depending on how the project advances
Whitehot
CC BY-SA 4.0
2021-02-26T05:40:30.537Z
I do this for fun, but thanks. Feel free to accept the answer if everything worked as expected.
The Unfun Cat
CC BY-SA 4.0
2021-02-26T11:53:45.050Z
So the `max_disjoint()` works great but when I combine it with the centering you gave, it throws an error: `ValueError: Buffer dtype mismatch, expected 'const int32_t' but got 'long long'`. Do you want the full traceback?
Whitehot
CC BY-SA 4.0
2021-02-26T19:54:25.170Z
Will fix. Could you post an issue on the pyranges github?
The Unfun Cat
CC BY-SA 4.0
2021-03-01T15:41:22.973Z
Thank you for all your help, but `max_disjoint()` has a behaviour that could be improved on I think: for now it seems to start at the first range on the chromosome, whereas it would be useful to be able to specify a column for it to sort on, like the score, before searching for overlaps.
Whitehot
CC BY-SA 4.0
2021-03-01T16:01:00.537Z
Actually I'm not entirely sure what it uses to sort the pyranges on, seeing as even in your example, it's kept the `10001;10002;1` range as the centre of the second cluster
Whitehot
CC BY-SA 4.0
2021-03-01T18:02:48.830Z
max_disjoint() is actually a well-known greedy algorithm: https://en.wikipedia.org/wiki/Maximum_disjoint_set#1-dimensional_intervals:_exact_polynomial_algorithm It starts by sorting on the end column. I will continue thinking about the problem
The Unfun Cat
CC BY-SA 4.0
2021-02-24T17:41:08.803Z
2021-02-25T20:52:27.593Z
2021-02-25T20:52:27.593Z
The Unfun Cat
<p>I have a pyranges object that I am looking to iterate over. For each range, I would like to &quot;resize&quot; it to its centre (so <code>start+end/2</code>) and save that into a separate pyranges object. There will be additional processing associated based on what ranges already exist nearby, so I can't just use a <code>resize</code> function, I have to iterate over the object.</p>
<p>First, how would I go about manually creating a new pyranges object, then inserting one by one new ranges into it?</p>
<p>Second, when I try to iterate over my current pyranges object (using this line <code>for myspot in pr.itergrs(hotspots, strand=False) :</code>), it gives me this error: <code>AttributeError: 'tuple' object has no attribute 'stranded'</code>. This error occurs whether or not I have strands in my pyranges, and whether or not I specify <code>strand=False</code>. Why does this occur?</p>
<hr />
<p>The data I am working with: a pyranges covering hotspots in the human genome (all chr1 for the moment)</p>
<pre><code>+--------------+-----------+-----------+-----------+
| Chromosome   | Start     | End       | Score     |
| (category)   | (int32)   | (int32)   | (int64)   |
|--------------+-----------+-----------+-----------|
| chr1         | 154421820 | 154421821 | 26413     |
| chr1         | 30749469  | 30749470  | 14833     |
| chr1         | 55259503  | 55259504  | 9900      |
| chr1         | 14985097  | 14985098  | 7735      |
| ...          | ...       | ...       | ...       |
| chr1         | 247975751 | 247976497 | 1         |
| chr1         | 247977599 | 247981391 | 1         |
| chr1         | 248018012 | 248019454 | 1         |
| chr1         | 248177298 | 248178490 | 1         |
+--------------+-----------+-----------+-----------+
</code></pre>
<p>The hotspots are sorted by their score, and I'd like to go through and for each one, resize it to its centre, keeping the score, but if there is already a hotspot within a certain range (say +/-5kb), to discard it instead, so that I end up with clear chunks between hotspots. This is the &quot;pseudocode&quot; I was planning to do to implement this:</p>
<pre><code>add centre of hotspot 1 to a new pyranges object
create ranges with score = -1 for 5kb flanking each side
for myhotspot in hotspots:
  resize to centre
  if not overlap with new pyranges object
    add centre
    create flanking regions
discard ranges with score&lt;0
</code></pre>
<p>This may not be the most elegant solution, but I'll never have more than ~2000 ranges to process so it should work fine. Is this implementable with pyranges?</p>

<python><pyranges>
Create new range based from existing range (python)
# Create new range based from existing range (python)

(Asked by: username_0 on 2021-02-24)

I have a pyranges object that I am looking to iterate over. For each range, I would like to "resize" it to its centre (so start+end/2) and save that into a separate pyranges object. There will be additional processing associated based on what ranges already exist nearby, so I can't just use a resize function, I have to iterate over the object.
First, how would I go about manually creating a new pyranges object, then inserting one by one new ranges into it?
Second, when I try to iterate over my current pyranges object (using this line for myspot in pr.itergrs(hotspots, strand=False) :), it gives me this error: AttributeError: 'tuple' object has no attribute 'stranded'. This error occurs whether or not I have strands in my pyranges, and whether or not I specify strand=False. Why does this occur?

The data I am working with: a pyranges covering hotspots in the human genome (all chr1 for the moment)
+--------------+-----------+-----------+-----------+
| Chromosome   | Start     | End       | Score     |
| (category)   | (int32)   | (int32)   | (int64)   |
|--------------+-----------+-----------+-----------|
| chr1         | 154421820 | 154421821 | 26413     |
| chr1         | 30749469  | 30749470  | 14833     |
| chr1         | 55259503  | 55259504  | 9900      |
| chr1         | 14985097  | 14985098  | 7735      |
| ...          | ...       | ...       | ...       |
| chr1         | 247975751 | 247976497 | 1         |
| chr1         | 247977599 | 247981391 | 1         |
| chr1         | 248018012 | 248019454 | 1         |
| chr1         | 248177298 | 248178490 | 1         |
+--------------+-----------+-----------+-----------+

The hotspots are sorted by their score, and I'd like to go through and for each one, resize it to its centre, keeping the score, but if there is already a hotspot within a certain range (say +/-5kb), to discard it instead, so that I end up with clear chunks between hotspots. This is the "pseudocode" I was planning to do to implement this:
add centre of hotspot 1 to a new pyranges object
create ranges with score = -1 for 5kb flanking each side
for myhotspot in hotspots:
  resize to centre
  if not overlap with new pyranges object
    add centre
    create flanking regions
discard ranges with score<0

This may not be the most elegant solution, but I'll never have more than ~2000 ranges to process so it should work fine. Is this implementable with pyranges?

**Tags:** python, pyranges

---

## 2 answers

### Answer by username_1 on 2021-03-01. Score: 1 | Accepted answer

Quick and dirty (not to mention atrociously bad O(n)) solution:
import pyranges as pr
import numpy as np
np.random.seed(42 * 10)

# create large df to test on
gr = pr.random(int(1e5), length=10000, chromsizes={"chr1": 249250621})
gr.Score = np.random.randint(250, size=len(gr))

def remove_worst_scores_until_no_overlap(gr):
    df = gr.df
    old_length = -1
    new_length = len(df)
    while old_length != new_length:
        df = df.drop(["Cluster", "Count"], axis=1, errors="ignore")
        df = pr.PyRanges(df).cluster(count=True).df.sort_values("Score", ascending=False)
        df = df[df.duplicated('Cluster', keep='last') | (df.Count == 1)]
        old_length = new_length
        new_length = len(df)

    return pr.PyRanges(df.drop("Cluster Count".split(), axis=1)).sort()

hotspots_no_overlaps = remove_worst_scores_until_no_overlap(gr)

> **Comments:**
> **username_0**: Quick and dirty, but it works!  

---

### Answer by username_1 on 2021-02-24. Score: 4

Answer to main question:
import pyranges as pr
assert int(pr.__version__.split(".")[2]) >= 93, "pip install pyranges==0.0.93"
import numpy as np
np.random.seed(42 * 10)

# create large df to test on
gr = pr.random(int(1e5), length=1, chromsizes={"chr1": 249250621})
gr.Score = np.random.randint(250, size=len(gr))

result = gr.max_disjoint(slack=5000)

# small example that is easy to validate
gr2 = pr.from_string("""Chromosome Start End Score
chr1 1 2 5
chr1 4999 5000 10
chr1 10001 10002 1
chr1 15000 15001 2
chr1 20100 20101 5
""")

gr2.max_disjoint(slack=5000)

# +--------------+-----------+-----------+-----------+
# | Chromosome   |     Start |       End |     Score |
# | (category)   |   (int32) |   (int32) |   (int64) |
# |--------------+-----------+-----------+-----------|
# | chr1         |         1 |         2 |         5 |
# | chr1         |     10001 |     10002 |         1 |
# | chr1         |     20100 |     20101 |         5 |
# +--------------+-----------+-----------+-----------+
# Unstranded PyRanges object has 3 rows and 4 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.


Could you please post an example of what data you have and what you want to end up with?

Inserting row after row is going to be slow, whether it is pandas or pyranges. It is best to collect all the data and then create a new dataframe or PyRanges.

pr.itergrs is for iterating over multiple pyranges at the same time. It takes a list of pyranges while you seem to have a collection of tuples. The function iterates over every chromosome/strand pair in the list and for it extracts the data from every pyranges. See the code example at the end.


To iterate over a pyranges you can do for k, df in gr: # ....
Anyways, to do what you want:
import pyranges as pr
import numpy as np

f1 = pr.data.f1()
f1.Start = ((f1.Start + f1.End)/2).astype(np.int32) # int is 64 bit
f1.End = f1.Start + 1


I understand I have not done the addititonal processing you speak of, but if you can post an example I'm sure I'd be able to help you :)

Code example pr.itergrs:
import pyranges as pr
f1 = pr.data.f1()
f2 = pr.data.f2()
for k, dfs in pr.itergrs([f1, f2], keys=True):
    print(k)
     for df in dfs:
         print("df:")
         print(df)

('chr1', '+')
df:
  Chromosome  Start  End       Name  Score Strand
0       chr1      3    6  interval1      0      +
2       chr1      8    9  interval3      0      +
df:
  Chromosome  Start  End Name  Score Strand
0       chr1      1    2    a      0      +
('chr1', '-')
df:
  Chromosome  Start  End       Name  Score Strand
1       chr1      5    7  interval2      0      -
df:
  Chromosome  Start  End Name  Score Strand
1       chr1      6    7    b      0      -

> **Comments:**
> **username_0**: Thanks for the pointers, this already helps a lot. I'll post my code in the morning (it's 7PM for me now)  
> **username_0**: Added the information you asked for :)  
> **username_0**: Thanks for the solution, do you think you could comment the function a bit more so I can be sure of what each part of it does? From what I can tell, you sort and cluster the df by score in f, then merge all ranges withing 5kb of each other, is that correct?  
> **username_1**: I collect all ranges within 5k and put them in separate clusters. Then the function f in apply merely picks out the row within each cluster with the highest score and discards the rest. I might have made a mistake in that if you have hotspots at 0, 4000, 8000 and 12000 then you would probably want hotspot 0 and 8000 while I discard all but one. Will think about the problem some more and update the answer tomorrow.  
> **username_0**: Yes I was wondering if it behaved like that. That's why I was thinking of creating "flanking" regions to better define which hotspots to discard, but if you find a more elegant solution, I'm all for it!  
> **username_1**: I've now added a method `max_disjoint` to PyRanges in 0.0.93. Thanks for making me realise there was a fundamental interval operation I was missing.  
> **username_0**: Surely "Dev adds new function to package because of you" should be a medal on bioinformatics SE, right? Thank you for your help in any case! I may be calling upon your services later depending on how the project advances  
> **username_1**: I do this for fun, but thanks. Feel free to accept the answer if everything worked as expected.  
> **username_0**: So the `max_disjoint()` works great but when I combine it with the centering you gave, it throws an error: `ValueError: Buffer dtype mismatch, expected 'const int32_t' but got 'long long'`. Do you want the full traceback?  
> **username_1**: Will fix. Could you post an issue on the pyranges github?  
> **username_0**: Thank you for all your help, but `max_disjoint()` has a behaviour that could be improved on I think: for now it seems to start at the first range on the chromosome, whereas it would be useful to be able to specify a column for it to sort on, like the score, before searching for overlaps.  
> **username_0**: Actually I'm not entirely sure what it uses to sort the pyranges on, seeing as even in your example, it's kept the `10001;10002;1` range as the centre of the second cluster  
> **username_1**: max_disjoint() is actually a well-known greedy algorithm: https://en.wikipedia.org/wiki/Maximum_disjoint_set#1-dimensional_intervals:_exact_polynomial_algorithm It starts by sorting on the end column. I will continue thinking about the problem  

---

