annotate rDiff/src/get_nonparametric_tests_caller.m @ 2:233c30f91d66

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