求一个聚类的perl程序

tanger_009

我有两个文件，第一个文件(posi_infor.txt)是基因的位置信息。1到15是基因的顺序编号，每个数字后面的是基因的名称。
第二个文件(pearson_cor.txt)是pearson检验相关的基因对。比如ENSG00000186092        ENSG00000235249，说明这两个基因具有相关性。

下面开始执行一种算法，叫做动态相关性聚类算法：
根据基因的位置信息(posi_infor.txt)，从基因1依次去比较后面的基因(即基因2-15)，如果基因间的两两比较存在于pearson相关对立面，就把这些基因聚类，归于同一个类别。但是这里面有个条件，同一个类别里面的基因编号顺序之间的gap不能大于3. 举个例子，假如基因1与基因2,3,4，9存在pearson相关对。那么基因1,2,3,4可以归为同一类别，但是不能把基因9归为此类，因为基因9与基因4之间的gap是4（即四个基因，基因5,6,7,8）。

还请大仙帮忙。我想了两天了。脑子太笨，实在搞不出来了。
谢谢啊

xiumu2280

只用1 比较 2-15
不用2 比较 3-15嘛？

xiumu2280

1. #!/usr/bin/perl
   
2. use strict;
   
3. use warnings;
   
4. use Data::Dumper;
   
5. 
   
6. 
   
7. my %hash_gene = map {chomp;(split)[1,0]}<DATA>;
   
8. open (IN,"C:\\Users\\yy\\Desktop\\pearson_cor.txt");
   
9. my %hash_new;
   
10. while (<IN>) {
    
11.         chomp;
    
12.         my @data = split;
    
13.         map {my $key1 = $_;map {push @{$hash_new{$key1}},$_."_".$hash_gene{$_}}@data}@data if @hash_gene{@data};
    
14. }
    
15. ############
    
16. foreach my $key (keys %hash_new) {
    
17.         my %hash_uniq;
    
18.         @hash_uniq{@{$hash_new{$key}}}=1;
    
19.         my @data = sort {$a->[1] <=> $b->[1]}map {[(split /_/)[0],(split /_/)[1]]}keys %hash_uniq;
    
20.         print "~~~~~~\n";
    
21.         print "$data[0][0]\t$data[0][1]\n";
    
22.         map {@data[$_,$_+1] && ($data[$_+1][1]-$data[$_][1]) <= 3 && print "$data[$_+1][0]\t$data[$_+1][1]\n"}(0..$#data);
    
23. }
    
24. 
    
25. 
    
26. __DATA__
    
27. 1   ENSG00000186891
    
28. 2   ENSG00000186092
    
29. 3   ENSG00000235249
    
30. 4   ENSG00000187634
    
31. 5   ENSG00000268179
    
32. 6   ENSG00000188976
    
33. 7   ENSG00000187961
    
34. 8   ENSG00000187583
    
35. 9   ENSG00000187642
    
36. 10  ENSG00000188290
    
37. 11  ENSG00000187608
    
38. 12  ENSG00000188157
    
39. 13  ENSG00000237330
    
40. 14  ENSG00000131591
    
41. 15  ENSG00000162571
    

复制代码

1. ~~~~~~
   
2. ENSG00000186092 2
   
3. ENSG00000268179 5
   
4. ~~~~~~
   
5. ENSG00000187961 7
   
6. ENSG00000188290 10
   
7. ~~~~~~
   
8. ENSG00000186092 2
   
9. ENSG00000235249 3
   
10. ENSG00000187634 4
    
11. ENSG00000268179 5
    
12. ~~~~~~
    
13. ENSG00000187961 7
    
14. ENSG00000187642 9
    
15. ENSG00000188290 10
    
16. ~~~~~~
    
17. ENSG00000186092 2
    
18. ENSG00000187634 4
    
19. ~~~~~~
    
20. ENSG00000186092 2
    
21. ENSG00000235249 3
    
22. ~~~~~~
    
23. ENSG00000187961 7
    
24. ENSG00000187642 9

复制代码

jzp520520

回复 3# xiumu2280


    大神!这是你自己写的？还是流程里面的？膜拜！{:2_171:}

xiumu2280

大神，这是我写的啊 {:2_171:} 回复 4# jzp520520


   

jzp520520

回复 5# xiumu2280 [/

碉堡了！我请教一下：if @hash_gene{@data};@hash_uniq{@{$hash_new{$key}}}=1;这两句不懂，可否解释下？

xiumu2280

这个是哈希切片的写法  回复 6# jzp520520


   

jzp520520

回复 7# xiumu2280

哈希切片没用过，也很少看到，我的理解是这样的，不知对不对，请大神赐教！
if @hash_gene{@data};
等效于
while (@data) {
        if (defined $hash_gene{$_} || exists $hash_gene{$_}) {
                BLOCK;
        }
}

@hash_uniq{@{$hash_new{$key}}}=1;
等效于
while (@{$hash_new{$key}}) {
        $hash_uniq{$_}=1;
}

rubyish

不知对不对?
2 类别
2, 3, 4, 5
7, 9, 10

rubyish

d e l ~ ~ ~