宝刀未老 | 写了一个 Ugly 的可能用于单细胞数据分析的脚本

写在前面

昨天傍晚，突然有个朋友联系到我，让帮忙写个 perl 脚本。想想是几年前认识的，后面也就基本没联系过了。那会我似乎是刚接手QQ社群“Bioinformatics中国”的管理。时间过得是真快。尽管确实没啥时间，我还是答应了，毕竟写个脚本嘛。

需求描述

有文本内容如下

0   40
1   6
1   11
1   12
1   17
1   19
1   20
1   31
1   38
2   4
2   5
2   6
2   10
2   31
2   38
4   6
4   10
4   31
4   38
5   6
5   38
6   12
6   17
6   20
6   31
6   38
8   18
9   13
9   23
9   30
12  20
15  26
17  19
17  20
17  25
17  38
19  20
19  25
19  31
19  38

其中每个数字按理是对应一个cluster，cluster和cluster有关联？或者overlap？就是一行，所以每行是两个数字。要求是把所有有之间关联或者间接关联的汇聚成一行。
Emmm，说的是perl脚本，但实际上我估计三年多来几乎没写过perl脚本。perl活在我的命令行one-liner世界。Anyway，我还是折腾折腾写了下。

脚本逻辑与具体实现

当时写了十来分钟，感觉快写完了，突然院楼就断电了。于是去打了会球，回来又折腾了大半个小时。原本是觉得挺简单的，不过还是花了一点时间。
整体实现逻辑简单

1. 读取整个文件，每两个有直接联系的cluster，存储起来
2. 遍历cluster，一旦发现可以合并的cluster，记录然后终止循环
3. 合并cluster，回到第二步，直到无法发现还能继续合并的cluster

于是得到脚本如下

#!/usr/bin/perl
use strict;
my $usage = "
    perl $0 in.file out.file
";
my $in = shift;
my $out = shift;
die $usage unless $in and $out;
open IN,'<',$in or die "Can't read input file $in\n";
open OUT,'>',$out or die "Can't write into file $out\n";
my @cluster;

while(<IN>){
    s/[\r\n]+//;
    push @cluster,[split /\s+/,$_];
}
close IN;
while(1){
    my $mergeCluster_1;
    my $mergeCluster_2;
    LOOP:for my $out(@cluster){
        my @out = @{$out};
        for my $in(@cluster){
            if($out == $in){
                next;
            }
            my @in = @{$in};
            for my $iter (@out){
                if((grep {$iter == $_} @in)){
                    $mergeCluster_1 = $out;
                    $mergeCluster_2 = $in;
                    # print "Hit\n";
                    last LOOP;
                }
            }
        }
    }
    if($mergeCluster_1){ # and $mergeCluster_2
        my @mergedCluster;
        my @updateCluster;
        for(@cluster){
            if($_!=$mergeCluster_1 && $_!=$mergeCluster_2){
                push @updateCluster,$_;
            }else{
                push @mergedCluster,@{$_};
            }
        }
        push @updateCluster,[@mergedCluster];
        @cluster = @updateCluster;
    }else{
        print "Finished...\n";
        last;
    }
}
for(@cluster){
    my %uniq;
    map {$uniq{$_}=1} @{$_};
    print OUT join qq{\t},sort {$a<=>$b} keys %uniq;
    print OUT "\n";
}

运行结果

0       40
8       18
15      26
9       13      23      30
1       2       4       5       6       10      11      12      17      19      20      25      31      38

Emmm，我还是觉得我挺厉害的。

写在最后

啊，宝刀未老。不过整个实现跟perl没啥关系，基本是C-style。
想想，写这些脚本，真是没意思。