|
0
|
1 function [P_VALUE, RET_STRUCT]= rDiff_parametric(CFG,gene,Counts_rDiff_parametric,Gene_expression,variance_function_parametric_1, variance_function_parametric_2)
|
|
|
2
|
|
|
3 % Calculates the p-Values of a negative binomial test on each
|
|
|
4 % alternative regions and combines the p-values using Bonferroni's correction
|
|
|
5
|
|
|
6
|
|
|
7 %Initialize gene.name
|
|
|
8 NR_OF_TRANS=size(gene.transcripts,2);
|
|
|
9 if NR_OF_TRANS<=1
|
|
|
10 RET_STRUCT='NR_OF_TRANS too small';
|
|
|
11 P_VALUE=1;
|
|
|
12 return
|
|
|
13 end
|
|
|
14
|
|
|
15
|
|
|
16 %get the samples that are expressed (have more than 10 reads)
|
|
|
17 TEMP_SAMPLE1=and(CFG.SAMPLES==1,Gene_expression>=10);
|
|
|
18 TEMP_SAMPLE2=and(CFG.SAMPLES==2,Gene_expression>=10);
|
|
|
19 SAMPLE1=find(TEMP_SAMPLE1);
|
|
|
20 SAMPLE2=find(TEMP_SAMPLE2);
|
|
|
21
|
|
|
22 %Check wether Counts_rDiff_parametric is nonempty
|
|
|
23 for j=1:length(TEMP_SAMPLE1)
|
|
|
24 TEMP_SAMPLE1(j)=and(not(isempty(Counts_rDiff_parametric{j})),TEMP_SAMPLE1(j));
|
|
|
25 end
|
|
|
26 for j=1:length(TEMP_SAMPLE2)
|
|
|
27 TEMP_SAMPLE2(j)=and(not(isempty(Counts_rDiff_parametric{j})),TEMP_SAMPLE2(j));
|
|
|
28 end
|
|
|
29
|
|
|
30 SAMPLE1=find(TEMP_SAMPLE1);
|
|
|
31 SAMPLE2=find(TEMP_SAMPLE2);
|
|
|
32
|
|
|
33
|
|
|
34 SAMPLE_LENGTH1=length(SAMPLE1);
|
|
|
35 SAMPLE_LENGTH2=length(SAMPLE2);
|
|
|
36
|
|
|
37 if min(SAMPLE_LENGTH1,SAMPLE_LENGTH2)==0
|
|
|
38 RET_STRUCT='SAMPLE_LENGTH too small';
|
|
|
39 P_VALUE=1;
|
|
|
40 return
|
|
|
41 end
|
|
|
42
|
|
|
43 % Get the region counts
|
|
|
44 region_counts_1=zeros(SAMPLE_LENGTH1,length(Counts_rDiff_parametric{1,1}));
|
|
|
45 for j=1:SAMPLE_LENGTH1
|
|
|
46 region_counts_1(j,:)=Counts_rDiff_parametric{1,SAMPLE1(j)};
|
|
|
47 end
|
|
|
48 region_counts_2=zeros(SAMPLE_LENGTH2,length(Counts_rDiff_parametric{1,1}));
|
|
|
49 for j=1:SAMPLE_LENGTH2
|
|
|
50 region_counts_2(j,:)=Counts_rDiff_parametric{1,SAMPLE2(j)};
|
|
|
51 end
|
|
|
52
|
|
|
53 % Get the gene expression
|
|
|
54 gene_expression_1=Gene_expression(SAMPLE1);
|
|
|
55 gene_expression_2=Gene_expression(SAMPLE2);
|
|
|
56
|
|
|
57 % compute the expected mean and the variance under the null
|
|
|
58 % hypothesis
|
|
|
59 [EXPECTED_MEAN,EXPECTED_VARIANCE]=get_mean_variance_seg(gene_expression_1,gene_expression_2,region_counts_1,region_counts_2,variance_function_parametric_1, variance_function_parametric_2);
|
|
|
60
|
|
|
61 %compute the p-values
|
|
|
62 OBSERVED_COUNTS=round([mean(region_counts_1,1);mean(region_counts_2,1)]);
|
|
|
63
|
|
|
64 P_LIST=ones(1,size(EXPECTED_MEAN,2));
|
|
|
65 %Iterate over the regions
|
|
|
66 SKIPPED_TESTS=0;
|
|
|
67 for i=1:length(P_LIST)
|
|
|
68 if sum(OBSERVED_COUNTS(:,i))==0
|
|
|
69 SKIPPED_TESTS=SKIPPED_TESTS+1;
|
|
|
70 continue
|
|
|
71 end
|
|
|
72
|
|
|
73 [P_VALUE,FL]=comp_nbin_p_value_mean_variance(EXPECTED_MEAN(1,i),EXPECTED_MEAN(2,i),EXPECTED_VARIANCE(1,i),EXPECTED_VARIANCE(2,i),OBSERVED_COUNTS(1,i),OBSERVED_COUNTS(2,i));
|
|
|
74 P_LIST(i)=P_VALUE;
|
|
|
75 end
|
|
|
76 if length(P_LIST)-SKIPPED_TESTS<=0
|
|
|
77 P_VALUE=1;
|
|
|
78 else
|
|
|
79 P_VALUE=min(P_LIST)*(length(P_LIST)-SKIPPED_TESTS);
|
|
|
80 end
|
|
|
81 RET_STRUCT={};
|
|
|
82 return |
