perl novel可变剪接识别（1）

想把之前做的可变剪接模型给大家说一下，看看有什么遗漏的没有，由于当时想法比较复杂，所以程序有点多，大致分三个部分来进行。

首先，拿到的结果是tophat给出的junction的数据，其次博主使用的数据库是ensembl的数据库，gencode也可以，先得到已知的参考junction：

#!/usr/bin/env perl
use warnings;
use strict;

die "perl $0 <junction.bed> <ensembl|gencode> <bp>\n" unless @ARGV eq 3; 
&showtime("Start...");
my $in = shift;
my %gene;
open GTF, $in or die $!;
while(<GTF>){
	chomp;
	my @tmp = split;
	next if(/^#/);
	my ($gid) = $_ =~ /gene_name "([^;]+)";/;
	my ($tid) = $_ =~ /transcript_id "([^;]+)";/;
	if($tmp[2] =~ /exon/){
		push @{$gene{$gid}{$tid}{$tmp[0]}}, "$tmp[3]\t$tmp[4]";
	}
}
close GTF;

my %hash;
foreach my $g(keys %gene){
	foreach my $t(keys %{$gene{$g}}){
		foreach my $chr(keys %{$gene{$g}{$t}}){
			my $flag = 0;
			my $end;
			foreach my $str(sort @{$gene{$g}{$t}{$chr}}){
				if($flag == 0){
					$end = (split /\t/, $str)[1];
					$flag = 1;
				}else{
					my $f_end = (split /\t/, $str)[0];
					$hash{$chr}{$end}{$f_end} = 0;
					$end = (split /\t/, $str)[1];
				}
			}
		}
	}
}
&showtime("GTF is done..."); #上面两步将所有的EXON连接方式给存起来

my $junc = shift;
my ($id) = $junc =~ /^(\w+)_alignment/;
open OUT1, ">$id.novel" or die $!;
open OUT2, ">$id.ref" or die $!;
my $bp = shift; #设置精确范围，可以在参考junction的范围内都识别成已知的剪接，博主设置为0为了精确
$bp ||= 0;
open JUNC, $junc or die $!;
while(<JUNC>){
	chomp;
	next if($. == 1);
	my @tmp = split;
	my $flag = 0;
	my ($block_s, $block_e) = (split /,/, $tmp[10])[0,1];
	my $intron_s = $tmp[1] + $block_s;
	my $intron_e = $tmp[2] - $block_e + 1;
	foreach my $s(keys %{$hash{$tmp[0]}}){
		foreach my $e(keys %{$hash{$tmp[0]}{$s}}){
			if(abs($intron_s - $s) <= $bp and abs($intron_e - $e) <= $bp){
				$flag = 1;
				print OUT2 "$tmp[0]\t$tmp[1]\t$tmp[2]\t$tmp[4]\t$tmp[5]\t$block_s\t$block_e\n";
				delete $hash{$tmp[0]}{$s}{$e};
				last;
			}
		}
	}
	if($flag eq 1)
	{
		next;
	}else{
		print OUT1 "$tmp[0]\t$tmp[1]\t$tmp[2]\t$tmp[4]\t$tmp[5]\t$block_s\t$block_e\n";
	}

}
close JUNC;
&showtime("Junctions is done...");


sub showtime(){
	my ($str) = @_;
	my $time = localtime;
	print STDERR "$str\t$time\n";
}

这样就得到两部分数据，已知与新的——新的将进行接下来的处理工作。