11
$\begingroup$

I'm examining some genomic coverage data which is basically a long list (a few million values) of integers, each saying how well (or "deep") this position in the genome is covered.

I would like to look for "valleys" in this data, that is, regions which are significantly "lower" than their surrounding environment.

Note that the size of the valleys I'm looking for may range from 50 bases to a few thousands.

What kind of paradigms would you recommend using to find those valleys?

UPDATE

Some graphical examples for the data: alt text alt text

UPDATE 2

Defining what is a valley is of course one of the question I'm struggling with. These are obvious ones for me: alt text alt text

but there some more complex situations. In general, there are 3 criteria I consider: 1. The (average? maximal?) coverage in the window with respect to the global average. 2. The (...) coverage in the window with respect to its immediate surrounding. 3. How large is the window: if I see very low coverage for a short span it is interesting, if I see very low coverage for a long span it's also interesting, if I see mildly low coverage for a short span it's not really interesting, but if I see mildly low coverage for a long span - it is.. So it's a combination of the length of the sapn and it's coverage. The longer it is, the higher I let the coverage be and still consider it a valley.

Thanks,

Dave

$\endgroup$
8
  • $\begingroup$ Could you provide a small data sample? $\endgroup$
    – Shane
    Sep 24, 2010 at 16:20
  • $\begingroup$ @Shane see update $\endgroup$
    – David B
    Sep 24, 2010 at 16:43
  • $\begingroup$ @David Thanks. As both the answers imply, time series analysis can be applied here since you have ordered observations. $\endgroup$
    – Shane
    Sep 24, 2010 at 16:50
  • $\begingroup$ This is kind of hard to answer without know exactly what you're looking for. Can you maybe circle the points on the plots that you're looking to capture? What do you consider a "valley"? how low does it have to go and what are you looking to return? It's hard to formulate a solution without knowing the question, ie thresholds and such. $\endgroup$
    – Falmarri
    Sep 24, 2010 at 17:48
  • $\begingroup$ @ Shane♦ Thank you. As I have no experience with time-series analysis also, could you leave a few pointers of where should I start? $\endgroup$
    – David B
    Sep 24, 2010 at 18:39

4 Answers 4

5
$\begingroup$

You could use some sort of Monte Carlo approach, using for instance the moving average of your data.

Take a moving average of the data, using a window of a reasonable size (I guess it's up to you deciding how wide).

Throughs in your data will (of course) be characterized by a lower average, so now you need to find some "threshold" to define "low".

To do that you randomly swap the values of your data (e.g. using sample()) and recalculate the moving average for your swapped data.

Repeat this last passage a reasonably high amount of times (>5000) and store all the averages of these trials. So essentially you will have a matrix with 5000 lines, one per trial, each one containing the moving average for that trial.

At this point for each column you pick the 5% (or 1% or whatever you want) quantile, that is the value under which lies only 5% of the means of the randomized data.

You now have a "confidence limit" (I'm not sure if that is the correct statistical term) to compare your original data with. If you find a part of your data that is lower than this limit then you can call that a through.

Of course, bare in mind that not this nor any other mathematical method could ever give you any indication of biological significance, although I'm sure you're well aware of that.

EDIT - an example

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)

This will just allow you to graphically find the regions, but you can easily find them using something on the lines of which(values>limits.5).

$\endgroup$
3
  • $\begingroup$ Obviously you can apply the same approach using something else than the moving average, this was just to give an idea. $\endgroup$
    – nico
    Sep 24, 2010 at 20:38
  • $\begingroup$ +1 Thank you very much, nico. Let me see if I got you right: at the end, this is basically like setting some global threshold and defining any point with value < threshold as a part of a valley. The sampling etc. is just used to get some meaningful measure (quantile) to set the threshold. Why can't we use a single threshold for the entire points, I mean, if we did enough simulations we would get straight (read and yellow) lines. Also, correct me if I'm mistaken, but this does not take into account the surrounding environment but examines the absolute value of each point. $\endgroup$
    – David B
    Sep 24, 2010 at 21:29
  • $\begingroup$ @David B: of course, you could use a global threshold and that would probably save you some calculation time. I guess choosing something like 1/3 of the global mean could be a start. This swapping process is probably more helpful if you use some other statistics than the moving average, it was mostly to give an idea. Anyway the moving average will take into account the surrounding, in the example it will take into account a window of 10 points. $\endgroup$
    – nico
    Sep 24, 2010 at 21:35
4
$\begingroup$

I'm completely ignorant of these data, but assuming the data are ordered (not in time, but by position?) it makes sense to make use of time series methods. There are lots of methods for identifying temporal clusters in data. Generally they are used to find high values but can be used for low values grouped together. I'm thinking here of scan statistics, cumulative sum statistics (and others) used to detect disease outbreaks in count data. Examples of these methods are in the surveillance package and the DCluster package.

$\endgroup$
1
  • $\begingroup$ @cxr Thank for your response. I has a look at surveillance and DCluster , but could you please be a bit more specific? They are both relatively large packages and their aim seems quite specific. I'm not sure where to begin. $\endgroup$
    – David B
    Sep 24, 2010 at 18:37
2
$\begingroup$

There are many options for this, but one good one: you can use the msExtrema function in the msProcess package.

Edit:

In financial performance analysis, this kind of analysis is often performed using a "drawdown" concept. The PerformanceAnalytics package has some useful functions to find these valleys. You could use the same algorithm here if you treat your observations as a time series.

Here are some examples of how you might be able to apply this to your data (where the "dates" are irrelevant but just used for ordering), but the first elements in the zoo object would be your data:

library(PerformanceAnalytics)
x <- zoo(cumsum(rnorm(50)), as.Date(1:50))
findDrawdowns(x)
table.Drawdowns(x)
chart.Drawdown(x)
$\endgroup$
6
  • $\begingroup$ Thank you Shane, but this seem to find local minima (or maxima) - i.e. a single point in a region. My data ( as any biological data) IS NOISY> I don't really care about point minima themselves but about larger regions which are low. $\endgroup$
    – David B
    Sep 24, 2010 at 15:53
  • $\begingroup$ If you have local maximum and minimum points, you can easily calculate the differences. So you want to know instances when the differences are both large in magnitude and in "duration"? Is this time series data? $\endgroup$
    – Shane
    Sep 24, 2010 at 15:59
  • $\begingroup$ @david Perhaps, you can iteratively use this function. Use the function to identify a minima. Drop that point and surrounding points (say x points within some tolerance level). You can choose a tolerance level (e.g., +- 10 counts) which would define a flat region for your application. Find a new minima on the new dataset. Will that work? $\endgroup$
    – user28
    Sep 24, 2010 at 16:01
  • $\begingroup$ @shane The analogy that comes to mind is that of valleys in a mountainous region. I think the goal is to identify all the valleys and the issue is some valleys are 'deeper' and some are 'shallow' relative to the mountains. $\endgroup$
    – user28
    Sep 24, 2010 at 16:02
  • $\begingroup$ @Shane It's not a time series, these are coordinate along the genome (chromosome). $\endgroup$
    – David B
    Sep 24, 2010 at 16:03
2
$\begingroup$

Some of the Bioconductor's packages (e.g., ShortRead, Biostrings, BSgenome, IRanges, genomeIntervals) offer facilities for dealing with genome positions or coverage vectors, e.g. for ChIP-seq and identifying enriched regions. As for the other answers, I agree that any method relying on ordered observations with some threshold-based filter would allow to isolate low signal within a specific bandwith.

Maybe you can also look at the methods used to identify so-called "islands"

Zang, C, Schones, DE, Zeng, C, Cui, K, Zhao, K, and Peng, W (2009). A clustering approach for identification of enriched domains from histone modification ChIP-Seq data. Bioinformatics, 25(15), 1952-1958.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.