Dear all,

I would like to filter the SNPs of a VCF file so to extract only SNPs distant of X kb between each other. I do not see any built-in function to perform this kind of filtering so aimed at writing a filter rule to achieve it and pass it through the filterVcf function of the VariantAnnotation package.

However, the kind of function I was considering would require a for loop which would compare for each SNP entry its position to the last_accepted_position so that the SNP is given a logical TRUE if last_accepted_position - position > X or a FALSE if not. At each iteration of the loop, last_accepted_position would be incremented to last_accepted_position + position if the SNP is TRUE, or would remain the same if SNP is FALSE.

How can I pass such loop using filterRules and filterVcf of the VariantAnnotation package?

Many thanks for your help,

Rémi

You should be able to do this without a loop.  Assuming that you have your SNPs in a VRanges or VCF object:

pos = start(vr)
inter_snp_distance = diff(pos)
idx = ( inter_snp_distance >= your_threshold )

In the context of using filterVcf(), and using Julian's approach, the idea is to create a 'closure' (in this case, a function that returns a function) that can remember state (last_pos) between invocations. Here's  my filter factory (notice the use of <<- to update last_pos in the environment in which the filter function is defined, rather than in the body of the function itself).

factory <- function(threshold) {
    last_pos <- 0L
    function(x) {
        dist <- diff(c(last_pos, start(x)))
        keep <- dist > threshold
        if (any(keep))
            last_pos <<- tail(start(x)[keep], 1)
        keep
    }
}
filt <- FilterRules(list(isDistant = factory(100)))

and application

fl <- system.file(package="VariantAnnotation", "extdata", "chr22.vcf.gz")
destination <- tempfile()
filtered <- filterVcf(fl, "hg19", destination, filters=filt)

A little more challenging is dealing with changes in chromosome, which requires that the last position be reset.

Ok so given that there are more than one way to pick up the SNPs we should do it in a really random fashion. The "memory-keeping initial_position" approach introduces a systematic bias. Imagine for example that you have SNPs at 100 consecutive positions, say at positions 101:200, then you're going to pick up SNPs at 101, 121, 141, 161, and 181, which is not good.

So here is an algo that is super fast and extracts a random subset of SNPS with a space of at least 'min.space' between them. The returned subset is "maximal" in the sense that no SNP can be added to it (i.e. the SNPs that didn't make it into it are too close):

library(IRanges)

clusterSNPs <- function(snp_pos, min.space=20)
{
   clusters_start_idx <-
       which(c(TRUE, diff(snp_pos) >= min.space))
   clusters_end_idx <-
       c((clusters_start_idx-1L)[-1L], length(snp_pos))
   PartitioningByEnd(clusters_end_idx)
}

selectOneFromEachCluster <- function(x)
{
    x_eltlens <- unname(elementLengths(x))
    i0 <- sample(max(x_eltlens), length(x_eltlens),
                 replace=TRUE) %% x_eltlens
    end(x) - i0
}

### Assumes that 'snp_pos' is sorted. Returns an index into the
### original set.
extractSpacedSNPs <- function(snp_pos, min.space=20)
{
    ans <- integer(0)
    remaining_idx <- seq_along(snp_pos)
    while (length(remaining_idx) != 0L) {
        remaining_pos <- snp_pos[remaining_idx]
        ## 1. Cluster the SNPs:
        clusters <- clusterSNPs(remaining_pos, min.space)
        ## 2. Pick up one SNP in each cluster and add it to
        ##    the result:
        selection <- selectOneFromEachCluster(clusters)
        ans <- c(ans, remaining_idx[selection])
        ## 3. Drop remaining SNPs that are at a distance less
        ##    than 'min.space' from the SNPs we just added to
        ##    the result:
        remaining_clusters <- relist(remaining_pos, clusters)
        keep <- which(unlist(abs(remaining_clusters -
                                 remaining_pos[selection]) >=
                                                  min.space))
        remaining_idx <- remaining_idx[keep]
    }
    sort(ans)
}

Then:

snp_pos <- 101:200
idx <- extractSpacedSNPs(snp_pos)
snp_pos[idx]
# [1] 109 136 163 191

Cheers,

H.

Dear Julian and Martin, thank you so much for your messages and propositions.

Julian, I wrote this script initially but the problem of the diff function is that it computes the difference between consecutive pairs of array elements. Thus, if I have for instance an array

0 50 60 70 100

And a threshold of 20, my SNP will be

TRUE FALSE FALSE TRUE

However using a memory-keeping initial_position variable, I would get

TRUE TRUE FALSE TRUE TRUE

I believe using the diff function would underestimate the number of relevant SNPs.

Best regards,

Rémi

Many thanks Martin and Hervé!

Wishing you a very happy new year,

Rémi