Mercurial > repos > vipints > rdiff

comparison rDiff/src/tools/detect_overlapping_regions.m @ 0:0f80a5141704

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
version 0.3 uploaded
author vipints
date Thu, 14 Feb 2013 23:38:36 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:0f80a5141704
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