Mercurial > repos > vipints > rdiff

comparison rDiff/src/get_nonparametric_tests_caller.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 []=get_nonparametric_tests_caller(PAR)
2
3
4 CFG = PAR.CFG;
5 genes = PAR.genes;
6 OUT_STR='';
7
8 % add paths
9 addpath(CFG.paths);
10 %load local variables
11 data_dir=CFG.data_dir;
12 OUT_STR=[];
13
14 variance_function_nonparametric_1=PAR.variance_function_nonparametric_1;
15 variance_function_nonparametric_2=PAR.variance_function_nonparametric_2;
16
17 Counts_rDiff_nonparametric=PAR.Counts_rDiff_nonparametric;
18 Gene_expression=PAR.Gene_expression;
19
20 %clear variabe PAR
21 clear PAR;
22
23 NUMBER_OF_TESTS_PER_GENE=(CFG.perform_nonparametric+CFG.perform_mmd);
24 if NUMBER_OF_TESTS_PER_GENE==0
25 return
26 end
27 P_VALS=cell(size(genes,2),NUMBER_OF_TESTS_PER_GENE+2);
28
29
30 %iterate over genes
31 for i=1:size(genes,2)
32
33 %TEMP_COUNT contains the counts for the current gene
34 TEMP_COUNT=cell(1,3);
35 gene = genes(i);
36
37
38 OLD_OUT_STR=OUT_STR;
39 OUT_STR=['Current gene: ' gene.name ];
40 %print progress
41 if CFG.use_rproc
42 fprintf([OUT_STR '\n'])
43 else
44 % Erase old progress
45 fprintf(repmat('\b',1,length(OLD_OUT_STR)));
46 fprintf([OUT_STR])
47 end
48
49 %set default return values
50 P_VALS{i,1}=gene.name;
51
52 %check that the gene has exons defined
53 if isempty(gene.exons)
54 P_VALS{i,4}='Exons field empty in gene structure';
55 continue;
56 end
57
58 %check that the gene is longer than the Reads. Otherwise the
59 %definition of regions does not makes sense
60 if gene.stop-gene.start<CFG.sequenced_length+3
61 continue;
62 end
63 %perform
64
65
66 %Get the reads
67 READ_SET={};
68 for bam_ix=1:length(CFG.BAM_FILES)
69 CFG.curr_bamfile=CFG.BAM_FILES{bam_ix};
70 READ_SET{end+1} = get_reads_for_gene(CFG,gene);
71 end
72 %merge the reads per sample
73 SAMPLE1=find(CFG.SAMPLES==1);
74 SAMPLE2=find(CFG.SAMPLES==2);
75
76
77 COUNTER=2;
78
79 if CFG.perform_mmd
80 reads1=[];
81 reads2=[];
82 SAMPLE_LENGTH1=length(SAMPLE1);
83 SAMPLE_LENGTH2=length(SAMPLE2);
84 for bam_ix=1:SAMPLE_LENGTH1
85 reads1=[reads1;READ_SET{SAMPLE1(bam_ix)}];
86 end
87
88 for bam_ix=1:SAMPLE_LENGTH2
89 reads2=[reads2;READ_SET{SAMPLE2(bam_ix)}];
90 end
91 [PV, INFO]= rDiff_mmd(CFG,reads1,reads2);
92 P_VALS{i,COUNTER}={PV, INFO};
93 clear reads1
94 clear reads2
95 COUNTER=COUNTER+1;
96 end
97
98 if CFG.perform_nonparametric
99 [PV, INFO]= rDiff_nonparametric(CFG,READ_SET(SAMPLE1),READ_SET(SAMPLE2),variance_function_nonparametric_1, variance_function_nonparametric_2);
100 P_VALS{i,COUNTER}={PV, INFO};
101 COUNTER=COUNTER+1;
102 end
103
104 end
105 fprintf('\n')
106 %Save the p-values
107 OUT_FILENAME=[CFG.outfile_prefix '.mat'];
108 save(OUT_FILENAME,'P_VALS')
109