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

分配第二个基因列表

$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

有什么建议么？

4

1 回答 1

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 编辑距离。然后，您可以将特定范围内的所有其他记录视为匹配项。这是相当昂贵的，但可能会产生更好的结果。

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

您的代码存在一些问题，尤其是在查找命中时。看起来它应该等同于：

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]。
于 2013-08-26T13:36:00.007 回答