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

最近需要下载千人基因组计划中的中国人数据作为参考,于是学习了下几个下载方法。发现还是有几个小技巧值得分享的,就记录分享一下!我找到的主要方法有两三个,分别是通过网页下载,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

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

发表评论