-
Notifications
You must be signed in to change notification settings - Fork 0
/
wilcox_test_spearman.R
125 lines (106 loc) · 5.05 KB
/
wilcox_test_spearman.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
##
## Applies a Wilcoxon rank sum of the gene-pair correlations between diseases genes and non-diseases genes for every cell-type.
##
## Loaded data:
## data/33k_correlation_spearman_celltype (created @ corr_spearman.R)
## data/immune_disease_genes (created @ immune_diseases_genes.R)
##
load("data/33k_correlation_spearman_celltype")
load("data/immune_disease_genes")
## Helper function to execute the test between disease and non-disease gene pairs.
wilcox.diseases.genes <- function(pairs, disease.genes) {
# Get the present disease and non disease genes pairs
disease.gene.pairs <- pairs[pairs$gene1 %in% disease.genes & pairs$gene2 %in% disease.genes,]
non.disease.gene.pairs <- pairs[!(pairs$gene1 %in% disease.genes & pairs$gene2 %in% disease.genes),]
# Do the test
test <- wilcox.test(disease.gene.pairs$rho^2, non.disease.gene.pairs$rho^2)
# Add the present disease genes to wilcoxon.test object
test$disease.genes <- unique(c(disease.gene.pairs$gene1, disease.gene.pairs$gene2))
test$disease.gene.pairs <- disease.gene.pairs
test$non.disease.gene.pairs <- non.disease.gene.pairs
return(test)
}
compare.distribution <- function(pairs, disease.genes) {
disease.gene.pairs <- pairs[pairs$gene1 %in% disease.genes & pairs$gene2 %in% disease.genes,]
non.disease.gene.pairs <- pairs[!(pairs$gene1 %in% disease.genes & pairs$gene2 %in% disease.genes),]
plot(density(non.disease.gene.pairs$rho^2))
#lines(density(disease.gene.pairs$rho^2))
}
## Do the actual test
# Diabetes 1 genes
nk.dia.wilcox <- wilcox.diseases.genes(pairs = nk.cells.cor.pairs, disease.genes = diabetes.1.genes)
b.dia.wilcox <- wilcox.diseases.genes(pairs = b.cells.cor.pairs, disease.genes = diabetes.1.genes)
th.dia.wilcox <- wilcox.diseases.genes(pairs = t.h.cells.cor.pairs, disease.genes = diabetes.1.genes)
tc.dia.wilcox <- wilcox.diseases.genes(pairs = t.c.cells.cor.pairs, disease.genes = diabetes.1.genes)
mono.dia.wilcox <- wilcox.diseases.genes(pairs = monocytes.cor.pairs, disease.genes = diabetes.1.genes)
# Celiac disease
nk.celiac.wilcox <- wilcox.diseases.genes(pairs = nk.cells.cor.pairs, disease.genes = celiac.genes)
b.celiac.wilcox <- wilcox.diseases.genes(pairs = b.cells.cor.pairs, disease.genes = celiac.genes)
th.celiac.wilcox <- wilcox.diseases.genes(pairs = t.h.cells.cor.pairs, disease.genes = celiac.genes)
tc.celiac.wilcox <- wilcox.diseases.genes(pairs = t.c.cells.cor.pairs, disease.genes = celiac.genes)
mono.celiac.wilcox <- wilcox.diseases.genes(pairs = monocytes.cor.pairs, disease.genes = celiac.genes)
# Rheuma
nk.rheuma.wilcox <- wilcox.diseases.genes(pairs = nk.cells.cor.pairs, disease.genes = rheuma.genes)
b.rheuma.wilcox <- wilcox.diseases.genes(pairs = b.cells.cor.pairs, disease.genes = rheuma.genes)
th.rheuma.wilcox <- wilcox.diseases.genes(pairs = t.h.cells.cor.pairs, disease.genes = rheuma.genes)
tc.rheuma.wilcox <- wilcox.diseases.genes(pairs = t.c.cells.cor.pairs, disease.genes = rheuma.genes)
mono.rheuma.wilcox <- wilcox.diseases.genes(pairs = monocytes.cor.pairs, disease.genes = rheuma.genes)
# IBD
nk.ibd.wilcox <- wilcox.diseases.genes(pairs = nk.cells.cor.pairs, disease.genes = ibd.genes)
b.ibd.wilcox <- wilcox.diseases.genes(pairs = b.cells.cor.pairs, disease.genes = ibd.genes)
th.ibd.wilcox <- wilcox.diseases.genes(pairs = t.h.cells.cor.pairs, disease.genes = ibd.genes)
tc.ibd.wilcox <- wilcox.diseases.genes(pairs = t.c.cells.cor.pairs, disease.genes = ibd.genes)
mono.ibd.wilcox <- wilcox.diseases.genes(pairs = monocytes.cor.pairs, disease.genes = ibd.genes)
## Save the results
save(list = ls(pattern=".+wilcox"), file = "data/spearman_wilcox_test_frac_0.05.RData")
load("data/spearman_wilcox_test_frac_0.05.RData")
## Explore the results
nk.dia.wilcox$p.value
b.dia.wilcox$p.value
b.dia.wilcox$disease.gene.pairs
th.dia.wilcox$p.value
tc.dia.wilcox$p.value
mono.dia.wilcox$p.value
mono.dia.wilcox$disease.genes
mono.dia.wilcox$disease.gene.pairs
nk.dia.wilcox$disease.gene.pairs
b.dia.wilcox$disease.gene.pairs
nk.celiac.wilcox$p.value
b.celiac.wilcox$p.value
th.celiac.wilcox$p.value
tc.celiac.wilcox$p.value
tc.celiac.wilcox$disease.gene.pairs
mono.celiac.wilcox$p.value
mono.celiac.wilcox$disease.gene.pairs
mono.celiac.wilcox$disease.genes
nk.rheuma.wilcox$p.value
b.rheuma.wilcox$p.value
th.rheuma.wilcox$p.value
tc.rheuma.wilcox$p.value
mono.rheuma.wilcox$p.value
nk.rheuma.wilcox$disease.genes
nk.rheuma.wilcox$disease.gene.pairs
nk.ibd.wilcox$p.value
b.ibd.wilcox$p.value
th.ibd.wilcox$p.value
tc.ibd.wilcox$p.value
mono.ibd.wilcox$p.value
nk.ibd.wilcox$disease.genes
nk.ibd.wilcox$disease.gene.pairs
mono.ibd.wilcox$disease.gene.pairs[1:10,]
mono.dia.wilcox$disease.gene.pairs[1:10,]
#### WORKSPACE
b.cells.cor.pairs
sqrt(7263766*2)
3812*3811/2
compare.distribution(b.cells.cor.pairs, diabetes.1.genes)
h = hist(nk.cells.cor.pairs$rho)
h$density = h$counts/sum(h$counts)*100
plot(h,freq=FALSE)
# Add a Normal Curve (Thanks to Peter Dalgaard)
x <- nk.cells.cor.pairs$rho
h<-hist(x, breaks=1000)
xfit<-seq(min(x),max(x),length=40)
yfit<-dnorm(xfit,mean=mean(x),sd=sd(x))
yfit <- yfit*diff(h$mids[1:2])*length(x)
lines(xfit, yfit, col="blue", lwd=2)