-
Notifications
You must be signed in to change notification settings - Fork 0
/
msviper_starter.r
79 lines (64 loc) · 2.9 KB
/
msviper_starter.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
args <- commandArgs(TRUE)
runDir <- args[1]
exprsFile <- paste(runDir, args[2], sep = "")
adjFile <- paste(runDir, args[3], sep = "")
phenoFile <- paste(runDir, args[4], sep = "")
context <- args[5]
caseGroups <- gsub('^"|"$', "", args[6])
controlGroups <- gsub('^"|"$', "", args[7])
gesFilter <- args[8]
minSize <- args[9]
bootStrapping <- args[10]
method <- args[11]
shadow <- args[12]
shadowValue <- args[13]
userLib <- args[14]
if (!is.na(userLib)) .libPaths(userLib)
library(viper)
exprs <- as.matrix(read.table(exprsFile, header = TRUE, sep = "\t", row.names = 1, as.is = FALSE, check.names = FALSE))
pData <- read.table(phenoFile, row.names = 1, header = TRUE, sep = "\t")
all(rownames(pData) == colnames(exprs))
metadata <- data.frame(labelDescription = context)
print(metadata)
phenoData <- new("AnnotatedDataFrame", data = pData, varMetadata = metadata)
print(phenoData)
dSet <- ExpressionSet(assayData = exprs, phenoData = phenoData)
print(dSet)
regul <- aracne2regulon(adjFile, dSet, verbose = FALSE)
print(bootStrapping)
cGroups <- c(unlist(strsplit(caseGroups, ",")))
ctlGroups <- c(unlist(strsplit(controlGroups, ",")))
if (bootStrapping == "TRUE") {
print(context)
signature <- bootstrapTtest(dSet, context, cGroups, ctlGroups, verbose = FALSE)
} else {
signature <- rowTtest(dSet, context, cGroups, ctlGroups)
signature <- (qnorm(signature$p.value / 2, lower.tail = FALSE) * sign(signature$statistic))[, 1]
}
nullmodel <- ttestNull(dSet, context, cGroups, ctlGroups, per = 1000, repos = TRUE, verbose = FALSE)
mrs <- msviper(signature, regul, nullmodel, minsize = as.numeric(minSize), ges.filter = as.logical(gesFilter == "TRUE"), verbose = FALSE)
print(mrs)
print(minSize)
print(gesFilter)
summary(mrs, length(mrs$regulon))
if (bootStrapping == "TRUE") {
mrs <- bootstrapmsviper(mrs, c(method))
}
mrs_le <- ledge(mrs)
if (class(signature) == "matrix") signature <- rowMeans(signature)
write.table(signature, file = paste(runDir, "signature.txt", sep = ""), sep = "\t")
slist <- mrs$signature
if (ncol(slist) > 0) slist <- rowMeans(slist)
write.table(slist, file = paste(runDir, "mrsSignature.txt", sep = ""), sep = "\t")
write.table(summary(mrs, length(mrs$regulon)), file = paste(runDir, "result.txt", sep = ""), sep = "\t")
write(paste(mrs_le$ledge, sep = "\t"), file = paste(runDir, "ledges.txt", sep = ""), sep = "\t")
write(names(mrs_le$ledge), file = paste(runDir, "masterRegulons.txt", sep = ""), sep = "\t")
write(paste(lapply(mrs$regulon, function(x) {
names(x$tfmode)
}), sep = "\t"), file = paste(runDir, "regulons.txt", sep = ""), sep = "\t")
if (shadow == "TRUE") {
mrshadow <- shadow(mrs, regulators = as.numeric(shadowValue), verbose = FALSE)
write.table(summary(mrshadow, length(mrshadow$regulon))$msviper.results, file = paste(runDir, "shadowResult.txt", sep = ""), sep = "\t")
write.table(mrshadow$shadow, file = paste(runDir, "shadowPair.txt", sep = ""), sep = "\t")
}
print("complete msviper process.")