如果只有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

參考

1. https://gist.github.com/martijnvermaat/09cfa3ec1aeaca9d6dec
2. https://www.biostars.org/p/349284/