+ removed several index bugs
[qpalma.git] / tools / calculateAlignmentQuality.m
1 %
2 % This script compares two alignments based on their EST overlaps
3 %
4 % The ground_truth and testrun
5 %
6 load /fml/ag-raetsch/home/fabio/tmp/zebrafish/confirmed_sequences.mat;
7 testrun = genes;
8 load /fml/ag-raetsch/share/projects/altsplicedata/zebrafish/confirmed_sequences.mat;
9 ground_truth = genes;
10 clear genes;
11
12 fh = fopen('overlapping.pos','w+')
13 for i = 1:length(ground_truth)
14 currentEntry = ground_truth(i);
15 currentExons = currentEntry.exons;
16 assert (length(currentEntry.transcripts) == length(currentEntry.exons));
17 numberOfEsts = length(currentEntry.transcripts);
18
19 for j = 1:length(testrun)
20 currentPred = testrun(j);
21
22 if strcmp(currentEntry.chr,currentPred.chr) && currentEntry.is_alt_spliced
23 continue;
24 end
25
26 start = currentEntry.start;
27 stop = currentEntry.stop;
28 predStart = currentPred.start;
29 predStop = currentPred.stop;
30
31 % if both entries are not overlapping at all
32 if predStop < start || predStart > stop
33 continue;
34 end
35
36 currentPredExons = currentPred.exons;
37 numberOfPredEsts = length(currentPredExons);
38
39 % Loop over all ESTs in ground truth and testrun
40 for estIdx = 1:numberOfEsts
41 %disp(sprintf('estIdx %i ',estIdx));
42 numberOfExons = size(currentExons{estIdx},1);
43
44 for exonIdx = 1:numberOfExons
45 %disp(sprintf('exonIdx %i ',exonIdx));
46 exonStart = currentExons{estIdx}(exonIdx,1);
47 exonStop = currentExons{estIdx}(exonIdx,2);
48
49 for predEstIdx = 1:numberOfPredEsts
50 %disp(sprintf('predEstIdx %i ',predEstIdx));
51 numberOfPredExons = size(currentPredExons{predEstIdx},1);
52 currentESTName = currentPred.transcripts{predEstIdx};
53
54 for predExonIdx = 1:numberOfPredExons
55 %%disp(sprintf('predExonIdx %i ',predExonIdx));
56 predExonStart = currentPredExons{predEstIdx}(predExonIdx,1);
57 predExonStop = currentPredExons{predEstIdx}(predExonIdx,2);
58
59 %disp('\n');
60 %%disp(sprintf('%i %i %i %i %i %i\n',i,j,estIdx,predEstIdx,exonIdx,predExonIdx));
61
62 if exonStart >= predExonStart && exonStart <= predExonStop
63 %%disp('Overlapping');
64 fprintf(fh,sprintf('%s before %d\n',currentESTName,exonStart-predExonStart))
65 % end of exonth is nested inside predExoniction
66 elseif exonStop >= predExonStart && exonStop <= predExonStop
67 %%disp('Overlapping');
68 fprintf(fh,sprintf('%s after %d\n',currentESTName,predExonStop-exonStop))
69 % predExoniction is nested inside exonth
70 elseif exonStart <= predExonStart && predExonStop <= exonStop
71 %disp('Overlapping');
72 fprintf('%s %d %d',fh)
73 else
74 d=1;
75 end
76 end
77 end
78 end
79 end
80 end
81 end