Get Rs Number Based On Chromosome and Position

Question

[aiteteji]

Hello Biostars,

I have a list of chromosome and coordinates in my vcf file, see examples:

1:10327
1:10434
1:10440
1:10469
...

Is there anyway to get the rs number for each coordinate using plink or bcftools? I found some old answers here, but I'll like to know how this is currently done.

[finswimmer]

The example you are showing is not a vcf file. So, is the example you've posted realy the data you have?

[zx8754]

Duplicate of: - Get Rs Number Based On Position

[cpad0112]

For annotating with dbsnp (or any other relevant vcf/bed), you can use bcftools annotate function or bedtools intersect function, provided variants are in vcf. OP is not in vcf.

Answer 1

[Alex Reynolds]

If you have disk space, you can set up a local copy of dbSNP for your assembly of interest, e.g. hg38:

$ VCF=ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh38p7/VCF/common_all_20180418.vcf.gz
$ wget -qO- ${VCF}
    | gunzip -c -
    | convert2bed --input=vcf --sort-tmpdir=${PWD} -
    | awk '{ print "chr"$0; }' -
    | starch --omit-signature -
    > common_20180418.starch

Then you can query your variants.vcf file against this dbSNP file any time you want:

$ bedops --element-of 1 common_20180418.starch <(vcf2bed < variants.vcf) | cut -f4 > answer.txt

If you want the BED records along with the rsId:

$ bedops --element-of 1 common_20180418.starch <(vcf2bed < variants.vcf) > answer.bed

If you have more disk space and more prep time, you can download a local copy of all SNPs:

$ VCF=ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh38p7/VCF/All_20180418.vcf.gz
$ wget -qO- ${VCF} | ...

The difference between "All" and "common" dbSNP sets is described in this document.

Answer 2

[Emily]

Run your variants through the Ensembl VEP. This will give you the known rsID if there is one, plus information about which genes were hit and the effect on those genes.

[davidgroenewoud]

Hi Emily. I am dealing with a similar issue as the OP. I tried your solution but Ensembl VEP does not give me any column with the rs identifier?

[Emily]

The column is called colocated variants. If there is not any data in that column it may be that your variants have not been observed before and do not have an rs number.

[qwang178]

I have tried vep for a vcf file intending to add rsID, in vain. Below is the command:

vep -i input.vcf -o output.vcf --vcf --cache --force_overwrite --cache_version 95 --fork 4

All the variant sites should have corresponding rsID, by double checking. Below is one of those lines taken from the output:

1       275654  .       T       C       .       PASS    CSQ=C|intron_variant&non_coding_transcript_variant|MODIFIER|AP006222.1|ENSG00000228463|Transcript|ENST00000424587|processed_transcript||2/3||||||||||-1||Clone_based_ensembl_gene|,C|intron_variant&non_coding_transcript_variant|MODIFIER|AP006222.1|ENSG00000228463|Transcript|ENST00000441866|processed_transcript||3/3||||||||||-1||Clone_based_ensembl_gene|

Do you know what might be the problem?

[zx8754]

Running online VEP gives: rs762923328, maybe check if reference files are uptodate?