-
Notifications
You must be signed in to change notification settings - Fork 20
/
GSA.listsets.revised.R
98 lines (88 loc) · 2.99 KB
/
GSA.listsets.revised.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
GSA.listsets.revised =
function (GSA.obj, geneset.names = NULL, maxchar = 20, FDRcut = 0.2)
{
if(is.null(geneset.names)){
geneset.names=rep("xxxxxx", length(GSA.obj$GSA.scores))
}
negflag= !(GSA.obj$resp.type=="Multiclass" | GSA.obj$resp.type=="taCorr")
r = GSA.obj$GSA.scores
rstar = GSA.obj$GSA.scores.perm
r[is.na(r)] = 0
rstar[is.na(rstar)] = 0
nperms = ncol(GSA.obj$GSA.scores.perm)
np = length(r)
geneset.names = substring(geneset.names, 1, maxchar)
pvalues.lo = GSA.obj$pvalues.lo
pvalues.hi = GSA.obj$pvalues.hi
m = sum(!is.na(pvalues.hi))
make.monotone.increasing = function(x) {
n = length(x)
if (n==0) return(NULL) ## added to prevent error when x is NULL
if (n==1) x[1] <- x[1] ## added to prevent error when length of x =1
else {
for (i in n:2) {
if (x[i - 1] > x[i]) {
x[i - 1] = x[i]
}
}
}
return(x)
}
oo = (1:length(r))[!is.na(pvalues.hi)]
fdr.lo = fdr.hi = rep(NA, length(r))
for (i in oo) {
if(negflag){ fdr.lo[i] = round(m * pvalues.lo[i]/sum(pvalues.lo[!is.na(pvalues.lo)] <= pvalues.lo[i]), 4) }
fdr.hi[i] = round(m * pvalues.hi[i]/sum(pvalues.hi[!is.na(pvalues.hi)] <= pvalues.hi[i]), 4)
}
fdr.lo=pmin(fdr.lo,1)
fdr.hi=pmin(fdr.hi,1)
res1=NULL
if(negflag){
oo1 = (1:length(r))[r < 0]
res1 = NULL
for (i in oo1) {
res1 = rbind(res1, c(i, geneset.names[i], round(GSA.obj$GSA.scores[i],
4), pvalues.lo[i], fdr.lo[i]))
}
if (!is.null(res1)) o1 = order(res1[, 4],decreasing = FALSE)
else o1 = NULL
res1 = res1[o1, , drop = F]
res1[, 5] = make.monotone.increasing(as.numeric(res1[, 5]))
}
oo2 = (1:length(r))[r > 0]
res2 = NULL
for (i in oo2) {
res2 = rbind(res2, c(i, geneset.names[i], round(GSA.obj$GSA.scores[i],
4), pvalues.hi[i], fdr.hi[i]))
}
if (!is.null(res2)) o2 = order(res2[, 4], decreasing = FALSE) else o2 = NULL
res2 = res2[o2, , drop = F]
res2[, 5] = make.monotone.increasing(as.numeric(res2[, 5]))
if (length(res1) == 0) {
res1 = NULL
}
if (length(res2) == 0) {
res2 = NULL
}
if ( (length(res1) > 0) & negflag) {
dimnames(res1) = list(NULL, c("Gene_set", "Gene_set_name",
"Score", "p-value", "FDR"))
}
if (length(res2) > 0) {
dimnames(res2) = list(NULL, c("Gene_set", "Gene_set_name",
"Score", "p-value", "FDR"))
}
nsets.neg=NULL
if(negflag){ res1 = res1[as.numeric(res1[, 5]) <= FDRcut,,drop=FALSE ]
nsets.neg = nrow(res1)
if (is.null(res1)) {
nsets.neg = 0
}}
res2 = res2[as.numeric(res2[, 5]) <= FDRcut,,drop=FALSE ]
nsets.pos = nrow(res2)
if (is.null(res2)) {
nsets.pos = 0
}
return(list(FDRcut = FDRcut, negative = res1, positive = res2,
nsets.neg = nsets.neg, nsets.pos = nsets.pos))
}