Mercurial > repos > jeremyjliu > region_motif_enrichment

comparison region_motif_lib/regions.cpp @ 0:5c044273554d draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
initial commit
author jeremyjliu
date Tue, 05 Aug 2014 13:56:22 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:5c044273554d
1 #include <iostream>
2 #include <vector>
3 #include <algorithm>
4 using namespace std;
5
6 extern "C" {
7 typedef pair <int,int> Se_t;
8 bool se_lt (const Se_t &l,const Se_t &r) { return l.first < r.first; }
9
10 void merge_regions(int *regions, int*nregionsR,int *merge_sepR) {
11 int nregs=nregionsR[0];
12 if(nregs==0) return;
13 int sep=merge_sepR[0];
14 if(sep<1) sep=1;
15 vector<Se_t> reg(nregs);
16 for(int ireg=0;ireg<nregs;ireg++) {
17 reg[ireg]=make_pair(regions[ireg*2],regions[ireg*2+1]);
18 }
19 sort(reg.begin(),reg.end(),se_lt);
20 int *reg_index = new int[nregs];
21 for(int ireg=1;ireg<nregs;ireg++) reg_index[ireg]=-1;
22 reg_index[0]=0;
23 int last_ireg=0;
24 int counter=1;
25 for(int ireg=1;ireg<nregs;ireg++) {
26 if(reg[ireg].first<=reg[last_ireg].second+sep) {
27 if(reg[ireg].second>reg[last_ireg].second) reg[last_ireg].second=reg[ireg].second;
28 } else {
29 last_ireg=ireg;
30 reg_index[counter]=ireg;
31 counter++;
32 }
33 }
34 for(int ireg=0;ireg<counter;ireg++) {
35 regions[ireg*2] = reg[reg_index[ireg]].first;
36 regions[ireg*2+1] = reg[reg_index[ireg]].second;
37 }
38 nregionsR[0]=counter;
39 delete [] reg_index;
40 }
41 void region_minus_region(int *regions, int*nregionsR,int *region2s, int*nregion2sR,int *updatedregions) {
42 int sep=1;
43 merge_regions(regions,nregionsR,&sep);
44 merge_regions(region2s,nregion2sR,&sep);
45 int nregs=nregionsR[0];
46 int nreg2s=nregion2sR[0];
47 for(int i=0;i<2*(nregs+nreg2s);i++) updatedregions[i]=-1;
48 if(nregs==0) return;
49 int ireg = 0;
50 int iregout = 0;
51 for(int ireg2=0; ireg2<nreg2s;ireg2++) {
52 if(ireg==nregs) break;
53 if(region2s[ireg2*2+1] < regions[2*ireg]) continue;
54 if(region2s[ireg2*2] > regions[2*ireg+1]) {
55 updatedregions[2*iregout] = regions[2*ireg];
56 updatedregions[2*iregout+1] = regions[2*ireg+1];
57 ireg++;
58 ireg2--;
59 iregout++;
60 continue;
61 }
62 int s = regions[ireg*2];
63 int e = regions[ireg*2+1];
64 int s2 = region2s[ireg2*2];
65 int e2 = region2s[ireg2*2+1];
66 if(s2<=s && e2>=e) {
67 ireg++;
68 ireg2--;
69 }
70 else if(s2<=s) {
71 regions[ireg*2] = e2+1;
72 continue;
73 } else if(e2>=e) {
74 updatedregions[2*iregout] = s;
75 updatedregions[2*iregout+1] = s2-1;
76 ireg2--;
77 iregout++;
78 ireg++;
79 } else {
80 updatedregions[2*iregout] = s;
81 updatedregions[2*iregout+1] = s2-1;
82 regions[ireg*2] = e2+1;
83 iregout++;
84 ireg2--;
85 }
86 }
87 while(ireg<nregs) {
88 updatedregions[2*iregout] = regions[2*ireg];
89 updatedregions[2*iregout+1] = regions[2*ireg+1];
90 ireg++;
91 iregout++;
92 }
93 }
94 void intersection_of_regions(int *regions, int*nregionsR,int *region2s, int*nregion2sR,int *updatedregions) {
95 int sep=1;
96 merge_regions(regions,nregionsR,&sep);
97 merge_regions(region2s,nregion2sR,&sep);
98 int nregs=nregionsR[0];
99 int nreg2s=nregion2sR[0];
100 for(int i=0;i<2*(nregs+nreg2s);i++) updatedregions[i]=-1;
101 if(nregs==0) return;
102 if(nreg2s==0) return;
103 int ireg2 = 0;
104 int iregout = 0;
105 for(int ireg=0; ireg<nregs;ireg++) {
106 if(ireg2==nreg2s) return;
107 if(regions[ireg*2+1] < region2s[2*ireg2]) continue;
108 if(regions[ireg*2] > region2s[2*ireg2+1]) {ireg2++; ireg--; continue;}
109
110 int s = regions[ireg*2];
111 if(s<region2s[ireg2*2]) s = region2s[ireg2*2];
112 int e = regions[ireg*2+1];
113 if(e>region2s[ireg2*2+1]) {
114 e = region2s[ireg2*2+1];
115 ireg2++;
116 ireg--;
117 }
118 updatedregions[2*iregout] = s;
119 updatedregions[2*iregout+1] = e;
120 iregout++;
121 }
122 }
123 }