将SNP坐标映射到基因坐标太慢了。

Question

我有两个由制表符分隔的文件，如下所示

File 1 (这些是单核苷酸多态性(SNP)位置)

Chr1 26690
Chr1 33667
Chr1 75049
.
.

Chr2 12342
Chr2 32642
Chr2 424421
.
.

文件2 (这些是基因起始坐标和结束坐标)

Chr1    2903    10817   LOC_Os01g01010
Chr1    2984    10562   LOC_Os01g01010
Chr1    11218   12435   LOC_Os01g01019
Chr1    12648   15915   LOC_Os01g01030
Chr1    16292   18304   LOC_Os01g01040
Chr1    16292   20323   LOC_Os01g01040
Chr1    16321   20323   LOC_Os01g01040
Chr1    16321   20323   LOC_Os01g01040
Chr1    22841   26971   LOC_Os01g01050
Chr1    22841   26971   LOC_Os01g01050
.
.

我想要的是将文件1中的SNP与文件2中的基因相匹配。脚本应该匹配文件第一列中的字符串，如果它们匹配，则应该找到文件2中的哪个基因包含相应的SNP，并从文件2的第四列返回轨迹ID。

这是我写的剧本

use strict;

my $i1 = $ARGV[0];            # SNP 
my $i2 = $ARGV[1];            # gene coordinate

open(I1, $i1);
open(I2, $i2);

my @snp  = ();
my @coor = ();

while( <I1> ) {
  push(@snp, $_);    
}

while ( <I2> ) {
  push(@coor, $_);   
}

for ( my $i = 0; $i <= $#snp; $i++ ) {

    my @snp_line = split "\t", $snp[$i];

    for ( my $j = 0; $j <= $#coor; $j++ ) {

        my @coor_line = split "\t", $coor[$i];   

        if ( $snp_line[0] eq $coor_line[0] ) {

            if ( $snp_line[1] >= $coor_line[1] && $snp_line[1] <= $coor_line[2] ) {
              print "$snp_line[0]\t$snp_line[1]\t$coor_line[3]\n";
              goto a;
            }
        }   
    }
a: 
}

问题是，这显然不是最好的方法，因为它对第1行中的每个SNP遍历了文件2中的所有60,000行。而且，它在一夜之间运行，没有通过Chr1；我们已经到了Chr12。

Answer 1

当将这些文件重新格式化为UCSC床格式时，您可以使用像BEDOPS这样的工具包对排序的床文件进行高效的设置操作。

将SNP的第一个文件转换为已排序的床文件：

$ awk -v OFS="\t" '{ print $1, $2, ($2+1); }' snps.txt | sort-bed - > snps.bed

对基因进行排序(“文件2")：

$ sort-bed genes.unsorted.txt > genes.bed

将SNP映射到基因：

$ bedmap --echo --echo-map-id-uniq --delim '\t' snps.bed genes.bed > answer.bed

如果需要，可以从答案中去掉SNP的末端位置：

$ cut -f1,2,4 answer.bed > answer.txt

这些工具将运行非常快，通常在几分钟内。

我不会使用Perl或Python来执行这类set操作，除非我正在做某种学术练习。

Answer 2

下面是一个工作脚本，上面发布的脚本有错误

use strict;

my $i1=$ARGV[0];            # SNP 
my $i2=$ARGV[1];            # gene coordinate

open(I1,$i1);
open(I2,$i2);

my @snp=();
my @coor=();

while(<I1>)
{
push(@snp,$_);  
}
while(<I2>)
{
push(@coor,$_); 
}

for(my $i=0;$i<=$#snp;$i++)
{
my @snp_line = split "\t",$snp[$i];

for(my $j=0;$j<=$#coor;$j++)
{
    my @coor_line = split "\t",$coor[$j];

    if ($snp_line[0] eq $coor_line[0])
    {   
        if ($snp_line[1] >= $coor_line[1] && $snp_line[1] <= $coor_line[2])
        {
            print "$snp_line[0]\t$snp_line[1]\t$coor_line[3]\n";                
        }

    }   

}

}

这个人做这件事。