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