|
0
|
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
|
