¿Cómo buscar valles en un gráfico?

Question

Estoy examinando algunos datos de cobertura genómica que es básicamente una larga lista (unos pocos millones de valores) de números enteros, cada uno de los cuales dice cuán bien (o "profundo") está cubierta esta posición en el genoma.

Me gustaría buscar "valles" en estos datos, es decir, regiones que son significativamente "más bajas" que su entorno.

Tengan en cuenta que el tamaño de los valles que busco puede variar desde 50 bases hasta unos pocos miles.

¿Qué clase de paradigmas recomendaría usar para encontrar esos valles?

ACTUALIZACIÓN

Algunos ejemplos gráficos de los datos:

ACTUALIZACIÓN 2

Definir lo que es un valle es, por supuesto, una de las cuestiones con las que estoy luchando. Estas son obvias para mí:

pero hay algunas situaciones más complejas. En general, hay 3 criterios que considero: 1. La (media? máxima?) cobertura en la ventana con respecto a la media global. 2. La (...) cobertura en la ventana con respecto a su entorno inmediato. 3. Qué tan grande es la ventana: si veo una cobertura muy baja para un período corto es interesante, si veo una cobertura muy baja para un período largo también es interesante, si veo una cobertura ligeramente baja para un período corto es no realmente interesante, pero si veo una cobertura ligeramente baja durante un largo período, es Así que es una combinación de la longitud de la savia y su cobertura. Cuanto más larga es, más alta dejo que sea la cobertura y sigo considerándola un valle.

Gracias,

Dave

Answer 1

Podría utilizar algún tipo de enfoque de Montecarlo, utilizando, por ejemplo, la media móvil de sus datos.

Haz una media móvil de los datos, utilizando una ventana de un tamaño razonable (supongo que depende de ti decidir cómo de ancha).

Los datos de la empresa se caracterizarán (por supuesto) por una media más baja, así que ahora hay que encontrar algún "umbral" para definir lo que es "bajo".

Para ello, intercambia aleatoriamente los valores de tus datos (por ejemplo, utilizando sample() ) y vuelva a calcular la media móvil para sus datos intercambiados.

Repite este último pasaje una cantidad razonablemente alta de veces (>5000) y guarda todos los promedios de estos ensayos. Así que esencialmente tendrás una matriz con 5000 líneas, una por ensayo, cada una de las cuales contiene la media móvil de ese ensayo.

En este punto, para cada columna, se elige el cuantil del 5% (o del 1% o lo que se quiera), es decir, el valor por debajo del cual se encuentra sólo el 5% de las medias de los datos aleatorios.

Ahora tiene un "límite de confianza" (no estoy seguro de que sea el término estadístico correcto) con el que comparar sus datos originales. Si encuentras una parte de tus datos que es inferior a este límite, entonces puedes llamarlo "a través".

Por supuesto, tenga en cuenta que ni este ni ningún otro método matemático podría darle ninguna indicación de la importancia biológica, aunque estoy seguro de que es muy consciente de ello.

EDITAR - un ejemplo

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)

Esto sólo le permitirá encontrar gráficamente las regiones, pero puede encontrarlas fácilmente utilizando algo parecido a which(values>limits.5) .

Answer 2

Desconozco por completo estos datos, pero asumiendo que los datos están ordenados (no en el tiempo, sino por posición ) tiene sentido hacer uso de métodos de series temporales. Hay muchos métodos para identificar grupos temporales en los datos. Generalmente se utilizan para encontrar valores altos, pero pueden utilizarse para valores bajos agrupados. Estoy pensando en los estadísticos de barrido, los estadísticos de suma acumulada (y otros) utilizados para detectar brotes de enfermedades en datos de recuento. Hay ejemplos de estos métodos en el paquete de vigilancia y en el paquete DCluster.

Answer 3

Algunos de los Bioconductor 's paquetes (por ejemplo, Lectura corta , Biostrings , BSgenome , IRanges , genomeIntervals ) ofrecen facilidades para tratar las posiciones del genoma o los vectores de cobertura, por ejemplo, para ChIP-seq e identificar las regiones enriquecidas. En cuanto a las demás respuestas, estoy de acuerdo en que cualquier método que se base en observaciones ordenadas con algún filtro basado en el umbral permitiría aislar la señal baja dentro de un ancho de banda específico.

Tal vez también puedas mirar los métodos utilizados para identificar las llamadas "islas"

Zang, C, Schones, DE, Zeng, C, Cui, K, Zhao, K, y Peng, W (2009). A enfoque de clustering para la identificación de dominios enriquecidos a partir de datos modificación de histonas ChIP-Seq . Bioinformática, 25(15) , 1952-1958.

Answer 4

Hay muchas opciones para esto, pero una buena: puedes usar el msExtrema en la función msProcess paquete .

Editar:

En el análisis de los resultados financieros, este tipo de análisis suele realizarse utilizando el concepto de "drawdown". El PerformanceAnalytics paquete tiene algunos funciones útiles para encontrar estos valles . Podría utilizar el mismo algoritmo aquí si trata sus observaciones como una serie temporal.

Estos son algunos ejemplos de cómo podría aplicar esto a sus datos (en los que las "fechas" son irrelevantes, sino que sólo se utilizan para ordenar), pero los primeros elementos del zoo objeto serían sus datos:

library(PerformanceAnalytics)
x <- zoo(cumsum(rnorm(50)), as.Date(1:50))
findDrawdowns(x)
table.Drawdowns(x)
chart.Drawdown(x)