|
0
|
1 function [new_genes]=detect_overlapping_regions(genes);
|
|
|
2 % this function determines regions in a gene which overlapp with
|
|
|
3 % other genes. Those regons are then saved in the field "non_unique_regions"
|
|
|
4
|
|
|
5 CHROMOSOMES={};
|
|
|
6 COUNTER=1;
|
|
|
7 for i=1:size(genes,2)
|
|
|
8 CHROMOSOMES{COUNTER}=genes(i).chr;
|
|
|
9 COUNTER=COUNTER+1;
|
|
|
10 end
|
|
|
11 CHROMOSOMES=unique(CHROMOSOMES);
|
|
|
12
|
|
|
13
|
|
|
14 INFO=zeros(size(genes,2),4);
|
|
|
15 for i=1:size(genes,2)
|
|
|
16 CHR_VAL=0;
|
|
|
17 for chr= 1:length(CHROMOSOMES)
|
|
|
18 if strcmp(genes(i).chr,CHROMOSOMES(chr))
|
|
|
19 CHR_VAL=chr;
|
|
|
20 end
|
|
|
21 end
|
|
|
22 INFO(i,:)=[i,genes(i).start,genes(i).stop, CHR_VAL];
|
|
|
23 end
|
|
|
24
|
|
|
25 COUNTER=1;
|
|
|
26 new_genes=genes;
|
|
|
27 for chr= 1:length(CHROMOSOMES)
|
|
|
28 GENES_ON_CHR=INFO(INFO(:,4)==chr,:);
|
|
|
29 [TEMP,POS]=sort(GENES_ON_CHR(:,2));
|
|
|
30 GENES_ON_CHR=GENES_ON_CHR(POS,:);
|
|
|
31 STARTS=GENES_ON_CHR(:,2);
|
|
|
32 STOPS=GENES_ON_CHR(:,3);
|
|
|
33 for i=1:(size(GENES_ON_CHR,1))
|
|
|
34 MIN_START=find(STOPS>=STARTS(i),1,'first');
|
|
|
35 MAX_STOP=find(STARTS<=STOPS(i),1,'last');
|
|
|
36 if MIN_START==i
|
|
|
37 MIN_START=[];
|
|
|
38 end
|
|
|
39 if MAX_STOP==i
|
|
|
40 MAX_STOP=[];
|
|
|
41 end
|
|
|
42 EXONS=[];
|
|
|
43 if not (isempty(MIN_START))
|
|
|
44 for CURR=MIN_START:(i-1)
|
|
|
45 if(not(isempty(genes(GENES_ON_CHR(CURR,1)).transcripts)))
|
|
|
46 for tra=1:size(genes(GENES_ON_CHR(CURR,1)).transcripts,2)
|
|
|
47 if(not(isempty(genes(GENES_ON_CHR(CURR,1)).exons)))
|
|
|
48 EXONS=[EXONS;genes(GENES_ON_CHR(CURR,1)).exons{tra}];
|
|
|
49 else
|
|
|
50 EXONS=[EXONS;genes(GENES_ON_CHR(CURR,1)).start,genes(GENES_ON_CHR(CURR,1)).stop];
|
|
|
51 end
|
|
|
52 end
|
|
|
53 else
|
|
|
54 EXONS=[EXONS;genes(GENES_ON_CHR(CURR,1)).start,genes(GENES_ON_CHR(CURR,1)).stop];
|
|
|
55 end
|
|
|
56 end
|
|
|
57 end
|
|
|
58 if not (isempty(MAX_STOP))
|
|
|
59 for CURR=(i+1):MAX_STOP
|
|
|
60 if(not(isempty(genes(GENES_ON_CHR(CURR,1)).transcripts)))
|
|
|
61 for tra=1:size(genes(GENES_ON_CHR(CURR,1)).transcripts,2)
|
|
|
62 if(not(isempty(genes(GENES_ON_CHR(CURR,1)).exons)))
|
|
|
63 EXONS=[EXONS;genes(GENES_ON_CHR(CURR,1)).exons{tra}];
|
|
|
64 else
|
|
|
65 EXONS=[EXONS;genes(GENES_ON_CHR(CURR,1)).start,genes(GENES_ON_CHR(CURR,1)).stop];
|
|
|
66 end
|
|
|
67 end
|
|
|
68 else
|
|
|
69 EXONS=[EXONS;genes(GENES_ON_CHR(CURR,1)).start,genes(GENES_ON_CHR(CURR,1)).stop];
|
|
|
70 end
|
|
|
71
|
|
|
72 end
|
|
|
73 end
|
|
|
74 if not (isempty([MAX_STOP,MIN_START]))
|
|
|
75 EXONS=EXONS(EXONS(:,2)>=STARTS(i),:);
|
|
|
76 EXONS=EXONS(EXONS(:,1)<=STOPS(i),:);
|
|
|
77 new_genes(GENES_ON_CHR(i,1)).non_unique_regions=EXONS;
|
|
|
78 else
|
|
|
79 new_genes(GENES_ON_CHR(i,1)).non_unique_regions=[];
|
|
|
80 end
|
|
|
81 end
|
|
|
82 COUNTER=COUNTER+1;
|
|
|
83 end
|
