-
Notifications
You must be signed in to change notification settings - Fork 1
/
ciceroFilter
executable file
·69 lines (53 loc) · 1.79 KB
/
ciceroFilter
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
#! /usr/bin/env Rscript
#
# Vivek Rai
# GPL-3
#
#
doc <- "
ciceroFilter.
Provides a few options to filter cicero connections output.
chrA_startA_endA chrB_startB_endB score
Usage: ciceroFilter [--keeppair] [--keepna] [--threshold=T] <cicero> <output>
Options:
<cicero> Cicero connections
<output> Output file
--keeppair Cicero outputs peak-pairs in both directions (A-B, B-A), keep as is. [default: False]
--keepna Keep peak-pairs with NA/invalid score [default: False]
--threshold=T Threshold to use for filtering peak-pairs [default: 0.05]
"
opts <- docopt::docopt(doc)
conns <- data.table::fread(opts$cicero,
header = F,
col.names = c("Peak1", "Peak2", "coaccess")
)
write(glue::glue("Found {nrow(conns)} peak-pairs.."), stdout())
#
# Given a peak, calculates the mid-point
#
calculate_mid <- function(peaks) {
tmp <- stringr::str_split(peaks, "_", simplify = T)
# mid point is START + offset
mid <- as.integer(tmp[, 2]) + as.integer((as.integer(tmp[, 3]) - as.integer(tmp[, 2])) / 2)
return(mid)
}
# Remove NAs from the co-access score otherwise `compare_connections` throws error
if (!opts$keepna) {
conns <- conns[!is.na(coaccess), ]
}
if (!opts$keeppair) {
conns$distance <- calculate_mid(conns$Peak2) - calculate_mid(conns$Peak1)
conns <- conns[distance > 0, ]
}
if (!is.null(opts$threshold)) {
threshold <- as.numeric(opts$threshold)
if (threshold < -1 | threshold > 1) {
stop("Error: Specify threshold between -1 to 1. Exit.")
}
conns <- conns[coaccess > threshold, ]
}
write(glue::glue("Left with {nrow(conns)} peak-pairs.."), stderr())
dir.create(dirname(opts$output), recursive = T, showWarnings = F)
data.table::fwrite(conns[, 1:3], file = opts$output, quote = F, row.names = F,
col.names = F, sep = "\t")