3
\$\begingroup\$

I frequently want to filter and select elements from a long (10,000s to millions of elements) vector in R. Often, this vector is a column in a data.frame. I am usually interested in rows of the vector that "match" any one of a separate list of 10 to 100 query values. However, since these are experimental data (mass spectrometry data for anyone interested), the target values and the query value(s) do not need to match exactly. They only need to match within a tolerance.

In the field I work in (mass spectrometry), the values are mass-to-charge ratios, or m/z values. There are two systems for specifying the tolerance or precision desired for comparing two numbers. The first is the ppm or "parts per million" tolerance, i.e., a relative measure of precision. The second is in absolute values, which invariably have units of atomic mass units or Daltons. This explains the variable names I use below.

filter_by_mz <- function(query_mzs, target_mzs, ppm = 10, mz_tol = NULL){
    # ensure that users only supply one mass tolerance argument (either ppm or mz_tol)
    if(!missing(ppm) & !missing(mz_tol)){
        stop("Only a single mz tolerance parameter {mz_tol or ppm} can be supplied, not both.")
    }

    # find min and max mzs
    if(missing(mz_tol)){
        min_mzs <- target_mzs * (1 - ppm / 1e6)
        max_mzs <- target_mzs * (1 + ppm / 1e6)
    } else{
        min_mzs <- target_mzs - mz_tol
        max_mzs <- target_mzs + mz_tol
    }

    # allocate boolean matrix
    n_rows <- length(query_mzs)
    n_cols <- length(target_mzs)
    bool_mat <- matrix(FALSE, nrow = n_rows, ncol = n_cols)

    # loop over target mzs to create columns of bool_mat
    for(idx in seq_along(target_mzs)){
        min_mz <- min_mzs[idx]
        max_mz <- max_mzs[idx]
        bool_mat[, idx] <- query_mzs >= min_mz & query_mzs <= max_mz
    }

    # apply any() to rows of matrix
    return(apply(bool_mat, FUN = any, MARGIN = 1))
}

Now my questions:

  1. How can I get rid of the for loop? It seems like a smart way of using apply or a similar function would obviate this loop.

  2. Is the argument handling OK? I wanted users of the function to be able supply either ppm or mz_tol, but not both. Am I using missing() correctly?

  3. Is the design correct? As I mentioned, I frequently will invoke this function on the columns of a data.frame. For example, I'd frequently need to do my_filtered_df <- filter_by_mz(my_df$mz, target_mz, ...). Would it be better to write the function to work on data.frames rather than just vectors?

\$\endgroup\$
1
\$\begingroup\$

Question 1

How can I get rid of the for loop? It seems like a smart way of using apply or a similar function would obviate this loop.

Use outer:

bool_mat <- outer(query_mzs, min_mz, ">=") &
            outer(query_mzs, max_mz, "<=")

You might also know that apply is a disguised for loop: while syntactically shorter, it is just as slow when it comes to computation times. Instead of apply(bool_mat, FUN = any, MARGIN = 1), you could just do:

rowSums(bool_mat) > 0L

where rowSums is a super fast (internal C-compiled) function.

Question 2

Is the argument handling OK? I wanted users of the function to be able supply either ppm or mz_tol, but not both. Am I using missing() correctly?

There are two schools for this kind of situation:

  1. Set the arguments to NULL, i.e., function(query_mzs, target_mzs, ppm = NULL, mz_tol = NULL), and uses is.null for testing if an argument has been set or not. This approach has the advantage that if the user called the function as follows: filter_by_mz(query_mzs, target_mzs, ppm = 10, mz_tol = NULL), then it would work as expected although both arguments have been passed explicitly.
  2. Do not set the argument, i.e. function(query_mzs, target_mzs, ppm, mz_tol) and use missing for testing. Keeping mz_tol = NULL instead of mz_tol alone is not dramatic, though I've not seen much people do that.

Under both designs, it is my humble opinion that you are making it confusing for the user when you make ppm = 10 the default when the user passes nothing. I would not provide any default.

A third approach might be:

function(query_mzs, target_mzs, tol, units = c("ppm", "amu")) {
   units <- match.arg(units)
   if (units == "ppm") {
      ...
   } else {
      ...
   }
   ...
}

Question 3

Is the design correct? As I mentioned, I frequently will invoke this function on the columns of a data.frame. For example, I'd frequently need to do my_filtered_df <- filter_by_mz(my_df$mz, target_mz, ...). Would it be better to write the function to work on data.frames rather than just vectors?

Your function returns a vector of booleans, not a data.frame. So you would have to do:

my_filtered_df <- my_df[filter_by_mz(my_df$mz, target_mz, ...), , drop = FALSE]

or maybe better

my_filtered_df <- subset(my_df, filter_by_mz(mz, target_mz, ...))

which can use multiple criteria separated by & and/or |:

my_filtered_df <- subset(my_df, filter_by_mz(mz1, target_mz1, ...) &
                                filter_by_mz(mz2, target_mz2, ...) |
                                ...)

I think your design is fine this way.

\$\endgroup\$
  • \$\begingroup\$ Thanks! I wound up using this solution. It is noticeably faster. I had no idea that rowSums() existed. \$\endgroup\$ – Curt F. Mar 27 '16 at 2:04
1
\$\begingroup\$

I do quite similar stuff. I find your approach quite lengthy. Especially the construction of the boolean matrix might be a bottleneck.

First I would consider to write a second function - say filter_by_mz_ which just does the filtering and have filter_by_mz just take care of the input and call the new function - just to keep things modular.

In your algorithm you basically compare all masses pairwise. A simpler approach might be to use expand.grid(q=query_mzs, t=target_mzs) this will calculate all pairs. This list can be then easy checked (let's call it gr) using abs(gr$t-gr$q) < ppm * gr$t - that would be fully vectorized.

A bool vector can usually be easy constructed using %in%.

Just FYI I find it much easier to write with dplyr - especially since it is possible to keep much additional info in the table. I often have the constraint that one mz can only be assigned once - extra data is often required to make such decisions.

Typically it looks like this:

Just the setup and some test data:

library(dplyr)

peaks <- data.frame(mz=c(21.01,28.00,30,30.1,35, 21.02), val=c(1,2,3,4,5,6))
info <- data.frame(mz=c(21.008, 28.003, 30.002, 32.02), str=c('wat', 'n2', 'aa' , 'bb'))

The magic happens here:

keyfun <- function(x) sprintf('%.1f',x)
peaks <- mutate(peaks, key=keyfun(mz))
info <- mutate(info, key=keyfun(mz))

inner_join(peaks, info, by='key') %>% filter( abs(mz.x - mz.y) < 0.0005 * mz.x)

Basically I create an temporary key that is too wide - just that I only get to work on combinations that could be. The inner_join will create all the permutations. And in a second step I filter out the ones that actually match.

Of course usually I give more informative names and spend some time labelling stuff. But this should suffice to give an idea of the method.

addition

It came to me that you could avoid the matrix if you could rewrite such that inside the loop you check for all query_mz instead of target_mz - then you could also put the any inside the loop. I suppose it would come down to just switching the variable names in a number of lines. After this transformation the for loop should be replacable by an sapply.

As R stores logicals quite expensive with 4 bytes per logical. That shoud give you extensive savings. (For me these numbers would go into the 100s so I would have >99% less memory usage)

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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