comparison deseq-hts_1.0/src/mask_dubl.m @ 0:94a108763d9e draft

deseq-hts version 1.0 wraps the DESeq 1.6.0
author vipints
date Wed, 09 May 2012 20:43:47 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:94a108763d9e
1 function [new_genes]=mask_dubl(genes,THRESH);
2
3 CHROMOSOMES={};
4 COUNTER=1;
5 for i=1:size(genes,2)
6 CHROMOSOMES{COUNTER}=genes(i).chr;
7 COUNTER=COUNTER+1;
8 end
9 CHROMOSOMES=unique(CHROMOSOMES);
10
11
12 INFO=zeros(size(genes,2),4);
13 for i=1:size(genes,2)
14 CHR_VAL=0;
15 for chr= 1:length(CHROMOSOMES)
16 if strcmp(genes(i).chr,CHROMOSOMES(chr))
17 CHR_VAL=chr;
18 end
19 end
20 INFO(i,:)=[i,genes(i).start,genes(i).stop, CHR_VAL];
21 end
22
23 COUNTER=1;
24 new_genes=genes;
25 for chr= 1:length(CHROMOSOMES)
26 GENES_ON_CHR=INFO(INFO(:,4)==chr,:);
27 [TEMP,POS]=sort(GENES_ON_CHR(:,2));
28 GENES_ON_CHR=GENES_ON_CHR(POS,:);
29 STARTS=GENES_ON_CHR(:,2);
30 STOPS=GENES_ON_CHR(:,3);
31 for i=1:(size(GENES_ON_CHR,1))
32 MIN_START=find(STOPS>=STARTS(i),1,'first');
33 MAX_STOP=find(STARTS<=STOPS(i),1,'last');
34 if MIN_START==i
35 MIN_START=[];
36 end
37 if MAX_STOP==i
38 MAX_STOP=[];
39 end
40 EXONS=[];
41 if not (isempty(MIN_START))
42 for CURR=MIN_START:(i-1)
43 if(not(isempty(genes(GENES_ON_CHR(CURR,1)).transcripts)))
44 for tra=1:size(genes(GENES_ON_CHR(CURR,1)).transcripts,2)
45 if(not(isempty(genes(GENES_ON_CHR(CURR,1)).exons)))
46 EXONS=[EXONS;genes(GENES_ON_CHR(CURR,1)).exons{tra}];
47 else
48 EXONS=[EXONS;genes(GENES_ON_CHR(CURR,1)).start,genes(GENES_ON_CHR(CURR,1)).stop];
49 end
50 end
51 else
52 EXONS=[EXONS;genes(GENES_ON_CHR(CURR,1)).start,genes(GENES_ON_CHR(CURR,1)).stop];
53 end
54 end
55 end
56 if not (isempty(MAX_STOP))
57 for CURR=(i+1):MAX_STOP
58 if(not(isempty(genes(GENES_ON_CHR(CURR,1)).transcripts)))
59 for tra=1:size(genes(GENES_ON_CHR(CURR,1)).transcripts,2)
60 if(not(isempty(genes(GENES_ON_CHR(CURR,1)).exons)))
61 EXONS=[EXONS;genes(GENES_ON_CHR(CURR,1)).exons{tra}];
62 else
63 EXONS=[EXONS;genes(GENES_ON_CHR(CURR,1)).start,genes(GENES_ON_CHR(CURR,1)).stop];
64 end
65 end
66 else
67 EXONS=[EXONS;genes(GENES_ON_CHR(CURR,1)).start,genes(GENES_ON_CHR(CURR,1)).stop];
68 end
69
70 end
71 end
72 if not (isempty([MAX_STOP,MIN_START]))
73 EXONS=EXONS(EXONS(:,2)>=STARTS(i),:);
74 EXONS=EXONS(EXONS(:,1)<=STOPS(i),:);
75 new_genes(GENES_ON_CHR(i,1)).non_unique_regions=EXONS;
76 else
77 new_genes(GENES_ON_CHR(i,1)).non_unique_regions=[];
78 end
79 end
80 COUNTER=COUNTER+1;
81 end