Jak szukać dolin na wykresie?

Badam niektóre dane pokrycia genomowego, które są w zasadzie długą listą (kilka milionów wartości) liczb całkowitych, z których każda mówi, jak dobrze (lub „głęboka”) pozycja w genomie jest objęta.

Chciałbym poszukać w tych danych „dolin”, czyli regionów znacznie „niższych” niż otaczające je środowisko.

Zauważ, że rozmiar dolin, których szukam, może wynosić od 50 zasad do kilku tysięcy.

Jakich paradygmatów poleciłbyś zastosować do znalezienia tych dolin?

AKTUALIZACJA

Kilka graficznych przykładów danych:

AKTUALIZACJA 2

Zdefiniowanie, czym jest dolina, jest oczywiście jednym z pytań, z którymi się zmagam. Są to dla mnie oczywiste:

ale istnieją bardziej złożone sytuacje. Zasadniczo rozważam 3 kryteria: 1. Zakres (średni? Maksymalny?) W oknie w stosunku do średniej globalnej. 2. Zasięg (...) w oknie w odniesieniu do jego bezpośredniego otoczenia. 3. Jak duże jest okno: jeśli widzę bardzo niski zasięg dla krótkiego zakresu, jest to interesujące, jeśli widzę bardzo niskie pokrycie dla dużego zakresu, jest również interesujące, jeśli widzę lekko niskie pokrycie dla krótkiego zakresu, to nie jest naprawdę interesujące , ale jeśli widzę lekko niski zasięg przez długi okres - to jest… Więc jest to kombinacja długości sapn i jego zasięgu. Im jest dłuższy, tym większy pozwalam na zasięg i nadal uważam go za dolinę.

Dzięki,

Dave

Możesz użyć pewnego rodzaju podejścia Monte Carlo, wykorzystując na przykład średnią ruchomą twoich danych.

Weź średnią ruchomą danych, korzystając z okna o rozsądnym rozmiarze (myślę, że to od Ciebie zależy, jak szerokie).

Przejścia w twoich danych będą (oczywiście) charakteryzować się niższą średnią, więc teraz musisz znaleźć jakiś „próg”, aby zdefiniować „niski”.

W tym celu losowo zamieniasz wartości swoich danych (np. Używając sample()) i ponownie obliczasz średnią ruchomą dla zamienionych danych.

Powtórz ten ostatni fragment stosunkowo dużo razy (> 5000) i zapisz wszystkie średnie z tych prób. Zasadniczo będziesz miał macierz zawierającą 5000 linii, po jednej na próbę, z których każda zawiera średnią ruchomą dla tej próby.

W tym momencie dla każdej kolumny wybierasz kwantyl 5% (lub 1% lub cokolwiek chcesz), czyli wartość, pod którą leży tylko 5% średnich danych losowych.

Masz teraz „limit ufności” (nie jestem pewien, czy jest to poprawny termin statystyczny), z którym możesz porównać swoje oryginalne dane. Jeśli znajdziesz część swoich danych, która jest niższa niż ten limit, możesz nazwać to poprzez.

Oczywiście, pamiętajcie, że ani ta, ani żadna inna metoda matematyczna nigdy nie dałaby wam żadnego wskazania o znaczeniu biologicznym, chociaż jestem pewien, że jesteście tego świadomi.

EDYCJA - przykład

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)

To pozwoli ci tylko graficznie znaleźć regiony, ale możesz je łatwo znaleźć, używając czegoś na linii which(values>limits.5).

Jestem całkowicie nieświadomy tych danych, ale zakładając, że dane są uporządkowane (nie w czasie, ale według pozycji?) Warto skorzystać z metod szeregów czasowych. Istnieje wiele metod identyfikowania klastrów czasowych w danych. Zasadniczo są one używane do znajdowania wysokich wartości, ale można ich używać do grupowania niskich wartości. Mam tu na myśli statystyki skanowania, statystyki sumy zbiorczej (i inne) używane do wykrywania wybuchów chorób w danych zliczania. Przykłady tych metod znajdują się w pakiecie nadzoru i pakiecie DCluster.

Jest na to wiele opcji, ale jedna dobra: możesz użyć msExtremafunkcji w msProcesspakiecie .

Edytować:

W analizie wyników finansowych tego rodzaju analiza jest często przeprowadzana przy użyciu koncepcji „wypłaty”. PerformanceAnalyticsPakiet ma kilka przydatnych funkcji, aby znaleźć tych dolin . Możesz użyć tego samego algorytmu tutaj, jeśli potraktujesz swoje obserwacje jako szereg czasowy.

Oto kilka przykładów tego, w jaki sposób możesz zastosować to do swoich danych (gdzie „daty” są nieistotne, ale służą jedynie do zamawiania), ale pierwszymi elementami w zooobiekcie byłyby twoje dane:

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

Niektóre pakiety Bioconductor (np. ShortRead , Biostrings , BSgenome , IRanges , GenomeIntervals ) oferują udogodnienia do radzenia sobie z pozycjami genomu lub wektorami pokrycia, np. Do sekwencji ChIP i identyfikacji wzbogaconych regionów. Co do innych odpowiedzi, zgadzam się, że każda metoda polegająca na uporządkowanych obserwacjach z pewnym filtrem opartym na progach pozwoliłaby na izolację niskiego sygnału w określonym paśmie.

Być może możesz także spojrzeć na metody stosowane do identyfikacji tak zwanych „wysp”

Zang, C, Schones, DE, Zeng, C, Cui, K, Zhao, K i Peng, W (2009). Podejście grupujące do identyfikacji wzbogaconych domen na podstawie danych ChIP-Seq modyfikacji histonów . Bioinformatics, 25 (15) , 1952–1958.