Mercurial > repos > vipints > rdiff

comparison rDiff/src/tests/rDiff_parametric.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 [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