Comment rechercher des vallées dans un graphique?

J'examine quelques données de couverture génomique qui sont essentiellement une longue liste (quelques millions de valeurs) d'entiers, chacun indiquant dans quelle mesure (ou "profondément") cette position dans le génome est couverte.

Je voudrais rechercher dans ces données des "vallées", c'est-à-dire des régions nettement "inférieures" à leur environnement.

Notez que la taille des vallées que je recherche peut aller de 50 bases à quelques milliers.

Quel genre de paradigmes recommanderiez-vous d'utiliser pour trouver ces vallées?

METTRE À JOUR

Quelques exemples graphiques pour les données:

MISE À JOUR 2

Définir ce qu'est une vallée est bien sûr l'une des questions avec lesquelles je me bats. Ce sont des évidences pour moi:

mais il y a des situations plus complexes. En général, il y a 3 critères que je considère: 1. La couverture (moyenne? Maximale?) Dans la fenêtre par rapport à la moyenne mondiale. 2. La couverture (...) dans la fenêtre par rapport à son environnement immédiat. 3. Quelle est la taille de la fenêtre: si je vois une couverture très faible pour une courte durée, c'est intéressant, si je vois une couverture très faible pour une longue durée, c'est aussi intéressant, si je vois une couverture légèrement faible pour une courte durée, ce n'est pas vraiment intéressant , mais si je vois une couverture légèrement faible pendant une longue période - c'est ... C'est donc une combinaison de la longueur du sapn et de sa couverture. Plus elle est longue, plus la couverture est élevée et je la considère toujours comme une vallée.

Merci,

Dave

Vous pouvez utiliser une sorte d'approche Monte Carlo, en utilisant par exemple la moyenne mobile de vos données.

Prenez une moyenne mobile des données, en utilisant une fenêtre d'une taille raisonnable (je suppose que c'est à vous de décider de la largeur).

Les traversées de vos données seront (bien sûr) caractérisées par une moyenne inférieure, vous devez donc maintenant trouver un certain "seuil" pour définir le "bas".

Pour ce faire, vous échangez au hasard les valeurs de vos données (par exemple en utilisant sample()) et recalculez la moyenne mobile pour vos données échangées.

Répétez ce dernier passage un nombre de fois raisonnablement élevé (> 5000) et enregistrez toutes les moyennes de ces essais. Donc, essentiellement, vous aurez une matrice de 5000 lignes, une par essai, chacune contenant la moyenne mobile de cet essai.

À ce stade, pour chaque colonne, vous choisissez le quantile de 5% (ou 1% ou tout ce que vous voulez), c'est-à-dire la valeur sous laquelle se trouve seulement 5% des moyennes des données randomisées.

Vous avez maintenant une "limite de confiance" (je ne sais pas si c'est le terme statistique correct) pour comparer vos données originales. Si vous trouvez une partie de vos données inférieure à cette limite, vous pouvez appeler cela un travers.

Bien sûr, gardez à l'esprit que ni cette méthode ni aucune autre méthode mathématique ne pourra jamais vous donner une indication de la signification biologique, bien que je suis sûr que vous en êtes bien conscient.

EDIT - un exemple

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)

Cela vous permettra simplement de trouver graphiquement les régions, mais vous pouvez facilement les trouver en utilisant quelque chose sur les lignes de which(values>limits.5).

J'ignore complètement ces données, mais en supposant que les données sont ordonnées (pas dans le temps, mais par position?), Il est logique d'utiliser des méthodes de séries chronologiques. Il existe de nombreuses méthodes pour identifier les clusters temporels dans les données. Généralement, ils sont utilisés pour trouver des valeurs élevées mais peuvent être utilisés pour des valeurs faibles regroupées. Je pense ici aux statistiques de scan, aux statistiques de somme cumulée (et autres) utilisées pour détecter les épidémies de maladie dans les données de comptage. Des exemples de ces méthodes se trouvent dans le package de surveillance et le package DCluster.

Il existe de nombreuses options pour cela, mais une bonne: vous pouvez utiliser la msExtremafonction dans le msProcesspackage .

Éditer:

Dans l'analyse de la performance financière, ce type d'analyse est souvent effectué en utilisant un concept de "drawdown". Le PerformanceAnalyticspackage a quelques fonctions utiles pour trouver ces vallées . Vous pouvez utiliser le même algorithme ici si vous traitez vos observations comme une série chronologique.

Voici quelques exemples de la façon dont vous pourriez être en mesure d'appliquer cela à vos données (où les "dates" ne sont pas pertinentes mais sont simplement utilisées pour la commande), mais les premiers éléments de l' zooobjet seraient vos données:

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

Une partie des Bioconductor paquets de (par exemple, ShortRead , Biostrings , BSgenome , IRanges , genomeIntervals ) offrent des installations pour traiter des positions du génome ou des vecteurs de couverture, par exemple pour ChIP-seq et l' identification des régions enrichies. En ce qui concerne les autres réponses, je conviens que toute méthode reposant sur des observations ordonnées avec un filtre basé sur un seuil permettrait d'isoler un signal bas dans une bande passante spécifique.

Peut-être que vous pouvez également regarder les méthodes utilisées pour identifier les soi-disant "îles"

Zang, C, Schones, DE, Zeng, C, Cui, K, Zhao, K et Peng, W (2009). Une approche de clustering pour l'identification des domaines enrichis à partir des données ChIP-Seq de modification des histones . Bioinformatics, 25 (15) , 1952-1958.