在做 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}};
这里使用了 ？：的判断方式，问号前面代表判断的条件，问号后面冒号前面代表条件判断如果为真则执行，否则就执行冒号后面的语句。
欢迎大家留言讨论！
