Obtaining Ld Values For A Snp Against All Others Within A Defined Region Using Ncbi2R / Plink

Question

[Agatha]

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.

Answer 1

[brentp]

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.

[Agatha]

Exactly what I needed. Thanks a lot.

Answer 2

[Peixe]

Try the SNAP!

It accepts a list of SNPs, and calculates all the pairwises between them in a specified region (up to 500 kb) for any HapMap/1KG Population

Hope it helps!

[Peixe]

...and in the output options you can include the MAF!

[Agatha]

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

[zx8754]

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

[Agatha]

Yes it's there..thanks. However i have checked the Biomart API and there is not way to get the information automatically, is it - for a higher number of snps?