如何反向推断基因型文件中的参考碱基（REF/ALT）？

目录

• 需求
• 解决
  • 方法一
  • 方法二

需求

客户随手丢来一个基因型文件，类似于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]	$F[1]";
    my $value = "$F[3]	$F[4]";
    $hash{$key}=$value;
}

open(ID,"<$ARGV[1]") or die $!; 
while(<ID>){
    chomp;
    my @F = split/s+/;
    my $key = "$F[0]	$F[1]";
    $pos{$key}=1;
}

foreach my $id(keys %pos){
    if(exists($hash{$id})){
        print "$id	$hash{$id}
";    
    }else{
        print "$id	-	-
";    
    }   
}
    
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

提供小麦的一个示例：
http://rest.ensembl.org/overlap/region/triticum_aestivum/4A:714193714-714193714?content-type=application/json;feature=variation

可能更慢更复杂些，这里不尝试了。

Ref：Question: How to get REF and ALT alleles from a genotype data?