Sie erfinden bedtools intersect (oder bedops) neu, für die es bereits ein praktisches Python-Modul gibt:

  von pybedtools importieren BedTools3 = BedTool ('s3.bed') s4 = BedTool ('s4.bed') print (s4.intersect (s3, wa = True, wb = True, F = 1) )  

Das wb = True entspricht -wb , wobei sich Bedtools in der Befehlszeile schneiden. In ähnlicher Weise ist F = 1 dasselbe wie -F 1 .

Bezüglich Ihres Python-Codes

Wenn Sie die experimentellen Bereiche möchten, die vollständig in einem der Referenzbereiche enthalten sind , müssen Sie die Koordinaten in der folgenden Reihenfolge haben:

  cst < = tst < ten < = cen  

Wenn Sie möchten, sind die experimentellen Bereiche, die einen der Referenzbereiche strong überlappen > muss entweder der Anfang oder das Ende des experimentellen Bereichs innerhalb des Referenzbereichs liegen:

  (cst < = tst < cen) oder (cst < ten < = cen)  

Ihr Code entspricht keiner der beiden Möglichkeiten:

  cst > = tst und cen < = ten  

Dies ist äquivalent zu:

  tst < = cst und cen < = zehn  

(oder tst < = cst < cen < = ten , da wir wissen, dass cst < cen , per Definition eines Bettes inter val).

Mit diesem Umschreiben können wir leichter erkennen, dass Sie tatsächlich die experimentellen Bereiche auswählen, die einen Referenzbereich enthalten.

Hier einige (python3) Code, der Ihnen die Ergebnisse in den beiden anderen Situationen liefert:

  #! / usr / bin / env python3ref_intervals = [] mit open ("s4.txt", "r") als f: für Zeile in f: (chr, start, end) = line.strip (). split ("\ t") ref_intervals.append ((chr, int (start), int (end))) exp_intervals = [ ] mit open ("s3.txt", "r") als f: für Zeile in f: (chr, start, end) = line.strip (). split ("\ t") exp_intervals.append ((chr, int (start), int (end))) enthalten = [] überlappend = [] für (r_chr, r_start, r_end) in ref_intervals: für (e_chr, e_start, e_end) in exp_intervals: wenn e_chr == r_chr: wenn r_start < = e_start < r_end oder r_start < e_end < = r_end: overlapping.append ((e_chr, e_start, e_end))
wenn r_start < = e_start < e_end < = r_end: enthält "\ t") print ("enthalten") für (chr, start, end) in enthalten: print (chr, start, end, sep = "\ t")  

Wenn ich Wenn Sie es ausführen, erhalte ich die folgenden Ergebnisse:

  $ ./overlap.py überlappend1 10 201 5 201 10 501 14 172 20 302 25 302 20 60contained1 10 201 14 172 20 302 25 30  

Es ist wahrscheinlich eine gute Übung, dies zu programmieren, aber wie die anderen Antworten zeigen, gibt es effiziente Tools, die in einer professionellen Umgebung eine bessere Lösung darstellen würden .

Sie können stattdessen BEDOPS verwenden:

  $ sort-bed s3.txt > s3.bed $ sort-bed s4.txt > s4.bed $ bedops --element-of 1 s3.bed s4.bed > answer.bed  

Wenn Sie es in Python ausführen müssen, können Sie subprocess.check_output () :

  Unterprozess importieren ... try: result = subprocess.check_output ("bedops --element-von 1% s% s >% s"% (set_a_fn, set_b_fn, answer_fn ), shell = True) außer subprocess.CalledProcessError als err: erhöhe SystemExit ("Bedops konnten nicht ausgeführt werden \ n") # mache Sachen mit 'result'  

pyranges Antwort:

  # pip install pyranges # oder # conda install -c bioconda pyrangesimport pyranges as prs3 = pr.read_bed ("s3.bed") s4 = pr.read_bed ("s4.bed") s3.intersect (s4, how = "Containment")  

Antwort mit Setup (Bearbeiten):

  Pandas als pd aus importieren StringIOimport pyranges als prc1 = "" 1 10 201 5 202 20 302 25 301 10 502 20 601 14 17 "" c2 = "" 1 10 202 20 30 "" column = "Chromosom Start End" .split () df1 = pd.read_table (StringIO (c1), sep = "\ s +", header = Keine, Namen = Spalten) df2 = pd.read_table (StringIO (c2), sep = "\ s +", header = Keine, Namen = Spalten) gr1 = pr.PyRanges (df1) gr2 = pr.PyRanges (df2) print (gr1.intersect (gr2, how = "Containment"))  

Ergebnis:

  + -------------- + ----------- + ----- ------ + | Chromosom | Start | Ende || (Kategorie) | (int32) | (int32) || -------------- + ----------- + ----------- || 1 | 10 | 20 || 1 | 14 | 17 || 2 | 20 | 30 || 2 | 25 | 30 | + -------------- + ----------- + ----------- + Unstranded PyRanges-Objekt hat 4 Zeilen und 3 Spalten aus 2 Chromosomen. Zum Drucken wurden die PyRanges nach Chromosomen sortiert.