1. 程式人生 > >perl novel可變剪接識別(1)

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";
}

這樣就得到兩部分資料,已知與新的——新的將進行接下來的處理工作。