forked from tdhock/PeakSegDP
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathPeakSegDP.R
113 lines (109 loc) · 4.01 KB
/
PeakSegDP.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
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
PeakSegDP <- structure(function
### Compute the PeakSeg model on a data.frame of compressed sequence
### reads.
(compressed,
### data.frame with columns chromStart, chromEnd, count.
maxPeaks
### maximum number of peaks to consider.
){
stopifnot(diff(compressed$chromStart) > 0)
count <- compressed$count
weight <- compressed$bases <- with(compressed, chromEnd-chromStart)
stopifnot(is.integer(count))
stopifnot(is.integer(weight))
stopifnot(count >= 0)
stopifnot(weight > 0)
stopifnot(is.integer(maxPeaks))
stopifnot(length(maxPeaks) == 1)
stopifnot(length(count) == length(weight))
maxSegments <- maxPeaks * 2 + 1
stopifnot(maxSegments > 0)
stopifnot(maxSegments <= nrow(compressed))
fit <- cDPA(count, weight, maxSegments)
segment.ends <- getPath(fit)
results <- list(breaks=list())
for(peaks in 0:maxPeaks){
peak.list <- list()
segments <- as.integer(peaks*2 + 1)
model.i <- peaks * 2 + 1
last.i <- as.integer(segment.ends[model.i, 1:model.i])
not.feasible <- any(last.i < 1 | is.na(last.i))
if(not.feasible){
s <- ifelse(peaks==1, "", "s")
warning("infeasible model with ", peaks, " peak", s)
}else{
break.after <- last.i[-model.i]
first.i <- as.integer(c(1, break.after+1))
model.error <- 0
if(length(break.after)){
results$breaks[[paste(peaks)]] <-
data.frame(peaks, segments, break.after,
chromEnd=compressed$chromEnd[break.after])
}
for(segment.i in seq_along(last.i)){
status <- ifelse(segment.i %% 2, "background", "peak")
first <- first.i[[segment.i]]
last <- last.i[[segment.i]]
seg.data <- compressed[first:last,]
count.num <- as.numeric(seg.data$count)
bases.num <- as.numeric(seg.data$bases)
seg.mean <- sum(count.num*bases.num)/sum(bases.num)
model.error <- model.error + with(seg.data, {
PoissonLoss(count, seg.mean, bases)
})
chromStart <- seg.data$chromStart[1]
chromEnd <- seg.data$chromEnd[nrow(seg.data)]
results$segments[[paste(peaks, segment.i)]] <-
data.frame(mean=seg.mean,
first,
last,
chromStart,
chromEnd,
status,
peaks,
segments)
if(status == "peak"){
peak.list[[paste(segment.i)]] <-
data.frame(first, last,
chromStart, chromEnd,
peaks, segments)
}
}#segment.i
results$peaks[[as.character(peaks)]] <- do.call(rbind, peak.list)
results$error[[as.character(peaks)]] <-
data.frame(segments=model.i, peaks, error=model.error)
}#!not.feasible
}#peaks
no.peaks <-
data.frame(first=integer(), last=integer(),
chromStart=integer(), chromEnd=integer(),
peaks=integer(), segments=integer())
results$peaks <- c(list("0"=no.peaks), results$peaks)
results$error <- do.call(rbind, results$error)
results$segments <- do.call(rbind, results$segments)
results$breaks <- do.call(rbind, results$breaks)
results
}, ex=function(){
library(PeakSegDP)
data(chr11ChIPseq, envir=environment())
one <- subset(chr11ChIPseq$coverage, sample.id=="McGill0002")[10000:12000,]
fit <- PeakSegDP(one, 3L)
if(interactive() && require(ggplot2)){
ggplot()+
geom_step(aes(chromStart/1e3, count), data=one)+
geom_segment(aes(chromStart/1e3, mean,
xend=chromEnd/1e3, yend=mean),
data=fit$segments, color="green")+
geom_segment(aes(chromStart/1e3, 0,
xend=chromEnd/1e3, yend=0),
data=subset(fit$segments, status=="peak"),
size=3, color="deepskyblue")+
theme_bw()+
theme(panel.margin=grid::unit(0, "cm"))+
facet_grid(peaks ~ ., scales="free", labeller=function(df){
s <- ifelse(df$peaks==1, "", "s")
df$peaks <- paste0(df$peaks, " peak", s)
df
})
}
})