0

给定一组基因和现有的一对基因，我想生成尚未存在的新基因对。

基因文件具有以下格式： 

123    
134   
23455  
3242  
3423  
...  
...

基因对文件具有以下格式： 

12,345    
134,23455   
23455,343  
3242,464452  
3423,7655  
...  
...

但是我仍然在 known_interactions 和 new_pairs 之间得到了一些共同的元素。我不确定错误在哪里。

对于参数，
perl generate_random_pairs.pl entrez_genes_file known_interactions_file 250000
我得到了一个共同的元素 15880。数字 250000 是告诉我希望程序生成多少个随机对。

#! usr/bin/perl

use strict;
use warnings;

if (@ARGV != 3) {
    die "Usage: generate_random_pairs.pl <entrez_genes> <known_interactions> <number_of_interactions>\n";
}
my ($e_file, $k_file, $interactions) = @ARGV;

open (IN, $e_file) or die "Error!! Cannot open $e_file\n";
open (IN2, $k_file) or die "Error!! Cannot open $k_file\n";

my @e_file = <IN>; s/\s+\z// for @e_file;
my @k_file = <IN2>; s/\s+\z// for @k_file;

my (%known_interactions);

my %entrez_genes;
$entrez_genes{$_}++ foreach @e_file;

foreach my $line (@k_file) {
    my @array = split (/,/, $line);
    $known_interactions{$array[0]} = $array[1];
}
my $count = 0;

foreach my $key1 (keys %entrez_genes) {
    foreach my $key2 (keys %entrez_genes) {
        if ($key1 != $key2) {
            if (exists $known_interactions{$key1} && ($known_interactions{$key1} == $key2)) {next;}
            if (exists $known_interactions{$key2} && ($known_interactions{$key2} == $key1)) {next;}
            if ($key1 < $key2) { print "$key1,$key2\n"; $count++; }
            else { print "$key2,$key1\n"; $count++; }
        }
        if ($count == $interactions) {
            die "$count\n";
        }
    }
}
4

2 回答 2

0

我看不出你的代码有什么问题。我想知道您的数据中是否有一些空格 - 无论是在逗号之后还是在行尾？例如，仅提取数字字段会更安全

my @e_file = map /\d+/g, <IN>;

此外，最好将对的两个元素都保留为哈希键，这样您就可以检查元素的存在。而且，如果您确保较低的数字始终是第一个，则无需进行两次查找。

这个例子应该适合你。它没有解决您要求的随机选择部分，但这不在您自己的代码中，也不是您的直接问题

use strict;
use warnings;

@ARGV = qw/ entrez_genes.txt known_interactions.txt 9 /;

if (@ARGV != 3) {
    die "Usage: generate_random_pairs.pl <entrez_genes> <known_interactions> <number_of_interactions>\n";
}

my ($e_file, $k_file, $interactions) = @ARGV;

open my $fh, '<', $e_file or die "Error!! Cannot open $e_file: $!";
my @e_file = sort { $a <=> $b } map /\d+/g, <$fh>;

open $fh, '<', $k_file or die "Error!! Cannot open $k_file: $!";
my %known_interactions;
while (<$fh>) {
  my $pair = join ',', sort { $a <=> $b } /\d+/g;
  $known_interactions{$pair}++;
}

close $fh;

my $count = 0;
PAIR:
for my $i (0 .. $#e_file-1) {
  for my $j ($i+1 .. $#e_file) {
    my $pair = join ',', @e_file[$i, $j];
    unless ($known_interactions{$pair}) {
      print $pair, "\n";
      last PAIR if ++$count >= $interactions;
    }
  }
}

print "\nTotal of $count interactions\n";
于 2012-08-14T02:44:43.840 回答
-1

首先，您不会从已知交互的文件中咀嚼（删除换行符）。这意味着给定一个文件，如：

1111,2222

您将构建此哈希：

 $known_interactions{1111} = "2222\n";

这可能就是您收到重复条目的原因。我的猜测是（如果没有您的实际输入文件无法确定）这些循环应该可以正常工作：

 map{
    chomp;
    $entrez_genes{$_}++ ;
 }@e_file;

和

map {
    chomp;
    my @array = sort(split (/,/));
    $known_interactions{$array[0]} = $array[1];
}@k_file;

此外，作为一般规则，如果我对交互对进行排序（生物信息学的乐趣:)），我会发现我的生活会更轻松。这样我就知道 111,222 和 222,111 将以相同的方式处理，并且我可以避免在代码中使用多个 if 语句。

然后你的下一个循环将是（恕我直言更具可读性）：

my @genes=keys(%entrez_genes);
for (my $i=0; $i<=$#genes;$i++) {
   for (my $k=$n; $k<=$#genes;$k++) {
     next if $genes[$n] == $genes[$k];
     my @pp=sort($genes[$n],$genes[$k]);
     next unless exists $known_interactions{$pp[0]};
     next if $known_interactions{$pp[0]} == $pp[1];
     print "$pp[0], $pp[1]\n";
     $count++;
     die "$count\n" if $count == $interactions;
  }
}
于 2012-08-14T01:36:00.077 回答