2013-08-26 139 views
0

我目前正在比较两个基因列表,目的是找到两个列表之间的重叠基因。哈希键的松散匹配?

此刻,我的基因的名字存储为一个哈希键两个列表(blast1和BLAST2),并找到钥匙(基因)同时存在于两个哈希值:

输入1:

XLOC_000157_6.21019:12.8196,_Change:1.04564,_p:0.04915,_q:0.999592  99.66 gi|475392713|dbj|AB759708.1|_Xenopus_laevis_PhyHd_mRNA_for_phytanoyl-CoA_dioxygenase_like_protein,_complete_cds 
XLOC_000159_636.025:343.104,_Change:-0.890436,_p:0.00575,_q:0.999592 99.47 gi|9909981|emb|AJ278067.1|_Xenopus_laevis_mRNA_for_putative_XIRG_protein 
XLOC_000561_31.1018:14.9273,_Change:-1.05905,_p:0.0073,_q:0.999592  91.57 gi|165973401|ref|NM_001113689.1|_Xenopus_(Silurana)_tropicalis_cytokine_inducible_SH2-containing_protein_(cish),_mRNA 

指定用于第一个基因列表...

$input1 = $ARGV[0]; 
open my $blast1, '<', $input1 or die $!; 

my $results1 = 0; 
my (@blast1ID, @blast1_info, @percent_id, @split); 
while (<$blast1>) { 
    chomp; 
    @split = split('\t'); 
    push @blast1_info, $split[0]; 
    push @percent_id, $split[1]; 
    push @blast1ID, $split[2]; 
    $results1++; 
} 

print "$results1 blast hits in '$input1'\n"; 

push @{$blast1{$blast1ID[$_]} }, [ $blast1_info[$_], $percent_id[$_] ] for 0 .. $#blast1ID; 

输入2:

XLOC_000561_31.1018:14.9273,_Change:-1.05905,_p:0.0073,_q:0.999592  91.57 gi|165973401|ref|NM_001113689.1|_Xenopus_(Silurana)_tropicalis_cytokine_inducible_SH2-containing_protein_(cish),_mRNA 
XLOC_000679_57.3461:29.2637,_Change:-0.970585,_p:0.03645,_q:0.999592 85.13 gi|51704135|gb|BC081195.1|_Xenopus_laevis_hypothetical_protein_LOC446937,_mRNA_(cDNA_clone_IMAGE:6640116),_partial_cds 
XLOC_000766_10.699:6.33756,_Change:-0.755473,_p:0.0384,_q:0.999592  99.04 gi|195972824|ref|NM_001130940.1|_Xenopus_laevis_interleukin_6_signal_transducer_(gp130,_oncostatin_M_receptor)_(il6st),_mRNA 

指定为第2个基因的列表

$input2 = $ARGV[1]; 
open my $blast2, '<', $input2 or die $!; 

my $results2 = 0; 
my (@blast2ID, @blast2_info, @percent_id); 
while (<$blast2>) { 
    chomp; 
    @split = split('\t'); 
    push @blast2_info, $split[0]; 
    push @percent_id, $split[1]; 
    push @blast2ID, $split[2]; 
    $results2++; 
} 
print "$results2 blast hits in '$input2'\n"; 

push @{$blast2{$blast2ID[$_]} }, [ $blast2_info[$_], $percent_id[$_] ] for 0 .. $#blast2ID; 

查找键同时存在于两个哈希值(基因):

my $intersect_count = 0; 
for my $key (sort keys %blast1) { 
    if (exists $blast1{$key} && $blast2{$key}) { 
     $intersect_count++; 
      for my $part1 (@ { $blast1{$key} }) { 
       ($hit1, $percent_id1) = @$part1; 
      } 
      for my $part2 (@ { $blast2{$key} }) { 
       ($hit2, $percent_id2) = @$part2; 
      } 
    push @intersect, "$key\tC1:$hit1 [$percent_id1]\tC2:$hit2 [$percent_id2]\n";    
    push @intersecting_list, "$key";     
    } 
} 

上面的代码会发现一个基因是存在于两个列表:

gi|165973401|ref|NM_001113689.1|_Xenopus_(Silurana)_tropicalis_cytokine_inducible_SH2-containing_protein_(cish),_mRNA 

我的问题是我该如何适应这个基因,输出中包含相似的名称?比如我希望看到:

gi|186928837|ref|NM_005982.3|_Homo_sapiens_SIX_homeobox_1_(SIX1),_mRNA 

找到一个匹配:

gi|154142326|ref|NM_001100275.1|_Xenopus_(Silurana)_tropicalis_SIX_homeobox_1_(six1),_mRNA 

有什么建议?

+0

'NM_005982.3' * like *'NM_001100275.1'如何?你想只匹配最初的两个字符吗? – Borodin

+0

不 - 这是一件困难的事情 - 我试图在基因名称上匹配 - 例如'SIX_homeobox_1_(SIX1)'和'SIX_homeobox_1_(six1)'。 – fugu

+1

那么,如果你想要独立匹配基因名称的* end *,那就没问题了。有没有像最后一个* 4 *字段(用下划线分隔),总是必须匹配的部分?你必须以某种方式定义“相似”。 – Borodin

回答

2

有两种策略可以使用

  1. 提取要使用实际的键,然后是精确匹配。

    原始密钥的某些部分可能对您没有任何用处 - 删除它们。根据输入,您可能还想要进行Unicode规范化,并执行大小写折叠。

    在你的情况下,

    gi|186928837|ref|NM_005982.3|_Homo_sapiens_SIX_homeobox_1_(SIX1),_mRNA 
    gi|154142326|ref|NM_001100275.1|_Xenopus_(Silurana)_tropicalis_SIX_homeobox_1_(six1),_mRNA 
    

    公用密钥可能看起来像

    gi|ref|nm_00|_six_homeobox_1_(six1),_mrna 
    
  2. 破除哈希值,并计算所有可能的记录之间的相似性指数。想知道这样的指数,你可能想看看Levenstein edit distance。然后,您可以将特定边界内的所有其他记录视为匹配。这相当昂贵,但可能会产生更好的结果。

我不知道你的问题域,所以我不能提出任何好的建议。


您的代码存在一些问题,特别是在查找匹配时。它看起来像它应该是相同的:

my $intersect_count = 0; 
for my $key (sort keys %blast1) { 
    if (exists $blast2{$key}) { 
     $intersect_count++; 
     my ($hit1, $percent_id1) = @{ $blast1{$key}[-1] }; 
     my ($hit2, $percent_id2) = @{ $blast2{$key}[-1] }; 
     push @intersect, "$key\tC1:$hit1 [$percent_id1]\tC2:$hit2 [$percent_id2]\n"; 
     push @intersecting_list, $key; 
    } 
} 

差异:

  1. exists $blast1{$key} && $blast2{$key}被解析为exists($blast1{$key}) && $blast2{$key},甚至是愚蠢的,因为我们知道,$blast1{$key}存在:我们刚刚通过keys获取吧!
  2. 在循环数组并将每个项目分配给变量时,变量将保留最后一项的值。即my $y; for my $x (@xs) { $y = $x }相当于my $y = $xs[-1],但效率较低。
+0

这是一个梦幻般的答案,我担心的是一个不可能的问题。我一定会考虑Levenstein的编辑距离。您是否介意详细阐述您在关键查找中所做的修改? – fugu

+1

@FlyingFrog我编辑了我的答案,并对我的修改进行了一些解释。我不建议你实际使用莱文斯坦距离。这对平滑拼写错误(或点突变)等是很有用的。您将不得不设计自己的度量标准,但Levenstein可能是其中的一部分。第一步是当两个键比较相等时为自己明确定义。 – amon