哈哈，这里快成为我吐槽的后花园了。

在帮老公做个工作，查找一个SNP及上下游25bp内所有SNP在Asian中的频率，他一直在dbSNP网页上一个一个找，但是SNP太多，花费时间过长，于是想让我写个脚本帮他实现。

ftp上下载的population的信息太多太杂（https://ftp.ncbi.nih.gov/snp/），不好筛选，他就想得到页面上的不同数据库的Population为Asian的信息：https://www.ncbi.nlm.nih.gov/snp/rs774321397#frequency_tab，好在这个信息可以download，文件下载链接可以按照rs号批量实现（https://www.ncbi.nlm.nih.gov/snp/rs774321397/download/frequency）。

这个工作主要分为三步：

1）下载最新版本154的vcf文件（https://ftp.ncbi.nih.gov/snp/latest_release/VCF/GCF_000001405.38.gz）根据rs号提取染色体位置，万万没想到居然在这一步花费的时间最多，主要是不想使用grep，因为速度比较慢，想找一下别的方法，结果一直报错，最后还是用的grep。。。

zgrep -m 1 -P "rs774321397\t" GCF_000001405.38.gz

-m Grep匹配到首条记录后停止

-P 使用perl的正则表达式语法

### 一开始使用的命令行是zgrep -m 1  'rs774321397'$'\t' GCF_000001405.38.gz，但是这个命令行在perl里不识别。

2）根据染色体位置提取上下游25bp内所有SNP，因为vcf文件有索引文件，所以直接用tabix提取就可以了。

https://www.jianshu.com/p/b6f885fc8b5f

tabix xx.vcf.gz 11:2343540-2343596

3）根据rs号提取frequency文件，读取Population为Asian的信息，本来想用python抓取页面信息，但是我的python太菜了，还是用wget把文件下载下来，然后用perl读取文件信息。

整个过程都用perl脚本自动实现了，开心。

#################################################################

根据rs号提取染色体位置所做的尝试：

https://www.biostars.org/p/59249/

首先尝试用代码连外网数据库genome-mysql.cse.ucsc.edu，但是报连接失败，ping是可以ping通的，查了半天原因没解决。

再尝试cruzdb的方法，报缺包sqlalchemy（No module named sqlalchemy），安装成功后，但是relationship不能使用（sqlalchemy.orm ImportError: cannot import name relationship），猜可能是下载的sqlalchemy的版本太新了，没有relationship，又失败了。。。

sudo apt-cache search sqlalchemy

sudo apt-get install python-sqlalchemy

第三次尝试rsidx，tabix是根据位置建立的索引，rsidx是根据rs建立的索引，通过rs得到position信息，但是我的python版本过低，没有安装成功。

https://bioinformatics.stackexchange.com/questions/8593/index-vcf-file-for-rapid-access-by-rsid