本文主要是介绍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";
}这样就得到两部分数据，已知与新的——新的将进行接下来的处理工作。
这篇关于perl novel可变剪接识别（1）的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！
