Get Rs Number Based On Chromosome and Position
2
1
Entering edit mode
3.5 years ago
aiteteji ▴ 10

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.

RNA-Seq snp vcf • 7.5k views
ADD COMMENT
2
Entering edit mode

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

ADD REPLY
1
Entering edit mode
ADD REPLY
0
Entering edit mode

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.

ADD REPLY
5
Entering edit mode
3.5 years ago

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.

ADD COMMENT
2
Entering edit mode
3.5 years ago
Emily 23k

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.

ADD COMMENT
0
Entering edit mode

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?

ADD REPLY
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

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?

ADD REPLY
0
Entering edit mode

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

ADD REPLY

Login before adding your answer.

Traffic: 593 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6