annotate rDiff/src/tests/rDiff_mmd.m @ 3:29a698dc5c7e default tip

Merge multiple heads.
author Dave Bouvier <dave@bx.psu.edu>
date Mon, 27 Jan 2014 14:15:36 -0500
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 [pval,info] = rDiff_mmd(CFG,reads1,reads2)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
2 % simple application of mmd to test for differential distributions
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
3 % of reads1, reads2
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
4 % reads1: N1 x L
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
5 % reads2: N2 x L
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
6
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
7
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
8 bootstraps=CFG.bootstraps;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
9
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
10
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
11 % ensure reads are sparse and remove zero collumns
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
12 %reads1temp = sparse(reads1(:,sum([reads1;reads2],1)>0));
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
13 %reads2 = sparse(reads2(:,sum([reads1;reads2],1)>0));
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
14 %reads1=reads1temp;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
15
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
16 statistic = eucl_dist(mean(reads1,1),mean(reads2,1))^2;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
17
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
18 allreads = [reads1;reads2];
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
19
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
20 N1 = size(reads1,1);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
21 N2 = size(reads2,1);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
22 N = N1 + N2;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
23
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
24
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
25 %Use the transpose to make the selection of columms faster
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
26 all_reads_trans=allreads';
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
27
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
28 %bootstraping
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
29 for i = 1:bootstraps
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
30
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
31 r = randperm(N);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
32
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
33 sample1 = all_reads_trans(:,r(1:N1));
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
34 sample2 = all_reads_trans(:,r(N1+1:N));
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
35
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
36 bootstrap_results(i) = eucl_dist(mean(sample1,2), mean(sample2,2))^2;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
37
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
38 end
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
39
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
40 %Calculate the p-value
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
41 pval = min(1,double(1+sum(bootstrap_results >= statistic)) / bootstraps);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
42 info = [];
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
43
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
44
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
45
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
46 function result = eucl_dist(A,B)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
47
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
48 result = sqrt(sum( (A - B) .^ 2 ));