千人基因组计划基因分型数据下载

最近需要下载千人基因组计划中的中国人数据作为参考，于是学习了下几个下载方法。发现还是有几个小技巧值得分享的，就记录分享一下！我找到的主要方法有两三个，分别是通过网页下载，API批量下载和ftp批量下载。

1.网页下载

这个比较简单，也已经有很多教程，列在这里，我就不写了。

千人基因组计划数据库下载某段区域SNP

传说中的千人基因组计划

2.ftp或者ascp批量下载

1)下载所有的数据和位点

如果你需要所有的数据和位点，那么ftp或者ascp批量下载是一个不错的选择，特别是后者，可以用完你的带宽，节约你的时间。

比如，这个教程提供的下载命令：【直播】我的基因组55：简单的PCA分析千人基因组的人群分布

#这个数据是分染色体下载的，还要合并，应该是几十G的数据
nohup wget  -c -r -nd -np -k -L -p  ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502 &

在1000 Genome Project找到一个ascp下载的命令：

ascp -i bin/aspera/etc/asperaweb_id_dsa.openssh -Tr -Q -l 100M -P33001 -L- fasp-g1k@fasp.1000genomes.ebi.ac.uk:vol1/ftp/release/20100804/ALL.2of4intersection.20100804.genotypes.vcf.gz ./

上面这个是linux系统下的命令，如果是Windows或者Mac参见我整理的命令，见下一篇推文。这个下载的数据看日期可能不如上面一个全面，60多G的数据，在100M的速度下，估计下载时间为两个小时左右。

2)下载特定个体或者snp

当你需要的是特定个体，区域或者位点的时候，官网的这个命令就派上了用场。

#首先，要下载个工具，从下面这个网页下载
https://codeload.github.com/vcftools/vcftools/legacy.zip/master
#然后解压
unzip vcftools-vcftools-v0.1.16-16-g954e607.zip
#安装
cd vcftools-vcftools-v0.1.16-16-g954e607
export PERL5LIB=/Users/zd200572/Downloads/vcftools-vcftools-954e607/src/perl/
./autogen.sh
#如果提示缺少autoconf安装之(ubuntu)
sudo apt install autoconf automake libtool
./configure
make
#在这里 -c 可以指定多个个体，逗号分隔，还可以是一个文件
#-c, --columns            File or comma-separated list of columns to keep in the vcf file. If file, one column per row
tabix -h ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20100804/ALL.2of4intersection.20100804.genotypes.vcf.gz 17:1471000-1472000 | perl vcftools-vcftools-954e607⁩/src⁩/⁨perl⁩/vcf-subset -c HG00098 | bgzip -c /tmp/HG00098.20100804.genotypes.vcf.gz

比较遗憾，总是报错，下载始终没有成功。报错如下：

[E::hts_open_format] Failed to open file ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20100804//ALL.2of4intersection.20100804.genotypes.vcf.gz
Could not read ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20100804//ALL.2of4intersection.20100804.genotypes.vcf.gz

可是文件明显是存在的，网页可访问，只好暂时作罢。

参考的官方网页地这：https://www.internationalgenome.org/faq/how-do-i-get-sub-section-vcf-file/

和https://www.internationalgenome.org/faq/can-i-get-haplotype-data-1000-genomes-individuals/

3.API批量下载

其实这个千人基因组计划的数据库是完全开放的，可以以匿名用户登陆的，如果你对数据库完全熟悉，完全可以以命令行或者图形客户端来查看导出，更为高效，无奈我水平不够，只有使用官方现成Perl API了。

Server  User    Password    Port
mysql-db.1000genomes.org    anonymous   -   4272
#https://www.internationalgenome.org/public-ensembl-mysql-instance/

其实作为一个API，它的安装还是极其新手不友好的，各种安装，编译。我尝试了下，没有成功，好在官方提供了虚拟机镜像，我想对于一般人，虚拟机够用了，直接下载，导入就可以使用，有图形界面，好操作。虚拟机的使用，官方给出了详细教程，基本上已经是一步步的教程，很详细，没必要翻译了，地址在这：

http://asia.ensembl.org/info/data/virtual_machine.html

#1.从这个地址下载虚拟机软件
#https://www.virtualbox.org/wiki/Downloads
#2.下载虚拟机镜像，用linux下axel 这个工具下载快点，其他多线程下载也可以，速度大概360k/s左右，下载完成要两个多小时
axel http://ftp.ensembl.org/pub/current_virtual_machine/Ensembl99.ova
#3.按上面的地址教程导入就好了

使用API测试一下

单个数据的下载

我要下载的是千人基因组某个snp的全部样本分型结果，官方有现成示例脚本，snp不多的话就直接下载了，多的话可以构建个循环进行，先看单个的下载。只做了一个改动，把$variation_adaptor->db->use_vcf(1);注释去掉了。大概5分钟一个结果，可能是数据库庞大，需要这么长时间，还是防止并发过多。反正速度不够快，不过不急的话也够用了。

perl get_genotype.pl rs123456 > rs123456.txt

#get_genotype.pl
use strict;
use warnings;
use Bio::EnsEMBL::Registry;

my $registry = 'Bio::EnsEMBL::Registry';

$registry->load_registry_from_db(
    -host => 'ensembldb.ensembl.org',
    -user => 'anonymous'
);

my $variation_adaptor = $registry->get_adaptor('homo_sapiens', 'variation', 'variation');

# OPTIONAL: uncomment the following line to retrieve 1000 genomes phase 3 data also
$variation_adaptor->db->use_vcf(1);

my $variation = $variation_adaptor->fetch_by_name($ARGV[0]);

my $genotypes = $variation->get_all_SampleGenotypes();

foreach my $genotype (@{$genotypes}) {
  print $genotype->sample->name, "\t", $genotype->genotype_string, "\n";
}

批量下载

这里我写了个python小脚本，perl曾经试图学过，没学会。很简单，就两句，应该shell更高效，估计一行就好，功力不到呀。

import os

li = [ "rs9469031", "rs6141383", "rs4809957", "rs41309931", "rs36600", "rs17728461" ]

for a in li:
    cmd = "perl ../test.pl %s > %s.txt" % (a, a)
    os.system(cmd)
    #break

好了，几番折腾后，终于获得了自己需要的数据，同样的还可以获得频率什么的。