Answers.Body: <p><a href="https://github.com/biocore-ntnu/pyranges" rel="nofollow noreferrer">pyranges</a> answer:</p>
  
  <pre><code># pip install pyranges
  # or 
  # conda install -c bioconda pyranges
  
  import pyranges as pr
  s3 = pr.read_bed("s3.bed")
  s4 = pr.read_bed("s4.bed")
  s3.intersect(s4, how="containment")
  </code></pre>
  
  <hr>
  
  <p>Answer with setup (Edit):</p>
  
  <pre><code>import pandas as pd
  from io import StringIO
  import pyranges as pr
  
  c1 = """1   10  20
  1   5   20
  2   20  30
  2   25  30
  1   10  50
  2   20  60
  1   14  17"""
  
  c2 = """1   10  20
  2   20  30"""
  
  columns = "Chromosome Start End".split()
  
  df1 = pd.read_table(StringIO(c1), sep="\s+", header=None, names=columns)
  df2 = pd.read_table(StringIO(c2), sep="\s+", header=None, names=columns)
  
  gr1 = pr.PyRanges(df1)
  gr2 = pr.PyRanges(df2)
  
  print(gr1.intersect(gr2, how="containment"))
  </code></pre>
  
  <p>Result:</p>
  
  <pre><code>+--------------+-----------+-----------+
  |   Chromosome |     Start |       End |
  |   (category) |   (int32) |   (int32) |
  |--------------+-----------+-----------|
  |            1 |        10 |        20 |
  |            1 |        14 |        17 |
  |            2 |        20 |        30 |
  |            2 |        25 |        30 |
  +--------------+-----------+-----------+
  Unstranded PyRanges object has 4 rows and 3 columns from 2 chromosomes.
  For printing, the PyRanges was sorted on Chromosome.
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-10-29T11:03:22.650Z
Answers.LastActivityDate: 2019-10-30T08:15:23.480Z
Answers.LastEditDate: 2019-10-30T08:15:23.480Z
Answers.UserDisplayName: The Unfun Cat
Answers.Body: <p>You're reinventing <code>bedtools intersect</code> (or bedops), for which there's <a href="https://daler.github.io/pybedtools/" rel="nofollow noreferrer">already a convenient python module</a>:</p>
  
  <pre class="lang-python prettyprint-override"><code>from pybedtools import BedTool
  
  s3 = BedTool('s3.bed')
  s4 = BedTool('s4.bed')
  
  print(s4.intersect(s3, wa=True, wb=True, F=1))
  </code></pre>
  
  <p>The <code>wb=True</code> is equivalent to <code>-wb</code> with <code>bedtools intersect</code> on the command line. Similarly, <code>F=1</code> is the same as <code>-F 1</code>.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-20T19:55:58.207Z
Answers.Text: `NotImplementedError: "intersectBed" does not appear to be installed or on the path, so this method is disabled.  Please install a more recent version of BEDTools and re-import to use this method.` used pip install pybedtools
Answers.UserDisplayName: novicebioinforesearcher
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-20T19:57:03.333Z
Answers.Text: You should be able to read and understand how to deal with that error message.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-20T20:15:13.440Z
Answers.Text: Hi i installed bedtools  made is environment variable i can use it in bash from any directory i still get the same error, is there a different way than this? thank you
Answers.UserDisplayName: novicebioinforesearcher
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-20T20:18:34.027Z
Answers.Text: You can always install it via conda, though it works for me with the outdated version of bedtools that I have installed (2.25.0).
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-20T20:24:06.300Z
Answers.Text: i was able to get it working thanks, but the result seems to be ambiguous have added it to my original post, 5-20 seems to be missing
Answers.UserDisplayName: novicebioinforesearcher
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-20T20:26:10.050Z
Answers.Text: Whether you want `F=1` or `f=1` will depend on your exact needs and the order you do things.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-20T19:38:59.547Z
Answers.LastActivityDate: 2017-06-22T11:22:23.687Z
Answers.LastEditDate: 2017-06-22T11:22:23.687Z
Answers.UserDisplayName: Devon Ryan
Answers.Body: <h3>Regarding your python code</h3>
  
  <p>If you want the <strong>experimental ranges that are entirely contained in one of the reference ranges</strong>, you need to have the coordinates in the following order:</p>
  
  <pre><code>cst &lt;= tst &lt; ten &lt;= cen
  </code></pre>
  
  <p>If what you want are the <strong>experimental ranges that overlap one of the reference ranges</strong>, you need to have either the start or the end of the experimental range fall within the reference range:</p>
  
  <pre><code>(cst &lt;= tst &lt; cen) or (cst &lt; ten &lt;= cen)
  </code></pre>
  
  <p>Your code is equivalent to neither possibilities:</p>
  
  <pre><code>cst &gt;= tst and cen &lt;= ten
  </code></pre>
  
  <p>This is equivalent to:</p>
  
  <pre><code>tst &lt;= cst and cen &lt;= ten
  </code></pre>
  
  <p>(Or <code>tst &lt;= cst &lt; cen &lt;= ten</code>, since we know that <code>cst &lt; cen</code>, by definition of a bed interval).</p>
  
  <p>With this rewriting, we can more easily see that you are actually selecting the <strong>experimental ranges that contain a reference range</strong>.</p>
  
  <p>Here is some (python3) code that gives you the results in the other two situations:</p>
  
  <pre><code>#!/usr/bin/env python3
  
  ref_intervals = []
  with open("s4.txt", "r") as f:
      for line in f:
          (chr, start, end) = line.strip().split("\t")
          ref_intervals.append((chr, int(start), int(end)))
  
  exp_intervals = []
  with open("s3.txt", "r") as f:
      for line in f:
          (chr, start, end) = line.strip().split("\t")
          exp_intervals.append((chr, int(start), int(end)))
  
  contained = []
  overlapping = []
  for (r_chr, r_start, r_end) in ref_intervals:
      for (e_chr, e_start, e_end) in exp_intervals:
          if e_chr == r_chr:
              if r_start &lt;= e_start &lt; r_end or r_start &lt; e_end &lt;= r_end:
                  overlapping.append((e_chr, e_start, e_end))
              if r_start &lt;= e_start &lt; e_end &lt;= r_end:
                  contained.append((e_chr, e_start, e_end))
  
  print("overlapping")
  for (chr, start, end) in overlapping:
      print(chr, start, end, sep="\t")
  
  print("contained")
  for (chr, start, end) in contained:
      print(chr, start, end, sep="\t")
  </code></pre>
  
  <p>If I run it, I obtain the following results:</p>
  
  <pre><code>$ ./overlap.py 
  overlapping
  1   10  20
  1   5   20
  1   10  50
  1   14  17
  2   20  30
  2   25  30
  2   20  60
  contained
  1   10  20
  1   14  17
  2   20  30
  2   25  30
  </code></pre>
  
  <p>It is probably an good exercise to program this, but as the other answers point out, <strong>there are efficient tools that would be a better solution in a professional setting</strong>.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-21T15:57:15.767Z
Answers.Text: this makes sense, I guess I need to start thinking and talking like a programmer which would help me write the questions like one.  Sure there are other tools but doing it myself helps me to learn and understand whats exactly going behind the tools.
Answers.UserDisplayName: novicebioinforesearcher
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-21T15:27:07.430Z
Answers.LastActivityDate: 2017-06-21T15:34:15.073Z
Answers.LastEditDate: 2017-06-21T15:34:15.073Z
Answers.UserDisplayName: bli
Answers.Body: <p>You could use <a href="http://bedops.readthedocs.io/en/latest/#" rel="nofollow noreferrer">BEDOPS</a>, instead:</p>
  
  <pre><code>$ sort-bed s3.txt &gt; s3.bed
  $ sort-bed s4.txt &gt; s4.bed
  $ bedops --element-of 1 s3.bed s4.bed &gt; answer.bed
  </code></pre>
  
  <p>If you need to run it from within Python, you could use <code>subprocess.check_output()</code>:</p>
  
  <pre><code>import subprocess
  ...
  try:
      result = subprocess.check_output("bedops --element-of 1 %s %s &gt; %s" % (set_a_fn, set_b_fn, answer_fn), shell=True)
  except subprocess.CalledProcessError as err:
      raise SystemExit("Could not run bedops\n")
  # do stuff with 'result'
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-20T19:37:07.933Z
Answers.Text: thank you I would certainly try it , but can you please help me understand what is wrong with my logic/code
Answers.UserDisplayName: novicebioinforesearcher
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-20T19:38:42.280Z
Answers.Text: You're better off not reinventing the wheel, is my advice. Python is not the right tool for set operations.
Answers.UserDisplayName: Alex Reynolds
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-20T19:39:34.730Z
Answers.Text: If you really want to do this, use debug statements (`sys.stderr.write(...)`) so you can see what's happening at each stage.
Answers.UserDisplayName: Alex Reynolds
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-20T19:46:52.830Z
Answers.Text: If this is a homework assignment and you are required to write your own solution from scratch, please add a "homework" tag to your question to make this explicit. We're happy to help with homework questions, but the responses will likely be different (and hopefully accommodating).
Answers.UserDisplayName: Daniel Standage
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-20T19:47:18.783Z
Answers.Text: If it's not a homework assignment, Alex is right: better to use well-tested tools for interval arithmetic.
Answers.UserDisplayName: Daniel Standage
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-20T20:25:26.347Z
Answers.Text: @DanielStandage yes it is homework will make changes to tag, tried doing it it would not allow me to create a tag
Answers.UserDisplayName: novicebioinforesearcher
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-20T20:30:21.943Z
Answers.Text: @novicebioinforesearcher I've added the homework tag to your post.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-20T20:40:02.237Z
Answers.Text: @AlexReynolds does the method suggested account for strand as well?
Answers.UserDisplayName: novicebioinforesearcher
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-20T20:55:14.067Z
Answers.Text: I'm not sure. I don't see anything in your code that separates intervals by strand. You could use `awk ($6=="+")` etc. to do this.
Answers.UserDisplayName: Alex Reynolds
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-20T21:48:41.443Z
Answers.Text: @AlexReynolds sorry for late reply i was curious you are right my code does not handle strand, bedops is new for me hence wanted to ask you upfront will install and update you with results
Answers.UserDisplayName: novicebioinforesearcher
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-20T19:35:59.390Z
Answers.LastActivityDate: 2017-06-20T19:35:59.390Z
Answers.UserDisplayName: Alex Reynolds
Body: <p>I have two files </p>
  
  <pre><code>s3.txt :
  1   10  20
  1   5   20
  2   20  30
  2   25  30
  1   10  50
  2   20  60
  1   14  17
  
  s4.txt:
  1   10  20
  2   20  30
  </code></pre>
  
  <p>I am trying to match col0 of both the files and get rows that fall between range(inclusive of themselves) 10-20 and 20-30 as seen in s4 file.
  file s4 has co ordinates which can be used as reference range (chrom start and end) and s3 has list of co ordinates from an experimental condition what I am trying to achieve is to which co-ordinates from my file s3 fall on or between my reference co ordinates in s4.  </p>
  
  <p>code so far:</p>
  
  <pre><code>containing_ranges = []
  with open('s4.txt', 'r') as f:
      for line in f:
          fields = line.strip().split('\t')
          containing_ranges.append(fields)
  
  tested_ranges = []       
  with open('s3.txt', 'r') as f:
      for line in f:
          fields = line.strip().split('\t')
          tested_ranges.append(fields)
  
  for c_range in containing_ranges:
      for t_range in tested_ranges:
          tst = int(t_range[1])
          ten = int(t_range[2])
          cst = int(c_range[1])
          cen = int(c_range[2])
          if  c_range[0] == t_range[0]:
              included = cst &gt;= tst and cen &lt;= ten
              if included == True:
                 print t_range
  </code></pre>
  
  <p>Output with missing row(1   14  17) : </p>
  
  <pre><code>['1', '10', '20']
  ['1', '5', '20']
  ['1', '10', '50']
  ['2', '20', '30']
  ['2', '20', '60']
  </code></pre>
  
  <p>Desired output:</p>
  
  <pre><code>1   10  20
  2   20  30
  2   25  30
  1   14  17
  </code></pre>
  
  <p>Not sure if my logic is wrong and why does it miss 14-17 as it falls between 10-20</p>
  
  <pre><code>[EDIT] using pybedtools
  
  &gt;&gt;&gt; print(s4.intersect(s3, wb=True))
  1   10  20  1   10  20
  1   10  20  1   5   20
  1   10  20  1   10  50
  1   14  17  1   14  17
  2   20  30  2   20  30
  2   25  30  2   25  30
  2   20  30  2   20  60
  
  &gt;&gt;&gt; print(s4.intersect(s3, wa=True, wb=True, F=1))
  1   10  20  1   10  20
  1   10  20  1   14  17
  2   20  30  2   20  30
  2   20  30  2   25  30
  
  
  using bedops 
  bin$ less answer.bed 
  1       5       20
  1       10      20
  1       10      50
  1       14      17
  2       20      30
  2       20      60
  2       25      30
  
  using @bli code(on python2.7)
  ('1', 10, 20)
  ('1', 14, 17)
  ('2', 20, 30)
  ('2', 25, 30)
   why can I not see the interval 1 5 20
  </code></pre>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-20T20:26:08.323Z
Comments.Text: What result do you get with `bedops --element-of`?
Comments.UserDisplayName: Alex Reynolds
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-20T23:04:06.730Z
Comments.Text: Please bring your logic/coding questions to Stack Overflow. The relationship of your question to the subject of Bioinformatics is merely coincidental.
Comments.UserDisplayName: Robert Cartaino
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-20T23:50:04.870Z
Comments.Text: @RobertC If OP adds a "bed" tag, this question will immediately look like a bioinformatics question. Also, see the answers. OP is much more likely to get such spot-on answers here. This question could be improved for sure, but it is not off-topic.
Comments.UserDisplayName: user172818
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-21T04:00:09.483Z
Comments.Text: Just use bedops, as indicated. Using wrappers to command-line tools is rarely a substitute for learning the tools.
Comments.UserDisplayName: Alex Reynolds
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-21T04:49:57.843Z
Comments.Text: Can you please add more story/context around this question? It has the appearance of a pure programming question (presumably why it was marked as off-topic). It would be nice if you could explain what the different numbers mean, and why you want to do this.
Comments.UserDisplayName: gringer
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-21T08:39:07.003Z
Comments.Text: You should use more meaningful names for your variables. It would make the code easier to read, for others but also for you.
Comments.UserDisplayName: bli
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-21T09:18:48.443Z
Comments.Text: It seems you want the ranges from `s3.txt` ("tested ranges") that are included in a range in `s4.txt` ("containing ranges"). In this case, I think your error is in the start and end coordinates comparisons. Let's note `t_start` and `t_end` the coordinates of the tested range and `c_start` and `c_end` the coordinates of the containing range. What you want is `c_start <= t_start and t_end <= c_end`.
Comments.UserDisplayName: bli
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-21T09:31:36.670Z
Comments.Text: I edited your question to use more meaningful variable names, and I also removed `1    5    20` from the expected output: If I understood correctly, this is not what you want because this is not included in any of the ranges defined in `s4.txt`
Comments.UserDisplayName: bli
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-21T09:41:13.227Z
Comments.Text: I can't post an answer since your question is "on hold", but here is a (hopefully) corrected version of your code, with minor some coding style improvements and using python3: http://paste.ubuntu.com/24915950/
  I hope this helps.
Comments.UserDisplayName: bli
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-21T14:09:21.757Z
Comments.Text: ohh I was away for a bit, so many suggestions thank you all. will edit my post once I go through each of your suggestions
Comments.UserDisplayName: novicebioinforesearcher
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-21T14:47:56.827Z
Comments.Text: @AlexReynolds added answer
Comments.UserDisplayName: novicebioinforesearcher
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-21T14:48:44.160Z
Comments.Text: @bli thank you for cleaning up code added answer
Comments.UserDisplayName: novicebioinforesearcher
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-21T14:56:40.477Z
Comments.Text: You say "which co-ordinates from my file s3 fall on or between my reference co ordinates in s4". If I interpret this correctly, this means that you accept partial overlaps too, not only complete inclusions. Then the desired output should be all the ranges in s3, and not the restricted list I mistakenly corrected.
Comments.UserDisplayName: bli
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-21T20:45:03.940Z
Comments.Text: @novicebioinforesearcher It looks like `bedops` was able to find your missing interval. If you need to deal with strand labels in the sixth column (per BED specification), you can split a BED file by strand via `awk '$6=="+"' in.bed > in.forward.bed` and `awk '$6=="-"' in.bed > in.reverse.bed`, and then run set operations on each of the strand-split files. If you need to reconstruct one file at the end, use `bedops -u` to do a multiset union of all input BED files.
Comments.UserDisplayName: Alex Reynolds
Tags: <bed><genomics><python><homework>
Title: Range overlap python error with genomic regions
ThreadText: # Range overlap python error with genomic regions
  
  (Asked by: username_0 on 2017-06-20)
  
  I have two files
  
  ```
  s3.txt :
  1   10  20
  1   5   20
  2   20  30
  2   25  30
  1   10  50
  2   20  60
  1   14  17
  
  s4.txt:
  1   10  20
  2   20  30
  
  ```
  
  I am trying to match col0 of both the files and get rows that fall between range(inclusive of themselves) 10-20 and 20-30 as seen in s4 file.
  file s4 has co ordinates which can be used as reference range (chrom start and end) and s3 has list of co ordinates from an experimental condition what I am trying to achieve is to which co-ordinates from my file s3 fall on or between my reference co ordinates in s4.
  
  code so far:
  
  ```
  containing_ranges = []
  with open('s4.txt', 'r') as f:
      for line in f:
          fields = line.strip().split('\t')
          containing_ranges.append(fields)
  
  tested_ranges = []       
  with open('s3.txt', 'r') as f:
      for line in f:
          fields = line.strip().split('\t')
          tested_ranges.append(fields)
  
  for c_range in containing_ranges:
      for t_range in tested_ranges:
          tst = int(t_range[1])
          ten = int(t_range[2])
          cst = int(c_range[1])
          cen = int(c_range[2])
          if  c_range[0] == t_range[0]:
              included = cst >= tst and cen <= ten
              if included == True:
                 print t_range
  
  ```
  
  Output with missing row(1 14 17) :
  
  ```
  ['1', '10', '20']
  ['1', '5', '20']
  ['1', '10', '50']
  ['2', '20', '30']
  ['2', '20', '60']
  
  ```
  
  Desired output:
  
  ```
  1   10  20
  2   20  30
  2   25  30
  1   14  17
  
  ```
  
  Not sure if my logic is wrong and why does it miss 14-17 as it falls between 10-20
  
  ```
  [EDIT] using pybedtools
  
  >>> print(s4.intersect(s3, wb=True))
  1   10  20  1   10  20
  1   10  20  1   5   20
  1   10  20  1   10  50
  1   14  17  1   14  17
  2   20  30  2   20  30
  2   25  30  2   25  30
  2   20  30  2   20  60
  
  >>> print(s4.intersect(s3, wa=True, wb=True, F=1))
  1   10  20  1   10  20
  1   10  20  1   14  17
  2   20  30  2   20  30
  2   20  30  2   25  30
  
  using bedops 
  bin$ less answer.bed 
  1       5       20
  1       10      20
  1       10      50
  1       14      17
  2       20      30
  2       20      60
  2       25      30
  
  using @username_5 code(on python2.7)
  ('1', 10, 20)
  ('1', 14, 17)
  ('2', 20, 30)
  ('2', 25, 30)
   why can I not see the interval 1 5 20
  
  ```
  
  **Tags:** bed, genomics, python, homework
  
  > **Comments:**
  > **username_1**: What result do you get with `bedops --element-of`?  
  > **username_2**: Please bring your logic/coding questions to Stack Overflow. The relationship of your question to the subject of Bioinformatics is merely coincidental.  
  > **username_3**: @RobertC If OP adds a "bed" tag, this question will immediately look like a bioinformatics question. Also, see the answers. OP is much more likely to get such spot-on answers here. This question could be improved for sure, but it is not off-topic.  
  > **username_1**: Just use bedops, as indicated. Using wrappers to command-line tools is rarely a substitute for learning the tools.  
  > **username_4**: Can you please add more story/context around this question? It has the appearance of a pure programming question (presumably why it was marked as off-topic). It would be nice if you could explain what the different numbers mean, and why you want to do this.  
  > **username_5**: You should use more meaningful names for your variables. It would make the code easier to read, for others but also for you.  
  > **username_5**: It seems you want the ranges from `s3.txt` ("tested ranges") that are included in a range in `s4.txt` ("containing ranges"). In this case, I think your error is in the start and end coordinates comparisons. Let's note `t_start` and `t_end` the coordinates of the tested range and `c_start` and `c_end` the coordinates of the containing range. What you want is `c_start <= t_start and t_end <= c_end`.  
  > **username_5**: I edited your question to use more meaningful variable names, and I also removed `1    5    20` from the expected output: If I understood correctly, this is not what you want because this is not included in any of the ranges defined in `s4.txt`  
  > **username_5**: I can't post an answer since your question is "on hold", but here is a (hopefully) corrected version of your code, with minor some coding style improvements and using python3: http://paste.ubuntu.com/24915950/ I hope this helps.  
  > **username_0**: ohh I was away for a bit, so many suggestions thank you all. will edit my post once I go through each of your suggestions  
  > **username_0**: @username_1 added answer  
  > **username_0**: @username_5 thank you for cleaning up code added answer  
  > **username_5**: You say "which co-ordinates from my file s3 fall on or between my reference co ordinates in s4". If I interpret this correctly, this means that you accept partial overlaps too, not only complete inclusions. Then the desired output should be all the ranges in s3, and not the restricted list I mistakenly corrected.  
  > **username_1**: @username_0 It looks like `bedops` was able to find your missing interval. If you need to deal with strand labels in the sixth column (per BED specification), you can split a BED file by strand via `awk '$6=="+"' in.bed > in.forward.bed` and `awk '$6=="-"' in.bed > in.reverse.bed`, and then run set operations on each of the strand-split files. If you need to reconstruct one file at the end, use `bedops -u` to do a multiset union of all input BED files.  
  
  ---
  
  ## 4 answers
  
  ### Answer by username_6 on 2019-10-29. Score: 1 | Accepted answer
  
  [pyranges](https://github.com/biocore-ntnu/pyranges) answer:
  
  ```
  # pip install pyranges
  # or 
  # conda install -c bioconda pyranges
  
  import pyranges as pr
  s3 = pr.read_bed("s3.bed")
  s4 = pr.read_bed("s4.bed")
  s3.intersect(s4, how="containment")
  
  ```
  
  ---
  
  Answer with setup (Edit):
  
  ```
  import pandas as pd
  from io import StringIO
  import pyranges as pr
  
  c1 = """1   10  20
  1   5   20
  2   20  30
  2   25  30
  1   10  50
  2   20  60
  1   14  17"""
  
  c2 = """1   10  20
  2   20  30"""
  
  columns = "Chromosome Start End".split()
  
  df1 = pd.read_table(StringIO(c1), sep="\s+", header=None, names=columns)
  df2 = pd.read_table(StringIO(c2), sep="\s+", header=None, names=columns)
  
  gr1 = pr.PyRanges(df1)
  gr2 = pr.PyRanges(df2)
  
  print(gr1.intersect(gr2, how="containment"))
  
  ```
  
  Result:
  
  ```
  +--------------+-----------+-----------+
  |   Chromosome |     Start |       End |
  |   (category) |   (int32) |   (int32) |
  |--------------+-----------+-----------|
  |            1 |        10 |        20 |
  |            1 |        14 |        17 |
  |            2 |        20 |        30 |
  |            2 |        25 |        30 |
  +--------------+-----------+-----------+
  Unstranded PyRanges object has 4 rows and 3 columns from 2 chromosomes.
  For printing, the PyRanges was sorted on Chromosome.
  
  ```
  
  
  ---
  
  ### Answer by username_7 on 2017-06-20. Score: 5
  
  You're reinventing `bedtools intersect` (or bedops), for which there's [already a convenient python module](https://daler.github.io/pybedtools/):
  
  ```python
  from pybedtools import BedTool
  
  s3 = BedTool('s3.bed')
  s4 = BedTool('s4.bed')
  
  print(s4.intersect(s3, wa=True, wb=True, F=1))
  
  ```
  
  The `wb=True` is equivalent to `-wb` with `bedtools intersect` on the command line. Similarly, `F=1` is the same as `-F 1`.
  
  > **Comments:**
  > **username_0**: `NotImplementedError: "intersectBed" does not appear to be installed or on the path, so this method is disabled.  Please install a more recent version of BEDTools and re-import to use this method.` used pip install pybedtools  
  > **username_7**: You should be able to read and understand how to deal with that error message.  
  > **username_0**: Hi i installed bedtools  made is environment variable i can use it in bash from any directory i still get the same error, is there a different way than this? thank you  
  > **username_7**: You can always install it via conda, though it works for me with the outdated version of bedtools that I have installed (2.25.0).  
  > **username_0**: i was able to get it working thanks, but the result seems to be ambiguous have added it to my original post, 5-20 seems to be missing  
  > **username_7**: Whether you want `F=1` or `f=1` will depend on your exact needs and the order you do things.  
  
  ---
  
  ### Answer by username_5 on 2017-06-21. Score: 3
  
  ### Regarding your python code
  
  If you want the **experimental ranges that are entirely contained in one of the reference ranges**, you need to have the coordinates in the following order:
  
  ```
  cst <= tst < ten <= cen
  
  ```
  
  If what you want are the **experimental ranges that overlap one of the reference ranges**, you need to have either the start or the end of the experimental range fall within the reference range:
  
  ```
  (cst <= tst < cen) or (cst < ten <= cen)
  
  ```
  
  Your code is equivalent to neither possibilities:
  
  ```
  cst >= tst and cen <= ten
  
  ```
  
  This is equivalent to:
  
  ```
  tst <= cst and cen <= ten
  
  ```
  
  (Or `tst <= cst < cen <= ten`, since we know that `cst < cen`, by definition of a bed interval).
  
  With this rewriting, we can more easily see that you are actually selecting the **experimental ranges that contain a reference range**.
  
  Here is some (python3) code that gives you the results in the other two situations:
  
  ```
  #!/usr/bin/env python3
  
  ref_intervals = []
  with open("s4.txt", "r") as f:
      for line in f:
          (chr, start, end) = line.strip().split("\t")
          ref_intervals.append((chr, int(start), int(end)))
  
  exp_intervals = []
  with open("s3.txt", "r") as f:
      for line in f:
          (chr, start, end) = line.strip().split("\t")
          exp_intervals.append((chr, int(start), int(end)))
  
  contained = []
  overlapping = []
  for (r_chr, r_start, r_end) in ref_intervals:
      for (e_chr, e_start, e_end) in exp_intervals:
          if e_chr == r_chr:
              if r_start <= e_start < r_end or r_start < e_end <= r_end:
                  overlapping.append((e_chr, e_start, e_end))
              if r_start <= e_start < e_end <= r_end:
                  contained.append((e_chr, e_start, e_end))
  
  print("overlapping")
  for (chr, start, end) in overlapping:
      print(chr, start, end, sep="\t")
  
  print("contained")
  for (chr, start, end) in contained:
      print(chr, start, end, sep="\t")
  
  ```
  
  If I run it, I obtain the following results:
  
  ```
  $ ./overlap.py 
  overlapping
  1   10  20
  1   5   20
  1   10  50
  1   14  17
  2   20  30
  2   25  30
  2   20  60
  contained
  1   10  20
  1   14  17
  2   20  30
  2   25  30
  
  ```
  
  It is probably an good exercise to program this, but as the other answers point out, **there are efficient tools that would be a better solution in a professional setting**.
  
  > **Comments:**
  > **username_0**: this makes sense, I guess I need to start thinking and talking like a programmer which would help me write the questions like one.  Sure there are other tools but doing it myself helps me to learn and understand whats exactly going behind the tools.  
  
  ---
  
  ### Answer by username_1 on 2017-06-20. Score: 1
  
  You could use [BEDOPS](http://bedops.readthedocs.io/en/latest/#), instead:
  
  ```
  $ sort-bed s3.txt > s3.bed
  $ sort-bed s4.txt > s4.bed
  $ bedops --element-of 1 s3.bed s4.bed > answer.bed
  
  ```
  
  If you need to run it from within Python, you could use `subprocess.check_output()`:
  
  ```
  import subprocess
  ...
  try:
      result = subprocess.check_output("bedops --element-of 1 %s %s > %s" % (set_a_fn, set_b_fn, answer_fn), shell=True)
  except subprocess.CalledProcessError as err:
      raise SystemExit("Could not run bedops\n")
  # do stuff with 'result'
  
  ```
  
  > **Comments:**
  > **username_0**: thank you I would certainly try it , but can you please help me understand what is wrong with my logic/code  
  > **username_1**: You're better off not reinventing the wheel, is my advice. Python is not the right tool for set operations.  
  > **username_1**: If you really want to do this, use debug statements (`sys.stderr.write(...)`) so you can see what's happening at each stage.  
  > **username_8**: If this is a homework assignment and you are required to write your own solution from scratch, please add a "homework" tag to your question to make this explicit. We're happy to help with homework questions, but the responses will likely be different (and hopefully accommodating).  
  > **username_8**: If it's not a homework assignment, username_1 is right: better to use well-tested tools for interval arithmetic.  
  > **username_0**: @username_8 yes it is homework will make changes to tag, tried doing it it would not allow me to create a tag  
  > **username_7**: @username_0 I've added the homework tag to your post.  
  > **username_0**: @username_1 does the method suggested account for strand as well?  
  > **username_1**: I'm not sure. I don't see anything in your code that separates intervals by strand. You could use `awk ($6=="+")` etc. to do this.  
  > **username_0**: @username_1 sorry for late reply i was curious you are right my code does not handle strand, bedops is new for me hence wanted to ask you upfront will install and update you with results  
  
  ---
  
  