交叉比对，取最大的，求思路或者程序。

54red

有两个输入 ：

file1：：：：：gene_list.txt
PEX19
NDUFS2

file2：：：：：NM.txt
NM_015840        ADAR        15        3603
NM_001111        ADAR        15        3681
NM_001193495        ADAR        14        2796
NM_001025107        ADAR        14        2796
NM_002857        PEX19        8        900
NR_036493        PEX19        1        2
NR_036492        PEX19        1        2
NM_001193644        PEX19        8        840
NM_004550        NDUFS2        14        1392
NM_001166159        NDUFS2        13        1374
NM_020435        GJC2        1        1320

NM中第二列（ADAR）在gene_list存在，比如PEX19有四个符合的，
NM_002857        PEX19        8        900
NR_036493        PEX19        1        2
NR_036492        PEX19        1        2
NM_001193644        PEX19        8        840
第三列取最大的，此处是8，有两个，取第四行最大的900，
print OUT NM_002857        PEX19        8        900；

我正在写程序，但是写的好复杂，请问各位大神，能提示一下么！

xiumu2280

简单的写了一下，本质应该是个排序问题，用了引用

这样的写法应该可以找到每个基因对应的数值最大条目

后续还需要你自己更改一下

1. use strict;
   
2. use warnings;
   
3. use Data::Dumper;
   
4. use File::Basename;
   
5. 
   
6. 
   
7. my %hash;
   
8. 
   
9. while (<DATA>) {
   
10.         chomp;
    
11.         my @data = split /\s+/,my $line = $_;
    
12.         if ($hash{$data[1]}) {
    
13.                 if (($data[2] == $hash{$data[1]}->[2] && $data[3] > $hash{$data[1]}->[3]) || $data[2] > $hash{$data[1]}->[2]) {
    
14.                         $hash{$data[1]} = \@data;
    
15.                 }
    
16.         }else{
    
17.                 $hash{$data[1]} = \@data;
    
18.         }
    
19.        
    
20. }
    
21. 
    
22. 
    
23. print Dumper \%hash;
    
24. 
    
25. 
    
26. __DATA__
    
27. NM_015840        ADAR        15        3603
    
28. NM_001111        ADAR        15        3681
    
29. NM_001193495        ADAR        14        2796
    
30. NM_001025107        ADAR        14        2796
    
31. NM_002857        PEX19        8        900
    
32. NR_036493        PEX19        1        2
    
33. NR_036492        PEX19        1        2
    
34. NM_001193644        PEX19        8        840
    
35. NM_004550        NDUFS2        14        1392
    
36. NM_001166159        NDUFS2        13        1374
    
37. NM_020435        GJC2        1        1320

复制代码

MMMIX

回复 1# 54red


    把你的代码贴上来看看

jason680

回复 1# 54red

>>...求思路

先读取file2，第二列当key，保存最大值，就好了...

第一笔，没有key资料(新资料),直接(最大值)保存
NM_015840        ADAR        15        3603
key: ADAR   
$gene{'ADAR'} = [( NM_015840        ADAR        15        3603 )]   # 思路,非程序代码

下一笔,先比第三列，再比第四列
NM_015840        ADAR        15      3683
key: ADAR   
$gene{'ADAR'} = [( NM_015840        ADAR        15        3683 )]   # 思路,非程序代码
保存最大值


再读取file1,依key值读取(最大值)
$gene{'PEX19'} <--- 保存最大值

54red

十分谢谢！我一直在看你的代码学习perl,可是引用这块始终没学好，再次谢谢！回复 2# xiumu2280


   

54red

谢谢jason680大神的思路，还把关键的引用怎么写都给我写好了 ^^回复 4# jason680


   

ba_du_co

1. 
   
2. my @genelist = qw/PEX19 NDUFS2/;
   
3. my %gene = map { $_, 1 } @genelist;
   
4. my %nm;
   
5. 
   
6. while (<DATA>) {
   
7.     my ( undef, $key, $v1, $v2 ) = split;
   
8.     next unless $gene{$key};
   
9. 
   
10.     if ( exists $nm{$key} ) {
    
11.         next if $v1 < $nm{$key}[0];
    
12. 
    
13.         if ( $v1 > $nm{$key}[0] or $v2 > $nm{$key}[1] ) {
    
14.             $nm{$key} = [ $v1, $v2, $_ ];
    
15.         }
    
16.     }
    
17.     else {
    
18.         $nm{$key} = [ $v1, $v2, $_ ];
    
19.     }
    
20. }
    
21. 
    
22. print $_->[2] for values %nm;
    
23. 
    
24. __DATA__
    
25. NM_015840        ADAR        15        3603
    
26. NM_001111        ADAR        15        3681
    
27. NM_001193495        ADAR        14        2796
    
28. NM_001025107        ADAR        14        2796
    
29. NM_002857        PEX19        8        900
    
30. NR_036493        PEX19        1        2
    
31. NR_036492        PEX19        1        2
    
32. NM_001193644        PEX19        8        840
    
33. NM_004550        NDUFS2        14        1392
    
34. NM_001166159        NDUFS2        13        1374
    
35. NM_020435        GJC2        1        1320
    

复制代码

54red

刚刚看见你的程序，谢谢啦！今天下班了，明天再仔细看看 ~回复 7# ba_du_co


   

sunzhiguolu

其实你的问题在于不要将解决问题的思路放在细节, 从整体看待这个问题就相对好解决一些.

RE_HASH

my @genelist = qw/PEX19 NDUFS2/;
my %gene = map { $_, 1 } @genelist;

while (<DATA>) {
    my ( undef, $key, $v1, $v2 ) = split;
    next unless $gene{$key};
    $OUT->{$key}->{$v1 * 1000000 + $v2} = $_;
}
map {
  print $OUT->{$_}->{(sort {$b<=>$a} keys %{$OUT->{$_}})[0]} ;
}sort keys %$OUT;

__DATA__
NM_015840        ADAR        15        3603
NM_001111        ADAR        15        3681
NM_001193495        ADAR        14        2796
NM_001025107        ADAR        14        2796
NM_002857        PEX19        8        900
NR_036493        PEX19        1        2
NR_036492        PEX19        1        2
NM_001193644        PEX19        8        840
NM_004550        NDUFS2        14        1392
NM_001166159        NDUFS2        13        1374
NM_020435        GJC2        1        1320


$> perl test.pl
NM_004550        NDUFS2        14        1392
NM_002857        PEX19        8        900