需求
客户随手丢来一个基因型文件，类似于hapmap格式，只是少了中间多余的那几列，像这种类hapmap格式文件，往往是芯片数据。
这样的数据因为缺乏等位基因：参考碱基和变异碱基信息，对应在vcf文件中就是REF和ALT，导致后续一些分析没法进行。
那么，问题来了：怎么根据这个基因型文件来推断参考和变异等位基因？
样本量大的时候是否能通过计算等位基因频率来判断？推断出来的结果不一定准确，鬼知道你的变异多不多？
解决
在网上查了下，不能只通过基因型文件来推断，还需要依赖一个参考变异文件，有两条途径：
方法一
在ensembl中下载参考变异文件：
http://ftp.ebi.ac.uk/ensemblgenomes/pub/plants/current/variation/vcf/
但愿有你的物种吧，记得注意版本。
国内根本访问不了，我游遍世界下了半天才下下来。
以玉米为例：
这其实相当于一个单倍型的参考文件，再次强调注意版本和你的基因型文件一致。
有了这个文件就可以和基因型文件的位置相匹配，然后得到参考和变异碱基了。
示例代码：
awk 'NR==FNR{line[$1" "$2]=$5" "$6; next} ($0 in line){print $0" "line[$0]; next} {print $0, "NA"}' zea_mays.vcf pos.txt
这个代码是错误的，awk数组的值不能连接两个字段，只能等于$5，而非想要的$5" "$6。还是不熟悉，放弃，希望有高手指点下。
写了个长长的垃圾perl代码：
#! /usr/bin/perl -w
use strict;
my %hash;
my %pos;
open(IN,"<$ARGV[0]") or die $!;
while(<IN>){
chomp;
next if /^#/;
my @F = split/\s+/;
my $key = "$F[0]\t$F[1]";
my $value = "$F[3]\t$F[4]";
$hash{$key}=$value;
}
open(ID,"<$ARGV[1]") or die $!;
while(<ID>){
chomp;
my @F = split/\s+/;
my $key = "$F[0]\t$F[1]";
$pos{$key}=1;
}
foreach my $id(keys %pos){
if(exists($hash{$id})){
print "$id\t$hash{$id}\n";
}else{
print "$id\t-\t-\n";
}
}
close IN;
close ID;
最后的结果要排下序：
perl map.pl zea_mays.vcf pos.txt >out
sort -nk 1 -nk 2 out >ref_res.txt
注意，因为是参考单倍型，不一定包含了基因型文件中的所有位点。后续要怎么搞？如果缺失不多，就删了那些位点吧。
如果你的基因型文件本身是vcf格式，那用vcftools就有这种取交集位点的功能，很方便。
方法二
Ensembl 有REST API 接口，需要准备好对应的json格式文件，进行调取。
GET overlap/region/:species/:region
http://rest.ensembl.org/documentation/info/overlap_region
可能更慢更复杂些，这里不尝试了。
Ref：Question: How to get REF and ALT alleles from a genotype data?
