在做 Small RNA 分析的過程中時(shí)遇到了一個(gè)有意思的事情，前期比對(duì)到不同的數(shù)據(jù)庫后需要將不同的短 tags 注釋到不同的類別上锹安，這個(gè)地方有點(diǎn)意思的事情是同一條 tags 可能被注釋到不同的類別上煮嫌，這里就需要對(duì)不同的類別進(jìn)行優(yōu)先級(jí)排序了器净。
一般按照慣例我們會(huì)根據(jù) 已知 miRNA -> rRNA -> tRNA -> ncRNA (snRNA, snoRNA) 的順序進(jìn)行優(yōu)先級(jí)注釋尖淘。
在寫程序的時(shí)候基本上有以下的思路：
首先我獲得了每個(gè)數(shù)據(jù)庫的比對(duì)結(jié)果，知道了哪些 tags 比對(duì)到了相應(yīng)的數(shù)據(jù)庫中倡勇，使用最基本的思路是從 總的 tags 中先剔除比對(duì)到相應(yīng)物種的成熟 miRNA逞刷，然后剔除比對(duì)到物種相應(yīng)的總的 hairpin ( miRNA 前體數(shù)據(jù)庫)，然后剔除比對(duì)到 GeneBank 數(shù)據(jù)庫中的 rRNA 以及 tRNA 妻熊， 如果有 snRNA 或者 snoRNA 等 tags 的話依次提取出來夸浅。
思路很簡(jiǎn)單，但是真正動(dòng)手寫腳本的時(shí)候才發(fā)現(xiàn)扔役，這個(gè)過程蘇是多么的冗余帆喇。因?yàn)椋@樣思路下我們?nèi)绻葘?duì)7個(gè)數(shù)據(jù)庫的話亿胸，最少就應(yīng)該用同樣的腳本過濾六次坯钦。好麻煩预皇，還好我比較喜歡偷懶，寫腳本的時(shí)候應(yīng)用了另一種思路：
直接處理所有的 tags 然后對(duì)不同的數(shù)據(jù)庫比對(duì)結(jié)果信息分別建立一個(gè)文件婉刀，通過不同的文件標(biāo)簽進(jìn)行 tags 分類標(biāo)記吟温，這樣就省事多了：
附上源代碼：( 遵循 GLP v3 開源協(xié)議 )
#!/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 真是一門強(qiáng)大的腳本們語言鲁豪。同樣的功能有多種多樣的方法可以將它演繹出來！
** 這里值得重點(diǎn)劃線的地方是 這一句代碼：**
my $t = @{$out{$k}}==1 ? $out{$k}->[0] :join ",",@{$out{$k}};
這里使用了 律秃？：的判斷方式爬橡，問號(hào)前面代表判斷的條件，問號(hào)后面冒號(hào)前面代表?xiàng)l件判斷如果為真則執(zhí)行棒动，否則就執(zhí)行冒號(hào)后面的語句糙申。
歡迎大家留言討論！
