-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy path02_EF.R
65 lines (58 loc) · 2.13 KB
/
02_EF.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
rm(list=ls());
source('supp_func.R');
pdf("2_EF.pdf", 20,10);
layout(matrix(c(1,2), nrow = 1, ncol = 2, byrow = TRUE))
par("mar"=c(7, 7, 7, 2))
clu.n <- c(seq(40, 100, length.out = 7), seq(200, 1000, length.out = 9));
cols <- rainbow(10, alpha=0.8);
dat0 <- read.table("../data/02/cor_glm_clu_CD4_n5/40.txt", header=T);
gen0 <- paste0(dat0[1:10,1], "_", dat0[1:10,2]);
out.s <- NULL;
for (n in clu.n){
inff <- paste0("../data/02/cor_glm_clu_CD4_n5/", n, ".txt");
dat2 <- read.table(inff, header=T);
dat2[,8] <- -log10(dat2[,8]);
gen2 <- paste0(dat2[,1], "_", dat2[,2]);
g2.i <- match(gen0, gen2);
df.s <- data.frame(clu=rep(n, length(gen0)), genes=gen0,
pearson.R=dat2[g2.i,3], spearman.R=dat2[g2.i,4], fdr=dat2[g2.i,8]);
out.s <- rbind(out.s, df.s);
}
gen.u <- unique(out.s[,2]);
i <- which(is.infinite(out.s[,5]));
out.s[i, 5] <- max(out.s[-i,5]);
ylim <- range(out.s[,5]);
i <- 0;
for (g in gen.u){
i <- i + 1;
if (i == 1){
g.j <- which(out.s[,2]==g);
plot(out.s[g.j,5], type="b", col=cols[i], xlab="", ylab="",
ylim=ylim, pch=19, lwd=4,cex.main=3,cex.lab=2.5, cex.axis=2.5,xaxt='n',cex = 4);
title(ylab="-log10 P Value", line=4, cex.lab =4.5, family = "sans")
title(xlab="Number of Clusters", line=5, cex.lab =4.5, family = "sans")
} else {
g.j <- which(out.s[,2]==g);
points(out.s[g.j,5], type="b", col=cols[i], pch=19, lwd=4,cex = 4);
}
}
axis(1, at=1:length(clu.n), labels=clu.n,cex.axis=2.5)
fig_label('E', pos='topleft',cex=5)
ylim <- range(out.s[,3]);
i <- 0;
for (g in gen.u){
i <- i + 1;
if (i == 1){
g.j <- which(out.s[,2]==g);
plot(out.s[g.j,3], type="b", col=cols[i], xlab="", ylab="",
ylim=ylim, pch=19, lwd=4,cex.main=3,cex.lab=2, cex.axis=2.5,cex = 4,xaxt='n');
title(ylab="Pearson R Value", line=4, cex.lab =4.5, family = "sans")
title(xlab="Number of Clusters", line=5, cex.lab =4.5, family = "sans")
} else {
g.j <- which(out.s[,2]==g);
points(out.s[g.j,3], type="b", col=cols[i], pch=19, lwd=4,cex = 4);
}
}
axis(1, at=1:length(clu.n), labels=clu.n,cex.axis=2.5)
fig_label('F', pos='topleft',cex=5)
dev.off()