annotate rDiff/src/rdiff_caller.m @ 3:29a698dc5c7e default tip

Merge multiple heads.
author Dave Bouvier <dave@bx.psu.edu>
date Mon, 27 Jan 2014 14:15:36 -0500
parents 0f80a5141704
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
1 function STAT = rDiff_caller(PAR)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
2 % genes = opt_transcripts_caller(PAR)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
3 %
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
4 % -- input --
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
5 % PAR contains
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
6 % CFG: configuration struct
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
7 % genes: struct defining genes with start, stops, exons etc.
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
8 % profile_weights: weights of profile functions
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
9 % intron_dists: distances to closest intron
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
10 %
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
11 % -- output --
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
12 % genes: struct with additional fields of eg. estimated transcript weights
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
13 addpath ../../tests
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
14 addpath ../../variance/final
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
15 addpath ../../experimental
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
16
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
17 samtools='/fml/ag-raetsch/home/drewe/svn/projects/RNASeq/difftest/samtools/';
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
18 CFG = PAR.CFG;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
19 SUFF_NAME=CFG.SUFF_NAME;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
20 %load(CFG.genes_path , 'genes');
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
21 %[genes]=mask_dubl(genes,0);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
22
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
23 load(CFG.genes_path , 'genes');
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
24 [genes]=mask_dubl(genes,0);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
25
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
26 %[genes]=remove_transcripts(genes);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
27 genes=genes(PAR.genes_idx);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
28 BAM_FILES=CFG.BAM_FILES;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
29 JB_NR=CFG.JB_NR;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
30
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
31
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
32
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
33 %clear PAR;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
34
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
35 %%%% paths
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
36 addpath(CFG.paths);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
37 addpath('/fml/ag-raetsch/home/drewe/svn/projects/RNASeq/difftest/variance/nbin');
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
38 addpath('/fml/ag-raetsch/home/drewe/svn/projects/RNASeq/difftest/experiments/');
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
39 addpath('/fml/ag-raetsch/home/drewe/svn/projects/RNASeq/difftest/tests/');
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
40 addpath('/fml/ag-raetsch/home/drewe/svn/projects/RNASeq/difftest/read_utils/');
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
41 addpath('/fml/ag-raetsch/home/drewe/svn/projects/RNASeq/difftest/tests/sequence_tools/');
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
42 addpath(genpath('/fml/ag-raetsch/home/drewe/svn/tools/chronux/'))
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
43 addpath('/fml/ag-raetsch/home/drewe/svn/projects/RNASeq/difftest/tests/functions/');
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
44
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
45 out_base=CFG.out_base;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
46 %sim_data_base=CFG.sim_data_base;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
47 real_data_base=CFG.real_data_base;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
48
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
49 expr12_bam=CFG.expr12_bam;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
50 expr11_bam=CFG.expr11_bam;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
51 expr22_bam=CFG.expr22_bam;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
52 expr21_bam=CFG.expr21_bam;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
53
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
54 CFG.SEQUENCED_LENGTH=80;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
55 SEQUENCED_LENGTH=80;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
56 %make map: gene->file
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
57
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
58 IDX1=CFG.IDX11;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
59 IDX2=CFG.IDX22;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
60
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
61 %BAM_FILES{IDX1}
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
62 %BAM_FILES{IDX2}
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
63
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
64
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
65 %genes=PAR.genes;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
66 clear PAR
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
67
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
68 STAT=cell(size(genes));
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
69 do_save = 0;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
70 out_base = '/fml/ag-raetsch/home/drewe/svn/projects/RNASeq/difftest/out';
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
71 save_dir=fullfile(out_base,'cache');
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
72
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
73 for i=1:size(genes,2)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
74 %try
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
75 RESULT=cell(1,7)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
76 gene = genes(i);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
77
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
78 RESULT{7}=JB_NR;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
79 RESULT{1}=gene.name
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
80 load_only = false;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
81
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
82 if or(isempty(gene.exons),gene.stop-gene.start<=SEQUENCED_LENGTH)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
83 PV=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
84 RESULT{2}={PV,''};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
85 RESULT{3}={PV,''};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
86 RESULT{4}={PV,''};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
87 RESULT{5}=[Inf,Inf];
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
88 RESULT{6}=[Inf,Inf];
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
89 RESULT{7}={'empty gene exons'};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
90 RESULT{8}={PV,''};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
91 continue;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
92 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
93
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
94 PV1=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
95 for l=1:size(gene.transcripts,2)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
96 EXONS=gene.exons{l}
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
97 if sum(EXONS(:,2)-EXONS(:,1))<=SEQUENCED_LENGTH
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
98 PV=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
99 RESULT{2}={PV,''};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
100 RESULT{3}={PV,''};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
101 RESULT{4}={PV,''};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
102 RESULT{5}=[Inf,Inf];
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
103 RESULT{6}=[Inf,Inf];
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
104 RESULT{7}={'empty gene exons'};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
105 RESULT{8}={PV,''};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
106 PV1=2;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
107
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
108 continue;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
109 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
110 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
111 if PV1==2
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
112 continue;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
113 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
114 gene.name
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
115
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
116
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
117
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
118 %adjust start and stop positions on the negative strand
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
119 if strcmp(gene.strand,'-')
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
120 gene.start=gene.start+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
121 gene.stop=gene.stop+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
122 [reads11,reads12,intron11,intron12] = get_reads_dual2_intron(expr11_bam,expr12_bam,gene,samtools,false,10);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
123 [reads11,FLAG]=remove_reads_from_other_genes(reads11,gene);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
124 [reads12,FLAG]=remove_reads_from_other_genes(reads12,gene);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
125
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
126
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
127 [reads21,reads22,intron21,intron22] = get_reads_dual2_intron(expr21_bam,expr22_bam,gene,samtools,false,10);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
128 [reads21,FLAG]=remove_reads_from_other_genes(reads21,gene);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
129 [reads22,FLAG]=remove_reads_from_other_genes(reads22,gene);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
130 gene.start=gene.start-1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
131 gene.stop=gene.stop-1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
132 else
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
133 [reads11,reads12,intron11,intron12] = get_reads_dual2_intron(expr11_bam,expr12_bam,gene,samtools,false,10);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
134 [reads11,FLAG]=remove_reads_from_other_genes(reads11,gene);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
135 [reads12,FLAG]=remove_reads_from_other_genes(reads12,gene);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
136
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
137 [reads21,reads22,intron21,intron22] = get_reads_dual2_intron(expr21_bam,expr22_bam,gene,samtools,false,10);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
138 [reads21,FLAG]=remove_reads_from_other_genes(reads21,gene);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
139 [reads22,FLAG]=remove_reads_from_other_genes(reads22,gene);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
140 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
141
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
142 reads11=reads11(sum(reads11,2)>70,:);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
143 reads12=reads12(sum(reads12,2)>70,:);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
144 reads21=reads21(sum(reads21,2)>70,:);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
145 reads22=reads22(sum(reads22,2)>70,:);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
146 %% TEST with read cliping
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
147
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
148 READS1={reads11,reads12};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
149 READS2={reads21,reads22};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
150 reads1=[reads11;reads12];
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
151 reads2=[reads21;reads22];
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
152
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
153 TOTAL_SIZE=(size(reads12,1)+size(reads12,1)+(size(reads21,1)+size(reads22,1)));
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
154 MIN_SIZE=min(size(reads12,1)+size(reads12,1),(size(reads21,1)+size(reads22,1)));
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
155 %keyboard
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
156
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
157 CLIP=0;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
158 if 1==0
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
159 if (size(reads1,1)>0)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
160 for i=1:size(reads1,1)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
161 reads1(i,find(reads1(i,:),2,'first'))=0;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
162 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
163 for i=1:size(reads1,1)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
164 reads1(i,find(reads1(i,:),2,'last'))=0;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
165 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
166 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
167 if (size(reads2,1)>0)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
168 for i=1:size(reads2,1)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
169 reads2(i,find(reads2(i,:),2,'first'))=0;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
170 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
171 for i=1:size(reads2,1)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
172 reads2(i,find(reads2(i,:),2,'last'))=0;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
173 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
174 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
175 CLIP=4;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
176 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
177
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
178 % keybloard
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
179 %% END TEST
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
180
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
181 COUNTER=3;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
182 FLAG=0;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
183 for k=1:size(gene.transcripts,2)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
184 if sum(gene.exons{k}(:,2)-gene.exons{k}(:,1))<75
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
185 FLAG=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
186 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
187 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
188
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
189 if FLAG==1
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
190 RESULT{2}={'gene shorter than 75 bp'};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
191 STAT{i}=RESULT;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
192 continue
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
193 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
194
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
195 if 1==0
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
196
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
197 for V_COUNTER=2%:length(CFG.VARIANCES1)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
198 if (TOTAL_SIZE>0)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
199 if(MIN_SIZE>0)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
200 PV=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
201 INFO='';
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
202 [PV,INFO] =diff_mmd([reads11;reads12],[reads21;reads22],gene);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
203 PV
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
204 RESULT{COUNTER}={PV,INFO};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
205 COUNTER=COUNTER+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
206 else
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
207 PV=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
208 RESULT{COUNTER}={PV,2};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
209 COUNTER=COUNTER+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
210 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
211 else
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
212 PV=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
213 RESULT{COUNTER}={PV,1};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
214 COUNTER=COUNTER+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
215 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
216 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
217 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
218 % keyboard
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
219
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
220 VAR_FACT=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
221 COV_POS=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
222 REWEIGHT=0;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
223 if 1==1
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
224 for V_COUNTER=2%:length(CFG.VARIANCES1)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
225 if (TOTAL_SIZE>0)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
226 if(MIN_SIZE>0)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
227 cen_ARR=0.1:0.1:1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
228 REWEIGHT_WEIGHTS=zeros(length( cen_ARR),size(reads1,2));
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
229 cen_count=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
230 for censor_frac= cen_ARR
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
231 %if strcmp(gene.name,'AT1G01073')
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
232 % keyboard
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
233 %end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
234 censor_frac
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
235 temp_reads1=reads1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
236 temp_reads2=reads2;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
237 %cut to relevant positions
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
238 COVERAGE=sum(temp_reads1,1)+sum(temp_reads2,1);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
239 NONZERO=COVERAGE>0;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
240 SORTED_COVERAGES=sort(COVERAGE(NONZERO));
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
241 NR_OF_NONZERO=sum(NONZERO);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
242 CHOSEN_POSITIONS=COVERAGE<=SORTED_COVERAGES(ceil(NR_OF_NONZERO*censor_frac));
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
243 REWEIGHT_WEIGHTS(cen_count,CHOSEN_POSITIONS)=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
244 cen_count=cen_count+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
245 %temp_reads1=temp_reads1(:,CHOSEN_POSITIONS);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
246 %temp_reads2=temp_reads2(:,CHOSEN_POSITIONS);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
247 %remove reads which have no coverage;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
248 %temp_reads1=temp_reads1(sum(temp_reads1,2)>0,:);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
249 %temp_reads2=temp_reads2(sum(temp_reads2,2)>0,:);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
250 %SIZE1=size(temp_reads1,1);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
251 %SIZE2=size(temp_reads2,1);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
252 %MIN_SIZE=min(SIZE1,SIZE2);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
253 %TOTAL_SIZE=SIZE1+SIZE2;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
254 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
255 if (TOTAL_SIZE>0)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
256 if(MIN_SIZE>0)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
257 PV=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
258 INFO='';
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
259 REWEIGHT=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
260 [PV,INFO] =diff_mmd_variance_subsample3({reads11,reads12},{reads21,reads22},gene,CFG.VARIANCES1{V_COUNTER},CFG.VARIANCES2{V_COUNTER},1,1,REWEIGHT, REWEIGHT_WEIGHTS)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
261 RESULT{COUNTER}={PV,INFO};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
262 COUNTER=COUNTER+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
263 else
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
264 PV=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
265 RESULT{COUNTER}={PV,2};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
266 COUNTER=COUNTER+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
267 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
268 else
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
269 PV=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
270 RESULT{COUNTER}={PV,1};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
271 COUNTER=COUNTER+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
272 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
273 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
274 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
275 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
276 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
277
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
278 if 1==0
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
279 VAR_FACT=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
280 COV_POS=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
281 REWEIGHT=0;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
282 for V_COUNTER=2%:length(CFG.VARIANCES1)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
283 for cov_fact_ix=1:length(COV_POS)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
284
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
285 if (TOTAL_SIZE>0)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
286 if(MIN_SIZE>0)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
287 PV=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
288 INFO='';
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
289 [PV,INFO] =diff_mmd_variance_subsample2({reads11,reads12},{reads21,reads22},gene,CFG.VARIANCES1{V_COUNTER},CFG.VARIANCES2{V_COUNTER},VAR_FACT(var_fact_ix),COV_POS,REWEIGHT);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
290 RESULT{COUNTER}={PV,INFO};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
291 COUNTER=COUNTER+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
292 else
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
293 PV=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
294 RESULT{COUNTER}={PV,2};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
295 COUNTER=COUNTER+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
296 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
297 else
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
298 PV=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
299 RESULT{COUNTER}={PV,1};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
300 COUNTER=COUNTER+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
301 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
302 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
303 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
304
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
305 if 1==1
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
306 VAR_FACT=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
307 COV_POS=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
308 REWEIGHT=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
309 for V_COUNTER=2%:length(CFG.VARIANCES1)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
310 for var_fact_ix=1:length(VAR_FACT)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
311
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
312 if (TOTAL_SIZE>0)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
313 if(MIN_SIZE>0)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
314 PV=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
315 INFO='';
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
316 [PV,INFO] =diff_mmd_variance_subsample2({reads11,reads12},{reads21,reads22},gene,CFG.VARIANCES1{V_COUNTER},CFG.VARIANCES2{V_COUNTER},VAR_FACT(var_fact_ix),COV_POS,REWEIGHT)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
317 RESULT{COUNTER}={PV,INFO};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
318 COUNTER=COUNTER+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
319 else
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
320 PV=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
321 RESULT{COUNTER}={PV,2};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
322 COUNTER=COUNTER+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
323 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
324 else
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
325 PV=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
326 RESULT{COUNTER}={PV,1};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
327 COUNTER=COUNTER+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
328 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
329 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
330 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
331 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
332 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
333
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
334 %keyboard
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
335 TOTAL_SIZE=(size(reads12,1)+size(reads12,1)+(size(reads21,1)+size(reads22,1)));
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
336 MIN_SIZE=min(size(reads12,1)+size(reads12,1),(size(reads21,1)+size(reads22,1)));
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
337 if 1==0
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
338 for V_COUNTER=1
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
339 if (TOTAL_SIZE>0)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
340 if(MIN_SIZE>0)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
341 PV=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
342 INFO='';
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
343 if isempty([intron11,intron12])|isempty([intron21,intron22])
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
344 [PV,INFO] =diff_mmd_variance_NB_NB_simple([reads11;reads12],[reads21;reads22],gene,CFG.VARIANCES1{V_COUNTER},CFG.VARIANCES2{V_COUNTER});
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
345 else
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
346 [PV,INFO] =diff_mmd_variance_splice([reads11;reads12],[reads21;reads22],0.5,[intron11,intron12],[intron21,intron22],gene,CFG.VARIANCES1{V_COUNTER},CFG.VARIANCES2{V_COUNTER});
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
347 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
348 RESULT{COUNTER}={PV,INFO};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
349 COUNTER=COUNTER+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
350 else
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
351 PV=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
352 RESULT{COUNTER}={PV,2};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
353 COUNTER=COUNTER+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
354 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
355 else
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
356 PV=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
357 RESULT{COUNTER}={PV,1};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
358 COUNTER=COUNTER+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
359 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
360 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
361 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
362 PV=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
363 %keyboard
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
364 for V_COUNTER=length(CFG.VARIANCES1)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
365 if (TOTAL_SIZE>0)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
366
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
367 if(MIN_SIZE>0)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
368 [P_VALUE, INFO]= diff_nbin7(READS1,READS2,gene,SEQUENCED_LENGTH,CFG.VARIANCES1{V_COUNTER},CFG.VARIANCES2{V_COUNTER});
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
369 RESULT{COUNTER}={P_VALUE,INFO};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
370 if not(isempty(INFO))
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
371 if iscell(INFO)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
372 RESULT{COUNTER}=INFO{5};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
373 COUNTER=COUNTER+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
374 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
375 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
376 else
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
377 PV=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
378 RESULT{COUNTER}={PV,2};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
379 COUNTER=COUNTER+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
380 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
381 else
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
382 PV=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
383 RESULT{COUNTER}={PV,1};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
384 COUNTER=COUNTER+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
385 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
386 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
387
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
388
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
389 [SPLICINGEVENTS,SEQUENCE,EXONSEQUENCE]=splicingsequence(gene);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
390 [UNIQUE_NEW_EXONS,GRAPHNODES,ORDER_OF_GRAPHNODE,EIRS_IN_SEQ]=transform_single_end_reads(SPLICINGEVENTS,SEQUENCE,EXONSEQUENCE,CFG.SEQUENCED_LENGTH-CLIP);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
391 [NEW_READS1,UNEXPLAINED_READS1,UNEXPLAINED_INDEX1]= convert_reads_to_region_indicators([reads11;reads12],UNIQUE_NEW_EXONS,GRAPHNODES,ORDER_OF_GRAPHNODE,EIRS_IN_SEQ,gene);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
392 [NEW_READS2,UNEXPLAINED_READS2,UNEXPLAINED_INDEX2]= convert_reads_to_region_indicators([reads21;reads22],UNIQUE_NEW_EXONS,GRAPHNODES,ORDER_OF_GRAPHNODE,EIRS_IN_SEQ,gene);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
393 % keyboard
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
394 TOTAL_SIZE=(size(NEW_READS1,1)+(size(NEW_READS2,1)));
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
395 MIN_SIZE=min(size(NEW_READS1,1),(size(NEW_READS2,1)));
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
396
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
397 if (TOTAL_SIZE>0)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
398 if(MIN_SIZE>0)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
399 [PV,INFO] = diff_poisson_bonf_3_unequal_segment(NEW_READS1,NEW_READS2,gene,CFG.SEQUENCED_LENGTH-CLIP);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
400 RESULT{COUNTER}={PV,INFO};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
401 COUNTER=COUNTER+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
402 else
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
403 PV=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
404 RESULT{COUNTER}={PV,2};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
405 COUNTER=COUNTER+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
406 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
407 else
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
408 PV=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
409 RESULT{COUNTER}={PV,1};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
410 COUNTER=COUNTER+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
411 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
412 PV=1
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
413 if 1==0
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
414 if (TOTAL_SIZE>0)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
415 if(MIN_SIZE>0)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
416 [PV,INFO] = diff_poisson_bonf_4_unequal_segment(NEW_READS1,NEW_READS2,gene,CFG.SEQUENCED_LENGTH-CLIP);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
417 RESULT{COUNTER}={PV,INFO};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
418 COUNTER=COUNTER+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
419 else
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
420 PV=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
421 RESULT{COUNTER}={PV,2};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
422 COUNTER=COUNTER+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
423 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
424 else
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
425 PV=1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
426 RESULT{COUNTER}={PV,1};
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
427 COUNTER=COUNTER+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
428 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
429 %keyboard
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
430 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
431 STAT{i}=RESULT;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
432
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
433 end;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
434 %keyboard
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
435 OUT_FILENAME=['/fml/ag-raetsch/home/drewe/svn/projects/RNASeq/difftest/out/analysis_artificial_variance_03_08_2012/' SUFF_NAME '_rep_mmd_07_07_2012_' int2str(JB_NR) '.mat'];
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
436 save(OUT_FILENAME,'STAT')