Take the 2-minute tour ×
Code Review Stack Exchange is a question and answer site for peer programmer code reviews. It's 100% free, no registration required.

I have this function containing some loops and a double for loop. What it does is set up a matrix first to store results in (the task at hand is comparing genomes) and then with sliding window, calculate the differences between sequences. It's very slow for very long sequences, I was hoping somebody know how to speed this up, the part causing the bottleneck is:

n1 <- 0

  for(i in seq(windowp1)) {

    n1 <- n1+1
    n2 <- 0
    setTxtProgressBar(prog, n1)
    for(i in seq(1:ncol(Distances))) {
      n2 <- n2+1

      Distances[n1,n2] <- sum(x[pairs[n2,1],c(windowp1[n1]:windowp2[n1])] != x[pairs[n2,2],c(windowp1[n1]:windowp2[n1])])

    }

    Distances2 <- (Distances/winsize)*100

  }

Where windowp1 is a vector of numbers defining the start of each window. Distances is the pre-made matrix for the results, it has as many cols as there are pairwise comparisons of sequences and as many row's as there are sliding window frames, determined by windowp1 and windowp2. The differences between the sequences are calculated simply by the hamming distance. The loop 1 ticks up the window, the second loop then ticks over to do all comparisons between sequences to do in the window, then the first loop ticks to select the next window. The final line makes the hamming distances in the matrix a % of the widowsize.

share|improve this question
add comment

1 Answer 1

Can't see exactly what's going with the sliding windows but it looks like you're calculating the full distance matrix: (i,j) and (j,i). If so, you probably only need half (usually lower) but you can make both with Distances[n1,n2] <- Distances[n2,n1] <- sum ... If this is correct, you'll basically double your speed.

But, if I misunderstood the above, you could still win a little time by reducing the number of index lookup, since these are really slow in R. In particular c(windowp1[n1]:windowp2[n1]) could be pulled out of the loop. Depending on what "x" looks like, you might try "xor(a, b)" rather than "a != b".

Adding a few other changes (eg, drop n1,n2 in favor of the loop variables) but leaving the possible "full matrix" thing I mentioned, here's a suggestion:

distances.ncol <- ncol(Distances)
for(i in seq_along(windowp1)) {
    wp12 <- c(windowp1[i]:windowp2[i])
    setTxtProgressBar(prog, i)

    for(j in seq_len(distances.ncol)) {
        Distances[i,j] <-
            sum(x[pairs[j,1],wp12] != x[pairs[j,2],wp12])
    }
    Distances2 <- (Distances/winsize)*100
}
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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