使用tabix快速查询rsID和染色体位置

在和遗传变异打交道的过程中，我们经常需要根据基因组坐标和allele来注释rsid，通常如果只有一两个的话我们可以直接查询相关的网站，如NCBI dbSNP或者myvariant，但是如果想批量注释上百万甚至上千万的rsid，甚至有些时候需要通过rsid查询位置和allele，这个时候我们可能会想到VEP或者annovar等注释软件，但是这些软件通常需要特定的输入格式，如果我们的数据不是VCF等常规格式就比较麻烦。

这里我介绍一种可以快速查询rsid和基因组坐标的方法，原理就是利用tabix快速查询的特性，自行构建完整的tabix查询文件就可以非常快速地查询rsid和基因组坐标，以下是具体操作步骤：

1. 下载原始文件

   一共要下载两个文件，以hg19 b156为例

   1. 记录位置和rsid信息的GCF文件，格式就是没有基因型的VCF，
   2. 记录之前版本被合并过的rsid的json文件

   如果只需要通过位置注释rsid的话只需要下载第一个文件就可以了。如果需要下载hg38的话可以下载[GCF_000001405.40.gz](https://ftp.ncbi.nlm.nih.gov/snp/archive/b156/VCF/GCF_000001405.40.gz) 文件``

2. 分割multi-allele SNP

   有很多SNP有多个ALT allele，为了后续方便比对allele进行注释，我们需要使用bcftools将有多个ALT allele的SNP拆分成多行。

   命令很简单：

   bcftools norm GCF_000001405.25.gz -m - -O z -o GCF_000001405.25.nomultiallelic.gz

3. 构建pos2rsid和rsid2pos文件

   接着我们需要从GCF文件再生成两个文件，一个是以位置为索引查询rsid的文件，一个是以rsid为索引查询位置和allele等信息的文件，查询实现都是使用tabix，以位置为索引很常见很好理解，但是以rsid为索引就不常见了，起始就是把每个rsid的第一位数字当成“染色体”，整个rsid的数字当作“坐标”，这样可以利用tabix的速度优势快速查询rsid。

   具体操作的时候还要把GCF里面的染色体编码替换成常用的1-22，X，Y

   我使用python完成这个操作：

   import json
   import bz2
   import numpy as np
   import pandas as pd
   import gzip
   from subprocess import call,check_output
   
   version = 'b156'
   
   # 先统计一下GCF文件的文件头有多少行
   n_header = 0
   with gzip.open('./GCF_000001405.25.nomultiallelic.gz','rb') as f:
       for line in f:
           if line.decode('utf-8').startswith('##'):
               n_header += 1
           else:
               break
   
   chr_name_map = pd.Series(data=list(range(1, 23)) + ['X', 'Y'],
                            index=check_output('tabix -l ./GCF_000001405.25.gz',
                                               shell=True).decode().split()[:24])
   
   for df in pd.read_csv('./GCF_000001405.25.nomultiallelic.gz',
                         sep='\t',
                         skiprows=n_header,
                         compression='gzip',
                         usecols=['#CHROM', 'POS', 'ID', 'REF', 'ALT'],
                         chunksize=100000):
       # replace chr
       df['#CHROM'] = df['#CHROM'].map(chr_name_map)
   
       # remove chrM
       df = df.dropna()
       # stop when df are all chrM
       if len(df) == 0:
           break
   
       # use first int of rsid as fake chr
       df['rsid'] = df['ID'].map(lambda rsid: rsid[2:])
       df['rsid_1st'] = df['ID'].map(lambda rsid: rsid[2])
   
       # write pos2snp file
       df[['#CHROM', 'POS', 'ID', 'REF',
           'ALT']].to_csv(f'./pos2snp_{version}.txt',
                          sep='\t',
                          index=False,
                          header=False,
                          mode='a')
       # write snp2pos file
       df[['rsid_1st', 'rsid', '#CHROM', 'POS', 'REF',
           'ALT']].to_csv(f'./snp2pos_{version}_unsorted.txt',
                          sep='\t',
                          index=False,
                          header=False,
                          mode='a')

4. 构建查询merged snp的文件

   需要解析json文件到表格文件，然后排序压缩再建索引。

   with open(f'./merged_{version}.txt', 'w') as f_out:
       with bz2.BZ2File('./refsnp-merged.json.bz2', 'rb') as f_in:
           for line in f_in:
               rs_obj = json.loads(line.decode('utf-8'))
               refsnp_id = rs_obj['refsnp_id']
               for merge_into in rs_obj['merged_snapshot_data']['merged_into']:
                   f_out.write(f"{refsnp_id[0]}\t{refsnp_id}\t{merge_into}\n")
                   for dbsnp1_merges in rs_obj['dbsnp1_merges']:
                       dbsnp1_merges = dbsnp1_merges['merged_rsid']
                       f_out.write(f"{dbsnp1_merges[0]}\t{dbsnp1_merges}\t{merge_into}\n")
   parsed_df = pd.read_csv(f'./merged_{version}.txt', sep='\t', names=['1', '2', '3'])
   parsed_df = parsed_df.sort_values(['1', '2'])
   parsed_df.to_csv(f'./merged_{version}.txt', sep='\t', index=False, header=False)

5. 排序，压缩，索引

   最后就是用bgzip压缩这两个文件，然后分别建索引

   snp2pos文件压缩之前还要先排序。

   # sort snp2pos file
   call(f'sort -k1,1 -k2,2n ./snp2pos_{version}_unsorted.txt > ./snp2pos_{version}.txt', shell=True)
   
   # bgzip and tabix
   call(f'bgzip ./pos2snp_{version}.txt', shell=True)
   call(f'tabix -s 1 -b 2 -e 2 ./pos2snp_{version}.txt.gz', shell=True)
   call(f'bgzip ./snp2pos_{version}.txt', shell=True)
   call(f'tabix -s 1 -b 2 -e 2 ./snp2pos_{version}.txt.gz', shell=True)
   
   call(f'bgzip ./merged_{version}.txt', shell=True)
   call(f'tabix -C -s 1 -b 2 -e 2 ./merged_{version}.txt.gz', shell=True)

6. 测试

   然后就可以用位置查询rsid或者用rsid查询位置了。

   # 通过位置查询rsid
   $ tabix pos2snp_b156.txt.gz 1:10019-10019
   1       10019   rs775809821     TA      T
   
   # 通过rsid查询位置, 例如rs
   $ tabix snp2pos_b156.txt.gz 7:775809821-775809821
   7       775809821       1       10019   TA      T
   
   # 有些rsid被合并到新的rsid里了，所以在查询rsid的时候需要先查询它是否被合并了。
   # 例如rs1086，直接查snp2pos没有结果。
   $ tabix snp2pos_b156.txt.gz 1:1086-1086
   
   # 查询merged文件就发现它被合并到了rs940
   $ tabix merged_b156.txt.gz 1:1086-1086
   1       1086    940
   
   # 最后查询snp2pos文件得到它的坐标是13:46627480
   $ tabix snp2pos_b156.txt.gz 9:940-940
   9       940     13      46627480        C       G

7. 删除中间文件

   # remove intermediate files
   call('rm ./GCF_000001405.25.gz', shell=True)
   call('rm ./GCF_000001405.25.gz.tbi', shell=True)
   call('rm ./GCF_000001405.25.nomultiallelic.gz', shell=True)
   call('rm ./refsnp-merged.json.bz2', shell=True)
   call(f'rm ./snp2pos_{version}_unsorted.txt', shell=True)