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.
•
written
by
aiteteji •
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.
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.
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?
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?
•
written
by
qwang178 •
Log in to answer this question.
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.
Duplicate of: - Get Rs Number Based On Position
The example you are showing is not a vcf file. So, is the example you've posted realy the data you have?