Prune GWAS data in R
[This article was first published on Getting Genetics Done, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
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.Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
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.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#------------------------------------------------------------------------------------- | |
# Select a set of SNPs > dist bp (default 100kb) apart | |
# Map is a matrix with colnames "chrom" and "position". | |
# Matrix MUST BE sorted by chrom,position | |
# Can have more columns but the colnames "chrom" and "position" can't be changed | |
# The function returns a vector of row indices corresponding | |
# to the rows you should pick for your subset. | |
#------------------------------------------------------------------------------------- | |
pickSNPs<-function(map, dist = 100000) { | |
t=as.data.frame(table(map$chrom)) | |
vec = map$position | |
subs = c(1,rep(NA,nrow(map)-1)); # length(subs) = nrow, but the 1st element is 1 => always select the 1st snp | |
for (k in 1:nrow(t)) { # t: count of SNPs per chr | |
if (k==1) i=1 else i=sum(t[1:(k-1),2])+1; # the 1st SNP on each ch | |
subs[i] = i | |
stop=sum(t[1:k,2]) | |
while (i<stop) { | |
for (j in (i+1):stop) { | |
if ((vec[j]-vec[i]) > dist) { | |
#cat(i, vec[i], j, vec[j],vec[j]-vec[i], x[j],'\n'); | |
subs[j]= j; | |
i=j; | |
next; # jump out of loop | |
} else if (j==stop) i=stop | |
} | |
} | |
} | |
subs[!is.na(subs)] # row number of selected SNPs | |
} | |
#------------------------------------------------------------------------------------- | |
map <- read.table("mapfile.txt",header=TRUE) | |
# mapfile might look something like this: | |
# snp chrom position major minor | |
# 1 rs9701055 1 555296 C T | |
# 2 rs9699599 1 558185 A G | |
# 3 rs12565286 1 711153 G C | |
# 4 rs28659788 1 713170 C G | |
# 5 rs11804171 1 713682 T A | |
# 6 rs2977656 1 719811 C T | |
# 7 rs12082473 1 730720 G A | |
# 8 rs12138618 1 740098 G A | |
# 9 rs3094315 1 742429 G A | |
# 1 rs3131972 1 742584 A G | |
# Select SNPs 500kb apart | |
keepRows <- pickSNPs(map,dist=500000) | |
mapsubset <- map[keepRows, ] |
To leave a comment for the author, please follow the link and comment on their blog: Getting Genetics Done.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.