Question: Obtaining Ld Values For A Snp Against All Others Within A Defined Region Using Ncbi2R / Plink
2
gravatar for Agatha
7.8 years ago by
Agatha350
Agatha350 wrote:

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.

ncbi ld snp maf R plink snps • 4.0k views
ADD COMMENTlink modified 7.8 years ago by brentp23k • written 7.8 years ago by Agatha350
2
gravatar for brentp
7.8 years ago by
brentp23k
Salt Lake City, UT
brentp23k wrote:

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.

ADD COMMENTlink written 7.8 years ago by brentp23k

Exactly what I needed. Thanks a lot.

ADD REPLYlink written 7.8 years ago by Agatha350
1
gravatar for Peixe
7.8 years ago by
Peixe600
Spain
Peixe600 wrote:

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!

ADD COMMENTlink written 7.8 years ago by Peixe600

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

ADD REPLYlink written 7.8 years ago by Peixe600

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

ADD REPLYlink modified 7.8 years ago • written 7.8 years ago by Agatha350

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

ADD REPLYlink written 7.8 years ago by zx87549.9k

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?

ADD REPLYlink written 7.8 years ago by Agatha350
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1138 users visited in the last hour