Geek asks and answers

Manhattan Peak Detection

The attached plot (Manhattan graph) contains the positions of the chromosome of the x-axis from the genome and the Y-axis (p), where p is the p-value associated with the points (variants) from this particular position. enter image description here

I used the following R code to create it (from a space package):

require(gap)
affy <-c(40220, 41400, 33801, 32334, 32056, 31470, 25835, 27457, 22864, 28501, 26273, 
     24954, 19188, 15721, 14356, 15309, 11281, 14881, 6399, 12400, 7125, 6207)
CM <- cumsum(affy)
n.markers <- sum(affy)
n.chr <- length(affy)
test <- data.frame(chr=rep(1:n.chr,affy),pos=1:n.markers,p=runif(n.markers))
oldpar <- par()
par(cex=0.6)
colors <- c("red","blue","green","cyan","yellow","gray","magenta","red","blue","green",          "cyan","yellow","gray","magenta","red","blue","green","cyan","yellow","gray","magenta","red")
mhtplot(test,control=mht.control(colors=colors),pch=19,bg=colors)
> head(test)
  chr pos          p
1   1   1 0.79296584
2   1   2 0.96675136
3   1   3 0.43870076
4   1   4 0.79825513
5   1   5 0.87554143
6   1   6 0.01207523

I'm interested in getting the coordinates of the graph peaks above a certain threshold (-log (p)).

+5
source share
2 answers

If you want indexes of values ​​above the 99th percentile:

# Add new column with log values
test = transform(test, log_p = -log10(test[["p"]]))
# Get the 99th percentile
pct99 = quantile(test[["log_p"]], 0.99)

... and get the values ​​from the source data test:

peaks = test[test[["log_p"]] > pct99,]
> head(peaks)
    chr pos           p    log_p
5     1   5 0.002798126 2.553133
135   1 135 0.003077302 2.511830
211   1 211 0.003174833 2.498279
586   1 586 0.005766859 2.239061
598   1 598 0.008864987 2.052322
790   1 790 0.001284629 2.891222

You can use this with any threshold value. Note that I did not calculate the first derivative, see this question for some pointers:

How to calculate the first time series derivative

, , () . , .

+1

, R

plot (x [, 1], x [, 2])

(x [, 1], x [, 2], labels = row.names(x))

, x [, 1] x ( x [, 2] #your -log10P value

, , enter, # , , #coordinate

coords < - locator (type = "l" )

Coords

0

More articles:


All Articles