Tuesday, March 29, 2011

Prune GWAS data in R

Hansong Wang, our biostats professor here at the Hawaii Cancer Center, generously gave me some R code that goes through a SNP annotation file (i.e. a mapfile) and selects SNPs that are at least a certain specified distance apart. You might want to do this if you're picking a subset of SNPs for PCA, for instance. Plink has an LD pruning feature, but if you can't load your data into PLINK, this poor-man's-pruning based on physical distance (not LD) is a quick solution.

Provide the function with a data frame containing containing column names "chrom" and "position," where the SNPs are ordered by chromosome and position. By default the function selects SNPs that are at least 100kb apart, but you can change this with the optional second argument. The function returns the row indices corresponding to the SNPs you want to keep. Then simply subset your dataset selecting only those row indices and all columns.

1 comment:

  1. Probably faster to pick (one chromosome) like

    pickSNPs1chr <- function(pos, dist=100000)
    {
    res <- logical(length(pos)); i <- 1L; d <- dist
    while ((i <- which.min((pos-d)^2)) != length(pos)) {
    if (pos[i] < d) i <- i + 1
    res[i] <- TRUE
    d <- pos[i] + dist
    }
    which(res)
    }

    ReplyDelete

Creative Commons License
Getting Genetics Done by Stephen Turner is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported License.