生信必备技能 -- 不可描述

在做 Small RNA 分析的过程中时遇到了一个有意思的事情，前期比对到不同的数据库后需要将不同的短 tags 注释到不同的类别上，这个地方有点意思的事情是同一条 tags 可能被注释到不同的类别上，这里就需要对不同的类别进行优先级排序了。
一般按照惯例我们会根据 已知 miRNA -> rRNA -> tRNA -> ncRNA (snRNA, snoRNA) 的顺序进行优先级注释。

在写程序的时候基本上有以下的思路：
首先我获得了每个数据库的比对结果，知道了哪些 tags 比对到了相应的数据库中，使用最基本的思路是从 总的 tags 中先剔除比对到相应物种的成熟 miRNA，然后剔除比对到物种相应的总的 hairpin ( miRNA 前体数据库)，然后剔除比对到 GeneBank 数据库中的 rRNA 以及 tRNA ， 如果有 snRNA 或者 snoRNA 等 tags 的话依次提取出来。

思路很简单，但是真正动手写脚本的时候才发现，这个过程苏是多么的冗余。因为，这样思路下我们如果要比对7个数据库的话，最少就应该用同样的脚本过滤六次。好麻烦，还好我比较喜欢偷懒，写脚本的时候应用了另一种思路：
直接处理所有的 tags 然后对不同的数据库比对结果信息分别建立一个文件，通过不同的文件标签进行 tags 分类标记，这样就省事多了：
附上源代码：( 遵循 GLP v3 开源协议 )

#!/usr/bin/perl -w
use strict;

die "Usage: perl $0 <list> <prefix> <indir> \n" if @ARGV < 3;

my $info_dir = $ARGV[2];
my $sample = $ARGV[1];

my ( %list, %out );

# deal with all tags
open IN,$ARGV[0] or die $!;
while (<IN>){
    chomp;
    next if /^#/;
    my ($id, $cnt) = split /\t/,$_;
    $list{$id}{"count"} = $cnt;
}
close IN;

# indir must contain files like: $sample.m2mature.info 
# t0023555        xxx_tRNA_AM087200_6:80:-        1       20      +
# t0031916        xxx_tRNA_AM087200_6:80:-        1       21      +

opendir DIR, $info_dir or die $!;
foreach my $file (sort grep(/$sample\.\S+\.info$/,readdir(DIR))){
    (my $type) = $file =~ /$sample\.(\S+)\.info/;
    if ($type eq "m2ncgb"){
        open IN,$file or die $!;
        while (<IN>)  {
            chomp;
            next if /^#/;
            next if /^$/;
            my @b = split /\t/,$_;
            my @c = split /_/,$b[1];
            push @{$out{$b[0]}},$c[1];
        }
        close IN;
    } else {
        open IN,$file or die $!;
        while (<IN>)  {
            chomp;
            next if /^#/;
            next if /^$/;
            my @b = split /\t/,$_;
            push @{$out{$b[0]}},$type;
        }
        close IN;
    }
}
close DIR;

foreach my $k (sort keys %list){
    if (exists $out{$k}){
        my $t = @{$out{$k}}==1 ? $out{$k}->[0] :join ",",@{$out{$k}};
        print "$k\t$t\n";
    }else {
        print "$k\t-\n";
    }
}

附上大牛的代码：行使同样的功能，代码量只有我的一半！

#! /usr/bin/perl -w
use strict;
die "perl $0  input.txt prefix  dir\n"  unless @ARGV == 3;

chomp(my @file = glob "$ARGV[2]/$ARGV[1].*info");
my %ha;
foreach(@file){
   my $flag = (split/\./,$_)[-2];
   ding($_,\%ha,$flag);
}
open IN,$ARGV[0]||die;
while(<IN>){
    chomp;
    my $id = (split/\t/,$_)[0];
    $ha{$id} ||= "NA,";
    chop $ha{$id};
    print "$id\t$ha{$id}\n";
}
close IN;

############################################
sub ding{
    my ($p,$q,$f) = @_;
    open IN,$p||die;
    while(<IN>){
        chomp;
        next if /^#|^$/;
        my $a = (split/\t/,$_)[0];
        $q->{$a} .="$f,";
    }
}

perl 真是一门强大的脚本们语言。同样的功能有多种多样的方法可以将它演绎出来！

** 这里值得重点划线的地方是 这一句代码：**

my $t = @{$out{$k}}==1 ? $out{$k}->[0] :join ",",@{$out{$k}};

这里使用了 ？：的判断方式，问号前面代表判断的条件，问号后面冒号前面代表条件判断如果为真则执行，否则就执行冒号后面的语句。

欢迎大家留言讨论！