Hello,

is it possible easily to get rs id's having genomic locations? Ex:
from this:

    Chr1	158669597	158669597
    Chr11	72946311	72946311

To this:

    Chr1	158669597	158669597	rs141159720
    Chr11	72946311	72946311	rs145119561

My idea was to dobnload all SNP's from UCSC or Ensembl and then compare two files by genomic coordinates in R it's simply:

    mydata$rs <- allsnps$rs[match(mydata$loc,allsnps$loc)]

However list of all human SNP's is huge, reference file has anout 20gb.
Do you know any easy solution to get result such like this?

Thanks.


If you're worried about disk usage, [BEDOPS][1] does set operations on compressed files in a format called [Starch][2], which offers a higher compression ratio than bzip2 and, as mentioned, allows direct set operations, without any extraction to an intermediate file.

With BEDOPS installed, you could create a Starch file of SNPs like so:

    $ wget -qO- http://hgdownload.cse.ucsc.edu/goldenpath/hg19/database/snp147.txt.gz \
      | gunzip -c \
      | awk -v OFS="\t" '{ print $$2,$$3,($$3+1),$$5 }' \
      | sort-bed - \
      | starch - \
      > hg19.snp147.starch

(Adjust as needed for the SNPs you need, for the reference genome you're working with.)

An uncompressed raw BED file of SNPs is roughly 5.5 GB. Compressed as a Starch file, as shown above, it should be about 5-10% of that, around 250-500 MB.

Given a sorted, tab-delimited genomic regions file `regions.bed` like this:

    chr1    158669597   158669598
    chr11   72946311    72946312

Then you could do a `bedmap` operation with the compressed SNP annotations like this:

    $ bedmap --echo --echo-map-id --delim '\t' regions.bed hg19.snp147.starch
    chr1    158669597       158669598    
    chr11   72946311        72946312    rs765534001

In this example `chr1:158669597-158669598` does not associate with a SNP ID, while `chr11:72946311-72946312` associates with SNP `rs765534001`.

So there are data hygiene issues with your input you should fix before using this toolkit:

1. Your regions should be half-open, zero-based. This means that the start and stop columns should not be equal. You can use `awk` to increment the stop field by 1, where start and stop are equal.
2. Your regions should be sorted with BEDOPS `sort-bed` for set operations to run quickly and correctly.
3. Your regions should start with `chr` and not `Chr`. This can be fixed with `awk`, as well.

  [1]: http://bedops.readthedocs.io/en/latest/
  [2]: http://bedops.readthedocs.io/en/latest/content/reference/file-management/compression/starch.html