forked from AstraZeneca-NGS/VarDict
-
Notifications
You must be signed in to change notification settings - Fork 0
/
testsomatic.R
executable file
·42 lines (36 loc) · 1.8 KB
/
testsomatic.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
#!/usr/bin/env Rscript
#args <- commandArgs(trailingOnly = TRUE)
d <- tryCatch( {
d <- read.table( file('stdin'), sep = "\t", header = F, colClasses=c("character", NA, NA, NA, NA, "character", "character", NA, NA, NA, NA, NA, NA, "character", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "character", NA, "character", NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "character", "character", "character", "character"))
}, error = function(e) {
return(NULL)
} )
if (!is.null(d)){
pvalues1 <- vector(mode="double", length=dim(d)[1])
oddratio1 <- vector(mode="double", length=dim(d)[1])
pvalues2 <- vector(mode="double", length=dim(d)[1])
oddratio2 <- vector(mode="double", length=dim(d)[1])
pvalues <- vector(mode="double", length=dim(d)[1])
oddratio <- vector(mode="double", length=dim(d)[1])
for( i in 1:dim(d)[1] ) {
h <- fisher.test(matrix(c(d[i,10], d[i,11], d[i,12], d[i,13]), nrow=2))
pvalues1[i] <- round(h$p.value, 5)
oddratio1[i] <- round(h$estimate, 5)
h <- fisher.test(matrix(c(d[i,28], d[i,29], d[i,30], d[i,31]), nrow=2))
pvalues2[i] <- round(h$p.value, 5)
oddratio2[i] <- round(h$estimate, 5)
tref <- if ( d[i,8] - d[i,9] < 0 ) 0 else d[i,8] - d[i,9]
rref <- if ( d[i,26] - d[i,27] < 0 ) 0 else d[i,26] - d[i,27]
h <- fisher.test(matrix(c(d[i,9], tref, d[i,27], rref), nrow=2), alternative="greater")
pv <- h$p.value
od <- h$estimate
h <- fisher.test(matrix(c(d[i,9], tref, d[i,27], rref), nrow=2), alternative="less")
if ( h$p.value < pv ) {
pv <- h$p.value
od <- h$estimate
}
pvalues[i] <- round(pv, 5)
oddratio[i] <- round(od, 5)
}
write.table(data.frame(d[,1:25], pvalues1, oddratio1, d[,26:43], pvalues2, oddratio2, d[, 44:dim(d)[2]], pvalues, oddratio), file = "", quote = F, sep = "\t", eol = "\n", row.names=F, col.names=F)
}