I am trying to request the LD information from NCBI in vicinity(+/- 100kb) of a list of snps via the NCBI2R package.
require(NCBI2R)
snp_info<-data.frame(c("rs17578796"),c("6"),c(52117256),c('T'),c('C'),c(6.730487))
colnames(snp_info)<-c("localname37","chromosome","position_37","allele1","allele2","log_p")
snp_info
localname37 chromosome position_37 allele1 allele2 log_p
1 rs17578796 6 52117256 T C 6.730487
kbs<-100^3
p_left<-pos-kbs
p_right<-pos+kbs
tryCatch( {get_left<-GetLDInfo(as.character(snp_info$chromosome),p_left,pos)},error=function(e){ print("\nErro in getting LD LEFT")
return })
tryCatch( {get_right<-GetLDInfo(as.character(snp_info$chromosome),pos,p_right)},error=function(e){ print("\nErro in getting LD RIGHT ")
return })
left<-subset(get_left, get_left$SNPB==snp_info$localname37)
right<-subset(get_right,get_right$SNPA==snp_info$localname37)
Although I am getting the LD info in the region, none of the pairs returned contain my query SNP . Is there another way to get the LD information and MAF for pairs of SNPs within a defined region? Or perhaps a pairwise LD computatuon between a SNP and a number of SNPs within a region ? I have also tried PLINK http://pngu.mgh.harvard.edu/~purcell/plink/ld.shtml
plink --file mydata
--r2
--ld-snp rs12345
--ld-window-kb 1000
--ld-window 99999
--ld-window-r2 0
But I have failed in generating the PLINK compatible file.
Any suggestions? Thank you.
2 answers
Not sure about the R code, but to use plink, I wrote a script to do this based on some code by Ryan D here on biostar:
If you want to know linkage for a single snp to all others, just use e.g.:
python linkage.py chr11:1240203-1247497 rs1234 > link.txt
otherwise, it will output all vs all for any SNP in that region.
It requires plink, vcftools and tabix.
I have tried SNAP too and it says that no information can be found for that specific snp. I have tried with the 1000 genomes and hapmap ...I have no idea why I cant get the info SNP Proxy Distance RSquared DPrime Arrays Chromosome Coordinate_HG18 rs17578796 WARNING Query snp not in 1000GenomesPilot1 rs17578796 WARNING No matching proxy snps found
I couldn't find rs17578796 in 1000G but you can use ensemble to get LD from HapMap population (CSHL-HAPMAP:HapMap-CEU [size: 180])
http://www.ensembl.org/Homo_sapiens/Share/9f935300377f2cf32a81397989e7e78b83110124
Log in to answer this question.