我正在检查一些基因组覆盖数据，这些数据基本上是一个很长的整数列表（几百万个值），每个整数都说明基因组中这个位置的覆盖程度（或“深度”）。
我想在这个数据中寻找“山谷”，即明显“低于”周围环境的区域。
请注意，我正在寻找的山谷的大小可能从 50 个碱基到数千个不等。
你会推荐使用什么样的范式来找到那些山谷？
更新
数据的一些图形示例：
更新 2
定义什么是山谷当然是我正在努力解决的问题之一。这些对我来说是显而易见的：
但还有一些更复杂的情况。一般来说，我考虑了 3 个标准： 1. 窗口中相对于全局平均值的（平均？最大？）覆盖率。2. 窗口中相对于其周围环境的 (...) 覆盖范围。3. 窗口有多大：如果我看到短时间内覆盖率很低，这很有趣，如果我看到很长一段时间内覆盖率很低，那也很有趣，如果我看到短时间内覆盖率很低，那不是很有趣，但是如果我看到长时间的覆盖率略低 - 它是......所以它是 sapn 的长度和它的覆盖率的组合。它越长，我让覆盖范围越高，仍然认为它是一个山谷。
谢谢，
戴夫
您可以使用某种蒙特卡罗方法，例如使用数据的移动平均值。
使用合理大小的窗口对数据进行移动平均（我想这取决于您决定多宽）。
您的数据中的通过率（当然）将以较低的平均值为特征，因此现在您需要找到一些“阈值”来定义“低”。
为此，您随机交换数据的值（例如使用sample()）并重新计算交换数据的移动平均值。
重复这最后一段相当高的次数（> 5000）并存储这些试验的所有平均值。所以基本上你将有一个包含 5000 行的矩阵，每个试验一个，每个包含该试验的移动平均值。
此时，您为每一列选择 5%（或 1% 或任何您想要的）分位数，即随机数据均值的 5% 的值。
您现在有一个“置信度限制”（我不确定这是否是正确的统计术语）来比较您的原始数据。如果您发现部分数据低于此限制，那么您可以称其为通过。
当然，请记住，这个或任何其他数学方法都不能给你任何生物学意义的迹象，尽管我相信你很清楚这一点。
编辑 - 一个例子
require(ares) # for the ma (moving average) function
# Some data with peaks and throughs
values <- cos(0.12 * 1:100) + 0.3 * rnorm(100)
plot(values, t="l")
# Calculate the moving average with a window of 10 points
mov.avg <- ma(values, 1, 10, FALSE)
numSwaps <- 1000
mov.avg.swp <- matrix(0, nrow=numSwaps, ncol=length(mov.avg))
# The swapping may take a while, so we display a progress bar
prog <- txtProgressBar(0, numSwaps, style=3)
for (i in 1:numSwaps)
{
# Swap the data
val.swp <- sample(values)
# Calculate the moving average
mov.avg.swp[i,] <- ma(val.swp, 1, 10, FALSE)
setTxtProgressBar(prog, i)
}
# Now find the 1% and 5% quantiles for each column
limits.1 <- apply(mov.avg.swp, 2, quantile, 0.01, na.rm=T)
limits.5 <- apply(mov.avg.swp, 2, quantile, 0.05, na.rm=T)
# Plot the limits
points(limits.5, t="l", col="orange", lwd=2)
points(limits.1, t="l", col="red", lwd=2)
这将只允许您以图形方式查找区域，但您可以使用which(values>limits.5).
我完全不知道这些数据，但假设数据是有序的（不是按时间排序，而是按位置排序？）使用时间序列方法是有意义的。有很多方法可以识别数据中的时间簇。通常，它们用于查找高值，但可用于组合在一起的低值。我在这里考虑的是用于检测计数数据中疾病爆发的扫描统计数据、累积总和统计数据（和其他数据）。这些方法的示例在监视包和 Dcluster 包中。
有很多选择，但有一个很好的选择：您可以使用packagemsExtrema中的函数。msProcess
编辑：
在财务业绩分析中，这种分析通常使用“回撤”概念进行。该PerformanceAnalytics软件包有一些有用的功能可以找到这些山谷。如果您将观察结果视为时间序列，则可以在此处使用相同的算法。
以下是一些示例，说明如何将其应用于数据（其中“日期”无关紧要，仅用于排序），但zoo对象中的第一个元素将是您的数据：
library(PerformanceAnalytics)
x <- zoo(cumsum(rnorm(50)), as.Date(1:50))
findDrawdowns(x)
table.Drawdowns(x)
chart.Drawdown(x)
一些Bioconductor的软件包（例如ShortRead、Biostrings、BSgenome、IRanges、genomeIntervals）提供了处理基因组位置或覆盖向量的工具，例如用于ChIP-seq和识别富集区域。至于其他答案，我同意任何依赖于有序观察和一些基于阈值的滤波器的方法都允许在特定带宽内隔离低信号。
或许你也可以看看用来识别所谓“孤岛”的方法
Zang, C, Schones, DE, Zeng, C, Cui, K, Zhao, K 和 Peng, W (2009)。一种从组蛋白修饰 ChIP-Seq 数据中识别富集域的聚类方法。生物信息学，25(15)，1952-1958。