如果只有SNP的染色體和物理位置信息，該如何批量轉換得到 rs ID？
思路非常簡單，只需要下載 dbSNP 的參考文件，根據位置信息從參考文件中獲取對應的 rs 編號即可。
下面列舉兩個例子。
重命名 PLINK 文件 SNP 名字
第一個例子是 PLINK 格式的文件，要把 .bim 文件中的 SNP 名字改為 rs id。
首先從 UCSC 下載純文本格式的 dbSNP release 151 並解壓，這里下載的是 hg19 版本：
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/snp151.txt.gz
gzip snp151.txt.gz -d
snp151.txt.gz 包含了所有 SNPs，總共有 12G 大。如果只需要常見 SNPs，或者說硬盤不夠大，則可以下載只有常見 SNPs 的文件，只有 748M：
注意，下載的文件的 chromStart 列是 0-based 的。0-based 指的是染色體坐標從 0 開始，第一個位置記為 0，而 1-based 則是從 1 開始算出。如果還是不明白 0-based 是什么意思，可以看看 Chromosome coordinate systems: 0-based, 1-based。
接下來，新建一個 Python 腳本，腳本名字為 rename_snps_shiyanhe.py：
# 歡迎訪問實驗盒www.shiyanhe.com
import sys
snps = {'snp_%s_%s' % (e[0][3:], e[1]): e[2] for e in (l.strip().split() for l in open(sys.argv[1]))}
bim = (l.strip().split() for l in open(sys.argv[2]))
new = open(sys.argv[3], 'w')
for e in bim:
e[1] = snps.get(e[1], e[1])
new.write('\t'.join(e) + '\n')
new.close()
bim.close()
按 rename_snps_shiyanhe.py SNP參考文件 原來的bim 新的bim文件 這樣的命令執行腳本，比如：
python rename_snps_shiyanhe.py snp151.txt original.bim new.bim
根據 VCF 文件查詢 rs id
如果要查詢 vcf 文件的 SNPs，根據前面下載的參考文件，自己寫腳本即可。不過，如果不會寫腳本，也可以用下面介紹的用 BEDOPS 的方法。
根據基因組坐標版本下載 dbSNP 數據，這里下載 hg38 版本:
vcf = ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh38p7/VCF/All_20180418.vcf.gz
wget -qO- ${VCF}
| gunzip -c -
| convert2bed --input=vcf --sort-tmpdir=${PWD} -
| awk '{ print "chr"$0; }' -
| starch --omit-signature -
> All_20180418.starch
如果硬盤空間不夠，也可下載常見 SNPs 的參考文件，路徑為 ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh38p7/VCF/common_all_20180418.vcf.gz 。
假設要重命名的 vcf 文件為 shiyanhe.com.vcf ，可以這么查詢：
bedops --element-of 1 All_20180418.starch <(vcf2bed < shiyanhe.com.vcf) | cut -f4 > rsid.txt
如果想得到 bed 文件：
bedops --element-of 1 all_20180418.starch <(vcf2bed < shiyanhe.com.vcf) > shiyanhe.com.recoded.bed
參考
